Early Results from GLASS-JWST. XXIV. The Mass–Metallicity Relation in Lensed Field Galaxies at Cosmic Noon with NIRISS^*

We utilize the Grism Redshift and Line Analysis software Grizli (Brammer 2023) to reduce NIRISS data using the standard JWST pipeline (version 1.11.1) and the latest reference file (under the jwst_1100.pmap context). The detailed procedures are largely described in Roberts-Borsani et al. (2022). Briefly, Grizli analyzes the paired direct imaging and grism exposures through forward modeling and yields contamination-subtracted 1D and 2D grism spectra, along with the best-fit spectroscopic redshifts.

For each source, the 1D spectrum is constructed using a linear superposition of a spectra from a library consisting of four sets of empirical continuum spectra covering a range of stellar population ages (Brammer et al. 2008; Erb et al. 2010; Conroy & van Dokkum 2012; Muzzin et al. 2013) and Gaussian-shaped nebular emission lines at the observed wavelengths given by the source redshift. The intrinsic 1D spectrum and the spatial distribution of flux measured in the paired direct image are utilized to generate a 2D model spectrum based on the grism sensitivity and dispersion function, similar to the "fluxcube" model produced by the aXe software (Walsh et al. 2009). This 2D forward-modeled spectrum is then compared to the observation by Grizli and a global χ² calculation is performed to determine the best-fit superposition coefficients for both the continuum templates and Gaussian amplitudes, the latter of which correspond to the best-fit emission line fluxes. In this way, our 2D forward-modeling practice not only determines the source redshift, but also measures the emission line fluxes, taking into account the morphological broadening effect. We refer the interested readers to Appendix A of Wang et al. (2019), for the full descriptions of the redshift fitting procedure.

We obtain a parent sample of 4756 sources with F150W apparent magnitudes between [18, 32] ABmag (the 5σ depth is 28.7 according to Treu et al. 2022), on which our Grizli analyses result in meaningful redshift constraints. Several goodness-of-fit criteria are implemented to ensure the reliability of our redshift fit: a reduced chi-square close to 1 (χ² < 2.2), a sharply peaked posterior of the redshift ${\left({\rm{\Delta }}z\right)}_{\mathrm{posterior}}/(1+{z}_{\mathrm{peak}})\lt 0.002$, high evidence of Bayesian information criterion compared to polynomials (BIC > 100). As a result, there are 348 sources in the redshift range z ∈ [0.05, 10], with secure grism redshift measurements according to the above joint selection criteria. A total of 86 sources with secure grism redshifts are at redshifts z ∈ [1.8, 2.3] ∪ [2.6, 3.4], ensuring that the slitless spectra cover several emission lines: [O ii], [Ne iii], Hδ, Hγ, Hβ, and [O iii] (also Hα, [S ii] for the former zone), with high sensitivity for our three NIRISS filters (F115W, F150W, and F200W). However, 6/86 sources of our NIRISS spectroscopy catalog do not match entries in the NIRCam photometric catalog (Paris et al. 2023) within 07 (5 × PSF).

The fluxes of the intrinsic nebular emission lines ([O ii], [Ne iii], Hδ, Hγ, Hβ, [O iii], Hα, and [S ii], the same as in Henry et al. 2021) are 2D forward modeled by Grizli as output. There are 57 sources with Hβ detection, to ensure the reliable measurement of SFR. No other emission line criteria (e.g., S/N [O iii]) are used for selection, to avoid potential metallicity bias. Then we visually inspect the 1D spectra of each galaxy individually, excluding seven of those that are heavily contaminated. The 50 galaxies showing prominent nebular emission features, with zero possible active galactic nucleus (AGN) exclusions in Section 3.4, will make up the final sample presented in Table A1. A "textbook case" of our samples (ID: 05184 in Table A1) has been carefully studied through spatial mapping in our recent work (Wang et al. 2022a). We show as an example 1D/2D spectra for six galaxies in our sample in Figure 1, annotated with their exposure times, best-fit grism redshifts, and stellar masses (which will be discussed in Section 3.3).

Since the 1D grism spectra are extracted by Grizli simultaneously, it allows us to directly fit it using several 1D Gaussian profiles to obtain line fluxes and errors, as detailed in Section 3.5. But we still use the previous 2D flux other than 1D as our default result for subsequent calculations. The comparison of the line flux measurements between this 1D line profile fitting and the 2D Grizli forward-modeling procedure, is discussed in Section 4.2.

We use these observed line fluxes $({f}_{i}^{{\rm{o}}},{\sigma }_{i}^{{\rm{o}}})$ to simultaneously estimate three parameters: jointly metallicity, nebular dust extinction, and dereddened Hβ line flux ($12+\mathrm{log}({\rm{O}}/{\rm{H}})$, A_v , f_Hβ ). We follow the previous series of work (Jones et al. 2015b; Wang et al. 2017, 2019, 2020, 2022b) by constructing a Bayesian inference method that uses multiple calibration relations to jointly constrain metallicity $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ and (A_v , f_Hβ ) simultaneously. Our method is more reliable than the conventional way of turning line flux ratios into metallicities, since it takes into account the intrinsic scatter in strong line O/H calibrations (${\sigma }_{{R}_{i}}$ in Equation (1)). And it combines multiple line flux measurements and properly marginalizes over the dust extinction correction. It also emphasizes bright lines (e.g., [O ii], [O iii]) with high signal-to-noise ratios (S/Ns) and marginalizes faint lines (e.g., Hβ) or even nondetection lines with low S/Ns quantitatively, (i.e., by assigning weights to each line according to its S/N in the likelihood function).

The Markov Chain Monte Carlo (MCMC) sampler Emcee software (Foreman-Mackey et al. 2013) is employed to sample the likelihood profile ${ \mathcal L }\propto \exp (-{\chi }^{2}/2)$ with

$$\begin{eqnarray}&&{\chi }^{2}:= \displaystyle \sum _{i}^{\mathrm{EL}}\displaystyle \frac{{\left({f}_{i}-{R}_{i}\cdot {f}_{{\rm{H}}\beta }\right)}^{2}}{{\left({\sigma }_{i}\right)}^{2}+{\left({\sigma }_{{R}_{i}}\right)}^{2}\cdot {f}_{{\rm{H}}\beta }^{2}},\quad {R}_{i}:= \displaystyle \frac{{f}_{i}}{{f}_{H\beta }}.\end{eqnarray} \tag{ 1 }$$

Here the summation i includes all emission lines, with their intrinsic scatters ${\sigma }_{{R}_{i}}:= {\sigma }_{i}^{\mathrm{cal}}\cdot {R}_{i}\cdot \mathrm{ln}10$. The inherent flux and uncertainty (f_i , σ_i ) for each line are corrected from observation $({f}_{i}^{{\rm{o}}},{\sigma }_{i}^{{\rm{o}}})$ for dust attenuation by parameter A_v using the Calzetti et al. (2000) extinction law. R_i refers to the line flux ratio, which is empirically calibrated by a polynomial as a function of metallicity: $\mathrm{log}R={\sum }_{j=0}^{n}{c}_{j}\cdot {(x)}^{j},x:= 12+\mathrm{log}({\rm{O}}/{\rm{H}})$, where (x)^j means jth power of x, with the coefficients summarized in Table 1. For flux ratio calibrations that do not use Hβ as the denominator (e.g., [Ne iii]/[O iii]), the terms f_Hβ in Equation (1) need to be replaced by the corresponding lines (e.g., f_{[O iii]}). And one more term of uncertainty (e.g., ${\sigma }_{{\rm{O}}3}^{2}\cdot {R}_{\mathrm{Ne}3}^{2}$) needs to be added to the denominator of χ².

Table 1. Coefficients for the Emission Line Flux Ratio Diagnostics Used in This Work

Diagnostic R and Notation	c0	c1	c2	c3	c4	c5	σcal	ref
O3:= [O iii]/Hβ	43.9836	−21.6211	3.4277	−0.1747	⋯	⋯	0.05	(Bian et al. 2018, B18)
O2:= [O ii]/Hβ	78.9068	−45.2533	7.4311	−0.3758	⋯	⋯	0.05
Hα/Hβ	0.45637	⋯	⋯	⋯	⋯	⋯	0.00	Balmer decrement
Hγ/Hβ	−0.32790	⋯	⋯	⋯	⋯	⋯	0.00
Ne3O3:= [Ne iii]/[O iii]	−1.11420	⋯	⋯	⋯	⋯	⋯	0.04	Jones et al. (2015a)
	−0.54571	0.45730	−0.82269	-0.02839	0.59396	0.34258	best
S2:= [S ii]/Hα	−0.43974	0.34034	−0.62850	-0.07077	0.47147	0.31767	upper	Jones et al. (2015a) b
	−0.65464	0.58976	−1.06047	0.01979	0.75382	0.37766	lower

Notes.

^a We note that the [O iii]/Hβ calibration reported in Bian et al. (2018) in fact refers to the flux ratio between [O iii] 4960,5008 and Hβ, i.e., a factor of (2.98/3.98) is needed (following Storey & Zeippen 2000) when we use the doublets let to calibrate pure [O iii] 5008. ^b The line flux ratio R_{[S ii ]} is calibrated by polynomial with coefficients given by the "best" row, and the uncertainty σ_{[S ii]} is given by the "upper" and "lower" rows, where the metallicity x is relative to solar $x:= 12+\mathrm{log}({\rm{O}}/{\rm{H}})-8.69$.

Download table as:  ASCIITypeset image

A wide range of strong line calibrations between line flux ratio and metallicity has been established (see Appendix C in Wang et al. 2019 for a summary; also see Kewley et al. 2019; Maiolino & Mannucci 2019 for recent reviews). Different choices can result in offsets as high as 0.7 dex (see, e.g., Kewley & Ellison 2008). In this work, we adopt mainly the diagnostics group "O₃–O₂" of calibrations prescribed by Bian et al. (2018, hereafter B18), for comparison with Sanders et al. (2021) and Wang et al. (2022b). The purely empirical calibrations in Bian et al. (2018, B18) are based on a sample of local analogs of high-z galaxies according to the location on the Baldwin, Phillips & Terlevich (BPT) diagram (Baldwin et al. 1981), with the notations and coefficients summarized in Table 1.

These calibrations are recommended for the metallicity range of $7.8\lt 12+\mathrm{log}({\rm{O}}/{\rm{H}})$, which is appropriate for our sample that does not reach metallicities as low as those found at higher redshift (Curti et al. 2023b; Heintz et al. 2023). As a sanity check, we computed metallicities using the calibrations from Sanders et al. (2023), and indeed we do not find galaxies with metallicities significantly lower than 7.8. In order to make complete use of emission lines of spectra, we also collect Ne₃O₃, S₂ diagnostics at the same time, even though the corresponding line fluxes are not so strong for our sample. We have tested that if they are removed, they do not significantly affect the metallicity estimation, which is dominated by the first two diagnostics O₃, O₂ in B18 and two Balmer decrements. We adopt the intrinsic Balmer decrement flux ratios assuming case B recombination with T_e ∼ 10,000 K. We neglect the line-blending effect, since they are likely small in most cases (see Figure 4 and Appendix C in Henry et al. 2021 for more information). This Bayesian method is used to derive properties ($12+\mathrm{log}({\rm{O}}/{\rm{H}})$, A_v , f_Hβ ) of galaxies both from our individual spectra sample here and from the stacked spectra presented in Section 3.5.

From the dereddened Hβ flux f_Hβ , we estimate the instantaneous SFR of our sample galaxies, based on Balmer line luminosities. This approach provides a valuable proxy of the ongoing star formation on a timescale of ∼10 Myr, highly relevant for galaxies displaying strong nebular emission lines. Assuming the Kennicutt (1998) calibration and the Balmer decrement ratio of Hα/Hβ = 2.86 from the case B recombination for typical H ii regions, we calculate

$$\begin{eqnarray}&&\mathrm{SFR}=4.65\times {10}^{-42}\displaystyle \frac{L({\rm{H}}\beta )}{[\mathrm{erg}\,{{\rm{s}}}^{-1}]}\times 2.86\,[{M}_{\odot }\,{\mathrm{yr}}^{-1}],\end{eqnarray} \tag{ 2 }$$

suitable for the Chabrier (2003) initial mass function. The total luminosity $L({\rm{H}}\beta )=4\pi {D}_{L}^{2}(z)\cdot {f}_{{\rm{H}}\beta }$ is corrected for lensing magnification according to Bergamini et al. (2023a). The corrected SFR values are given in Table A1.

In this section, we fit broadband photometry to obtain stellar mass M_* of target galaxies through SED fitting. We directly use the combined photometric catalog released by the GLASS-JWST team (Paris et al. 2023). The photometric fluxes measured within 2 × PSF FWHM apertures of all 16 bands are included if available. We match 2983/4756 galaxies of our NIRISS spectroscopy catalog in Section 3.1 to the 24,389 galaxies of the NIRCam photometric catalog with on-sky distances (d2d) lower than 07 (5 × FWHM in the F444W band, conservatively). As done in Section 3.1, the final selected sample of 50 galaxies yields accurate d2d match (<014, around the angular resolution of JWST/NIRISS), and visually crossmatching with the NIRCam image further validates our sources.

To estimate the stellar masses M_* of our sample galaxies, we use the Bagpipes software (Carnall et al. 2018) to fit the BC03 (Bruzual & Charlot 2003) models of SEDs to the photometric measurements derived above. We assume the Chabrier (2003) initial mass function, a metallicity range of Z/Z_⊙ ∈ (0, 2.5), and the Calzetti et al. (2000) extinction law with A_v in the range of (0, 3). We use the double power-law (DPL) model rather than the simple exponentially declining form to capture the complex star formation history (SFH) of our galaxies at Cosmic Noon (rather than local Universe), following Carnall et al. (2019). The nebular emission component is also added into the SED during the fit, since our galaxies are exclusively strong line emitters by selection. The redshifts of our galaxies are fixed to their best-fit grism values, with a conservative uncertainty of z_σ = 0.003. Note that we have obtained the entire redshift posterior from Grizli in Section 3.1 and set a criterion of ${\left({\rm{\Delta }}z\right)}_{\mathrm{posterior}}/(1+{z}_{\mathrm{peak}})\lt 0.002$ for secure redshift measurements. But here we still set a Gaussian prior centered on z_peak with z_σ = 0.003 for simplicity in SED fitting, following Momcheva et al. (2016). Actually, the minimum, median, and maximum values of Δz/(1 + z) for our sample are 1.4 × 10⁻⁴, 2.8 × 10⁻⁴, and 1.5 × 10⁻³, respectively.

Our mass estimates are in agreement with Santini et al. (2023), even though we stress that our results are more robust, because we use spectroscopic redshifts. After correcting magnification according to our recent lensing model (Bergamini et al. 2023a), we are allowed to take a glimpse of the loci of our galaxies in the SFR–M_* diagram as in Figure 2. We show the star-forming main sequence fitted by Speagle et al. (2014), which is extrapolated from $\mathrm{log}({M}_{* }/{M}_{\odot })\in [9.7,11.1]$ to the mass range of our sample with ±0.2 dex scatters. Sanders et al. (2021) give stacked results of field galaxies fairly close to their extrapolated best fit out to $\mathrm{log}({M}_{* }/{M}_{\odot })=9$. Our sample generally scatters around the main sequence at higher M_*. But at lower M_* high-SFR galaxies are dominant, especially for z ∼ 3 at M_*/M_⊙ ≲ 3 × 10⁸. It might account for the low metallicity at the low-mass region when assuming the FMR (Mannucci et al. 2010), which will be discussed in Section 4.1.

Zoom In Zoom Out Reset image size

The metallicity diagnostics used in this work are strictly for star-forming regions/galaxies, and the results will be incorrect if there is AGN emission. So the last step is to exclude the AGN contamination from purely star-forming galaxies, by using the mass-excitation (MEx) diagram as shown in Figure 3.

Zoom In Zoom Out Reset image size

AGNs leave strong signatures on nebular line ratios such as O III λ5007/Hβ and/or N II λ6584/Hα, which form the most traditional version of the BPT diagram (Baldwin et al. 1981). Due to the limited spectral resolution of JWST/NIRISS slitless spectroscopy (R ∼ 150), [N ii] is entirely blended with Hα, which precludes us from using the BPT diagram to remove AGN contamination.

Fortunately, Juneau et al. (2014) proposed an effective approach coined the MEx diagram, using M_* as a proxy for [N ii]/Hα, which functions well at z ∼ 0 (i.e., SDSS DR7). Coil et al. (2015) further modified the MEx demarcation by horizontally shifting these curves to high M_* by 0.75 dex, which is shown to be more applicable to the MOSDEF sample (Sanders et al. 2021) at z ∼ 2.3. We thus rely on this modified MEx to prune AGN contamination from our galaxy sample. As shown in Figure 3, the green and red curves mark the steep gradient of P(AGN) ∼ 0.3 and P(AGN) ∼ 0.8, respectively, which represent the probability that the galaxy hosts an AGN.

Most of the sources are clearly unlikely AGN, and some scattered around the critical line are ambiguous. There are only two galaxies slightly above the upper demarcation within 1σ. Because our analysis is based on stacking, a small minority of contaminating AGN will have a negligible impact. Given the limited sample size, we tend to retain more applicable data, and consequently, no possible AGN is eliminated and we preserve all 50 galaxies.

Robust emission lines are required to estimate metallicity for MZR measurement. So we need composite spectra obtained by stacking procedure to achieve higher S/N from low-resolution grism spectra. In the previous subsection, we have selected 50 spectroscopically confirmed galaxies in the A2744-lensed field that are undergoing active star formation. Then they are divided into two redshift bins (z ∈ [1.8, 2.3] and z ∈ [2.6, 3.4]), and three mass bins, respectively, as in Table 2. Our choice of binning aims to have a reasonable number of galaxies per bin. We tested that changing the mass bins does not significantly affect our conclusions. Approximately each mass bin contains ∼seven individual galaxies, and the S/N will be increased roughly by a factor of $\sqrt{7}=2.6$. The 1D/2D spectra of representative galaxies in each of the six bins are shown in Figure 1.

Table 2. Measured Properties of the Stacked Spectra

Group	Ngal	Mass Range	log${M}_{* }^{\mathrm{med}}$	[O iii]/Hβ	[O ii]/Hβ	[O iii]/[O ii]	Hγ/Hβ	[Ne iii]/[O iii]	Hα/Hβ	[S ii]/Hα	$12+\mathrm{log}({\rm{O}}/{\rm{H}})$
1.8 < zgrism < 2.3
11	7	[6.8,7.7)	7.42	7.21 ± 0.55	0.75 ± 0.07	9.62 ± 0.58	0.26 ± 0.04	0.07 ± 0.01	2.88 ± 0.26	0.01 ± 0.02	${8.00}_{-0.04}^{+0.05}$
12	10	[7.7,8.7)	8.20	7.27 ± 0.63	1.40 ± 0.16	5.20 ± 0.44	0.13 ± 0.07	0.01 ± 0.02	2.77 ± 0.29	0.08 ± 0.03	${8.15}_{-0.05}^{+0.05}$
13	11	[8.7,9.9)	9.09	4.84 ± 0.23	1.93 ± 0.15	2.51 ± 0.17	0.09 ± 0.05	0.02 ± 0.02	3.56 ± 0.24	0.09 ± 0.03	${8.37}_{-0.05}^{+0.04}$
2.6 < zgrism < 3.4
21	5	[7.1,8.2)	7.84	5.10 ± 0.64	0.32 ± 0.07	15.83 ± 2.94	...	0.04 ± 0.01	...	...	${7.98}_{-0.12}^{+0.19}$
22	9	[8.2,9.2)	8.84	7.48 ± 0.55	1.75 ± 0.15	4.29 ± 0.23	0.23 ± 0.15	0.03 ± 0.01	...	...	${8.25}_{-0.06}^{+0.05}$
23	8	[9.2,10.0)	9.57	3.91 ± 0.33	1.80 ± 0.24	2.17 ± 0.22	0.28 ± 0.16	0.05 ± 0.04	...	...	${8.47}_{-0.06}^{+0.05}$

Notes.

The multiple emission line flux ratios are measured from the stacked spectra shown in Figure 4. The mass range and the median stellar mass $\mathrm{log}{M}_{* }^{\mathrm{med}}$ are both logarithmic values $\mathrm{log}({M}_{* }/{M}_{\odot })$. The metallicity inference is derived from the measured line flux ratios in the stacked spectra presented in each corresponding row, using the method described in Section 3.2. Here we use the strong line calibrations prescribed by Bian et al. (2018, B18) and some others. See Table 1 for the relevant coefficients.

Download table as:  ASCIITypeset image

Then we adopt the following stacking procedures, similar to those utilized by Henry et al. (2021) and Wang et al. (2022b):

• 1.  
  Subtract continuum models from the extracted grism spectra. The continua are constructed by Grizli combining two orients. We apply a multiplicative factor to the continuum models to make sure there is no offset between the modeled and observed continuum levels around emission lines, to avoid continuum oversubtraction.
• 2.  
  Normalize the continuum-subtracted spectrum of each object using its measured [O iii] flux, to avoid excessive weighting toward objects with stronger line fluxes. Here the [O iii] fluxes we used are the results of 1D line profile fitting instead of 2D forward modeling by Grizli, for a more straightforward normalization.
• 3.  
  De-redshift each normalized spectrum to its rest frame, and resample on the same wavelength grid using SpectRes ²⁶ with the integrated flux preservation.
• 4.  
  Take the median and the variance of the normalized fluxes at each wavelength grid as the value and uncertainty of the stacked spectrum.

As shown in Figure 4, these key lines are more significant in stacked spectra. The (relative) emission line fluxes are measured by fitting a set of Gaussian profiles to the line in stacked spectra, as well as individual spectra. We simultaneously fit [O ii], [Ne iii], Hδ, Hγ, Hβ, [O iii], Hα, and [S ii]. The amplitude ratio of [O iii] λ λ4960, 5008 doublets is fixed to 1:2.98 following Storey & Zeippen (2000). The centroids of Gaussian profiles are allowed a small shift of the corresponding rest-frame wavelengths of emission lines, within ±10 Å, in order to accommodate systematic uncertainties. The FWHMs of each line are not required to be the same, but are set between [10, 25]Å, consistent with the rest-frame spectral resolution Δλ ≈ 7 Å corresponding to R ≈ 150 for NIRISS. We use the software LMFit ²⁷ to perform the nonlinear least-squares minimization, with the measured quantities summarized in Table 2. The stacked metallicity is estimated using the same methods as the individual galaxies outlined in Section 3.2. Our later discussion will mainly focus on the stacked results.

Zoom In Zoom Out Reset image size

Our key scientific result is the measurement of the gas-phase MZR in the low mass range of $\mathrm{log}({M}_{* }/{M}_{\odot })\in (6.9,10.0)$ at z ∈ (1.8, 3.4). The slope of the MZR has been shown to be a key diagnostic of galaxy chemical evolution and the cycling of baryons and metals through star formation and gas flows (see e.g., Maiolino & Mannucci 2019, and references therein). In particular, Sanders et al. (2021) argue that the shape of the MZR at z ∼ 2–3 is more tightly regulated by the efficiency of metal removal by gas outflows ζ_out, rather than by the change of gas fractions with stellar mass μ_gas(M_*). Henry et al. (2013a) observe a steepening of the MZR slope at z ∼ 2, suggesting a transition from momentum-driven winds to energy-driven winds as the primary prescription for galactic outflows in the low-mass end.

We find a clear correlation between metallicity and stellar mass for both individual galaxies and stacked spectra at z ∈ [1.8, 2.3] and z ∈ [2.6, 3.4], as shown in the left panel of Figure 5. The z ∼ 2 and z ∼ 3 individual galaxy samples have Spearman correlation coefficients of 0.788 and 0.688 with p-values of 6.36 × 10⁻⁷ and 3.98 × 10⁻⁴, respectively. We perform linear regression over the stacks to derive the MZR:

$$\begin{eqnarray}&&12+\mathrm{log}({\rm{O}}/{\rm{H}})=\beta \times \mathrm{log}({M}_{* }/{10}^{8}{M}_{\odot })+{Z}_{8},\end{eqnarray} \tag{ 3 }$$

where β is the slope and Z₈ is the normalization at M_* = 10⁸ M_⊙, as the blue and red solid lines with uncertainties at z ∼ 2 and 3 in both panels of Figure 5. We measure the MZR slope to be β = 0.223 ± 0.017 and β = 0.294 ± 0.010 for our galaxy samples at z_median = 1.90 and z_median = 2.88, respectively. We see moderate evolution in the MZR normalization from z ∼ 2 to z ∼ 3: ΔZ₈ = − 0.11 ± 0.02. The stacked MZRs demonstrate good agreement with the individual results (linear fits are shown in the shaded regions in the left panel of Figure 5). The large uncertainty of the stacked metallicity in the z ∼ 3 lowest-mass bin, comes from the limited number of galaxies. More importantly, all five galaxies within this bin are high-SFR galaxies (Figure 2), which might explain their low stacked metallicity, under the assumption that the star-forming main sequence (Speagle et al. 2014) and the FMR (Mannucci et al. 2010) are valid below M_* ≲ 8. A detailed study and characterization of incompleteness at the low-mass end is beyond the scope of this paper and is left for future work.

Zoom In Zoom Out Reset image size

We summarize our measurements in Table 3, along with other literature results. The right panel of Figure 5 shows the comparisons to other observations and two cosmological hydrodynamic simulations. In addition to z ∼ 2 and 3, we also include the three latest MZR measurements at a very high redshift from JWST/NIRSpec for comparison. We measure the slope of the MZR to be β ∼ 0.25 for both z ∼ 2 and z ∼ 3. Our slopes at low mass are slightly lower than those found by Sanders et al. (2021), but ours are in lower mass ranges. The shallower normalization could be accounted for the MZR evolution from ours z_median = 1.90 and 2.88 to theirs z ∼ 2.3 and 3.3. Furthermore, we follow their analytical model to understand what physical processes set the slope at the dwarf mass range. In the Peeples & Shankar (2011) model, the metallicity of the ISM is expressed as

$$\begin{eqnarray}&&{Z}_{\mathrm{ISM}}=\displaystyle \frac{y}{{\zeta }_{\mathrm{out}}-{\zeta }_{\mathrm{in}}+\alpha {\mu }_{\mathrm{gas}}+1}.\end{eqnarray} \tag{ 4 }$$

Following the assumption by Sanders et al. (2021) that the gas fraction ${\mu }_{\mathrm{gas}}={10}^{{\mu }_{0}}{{M}_{* }}^{-0.36}$ (μ₀ = 3.89 and 3.96 for z ∼ 2 and 3, respectively), the coefficient α = 0.7 · (0.64 + β), the nucleosynthetic stellar yield $y/{Z}_{\mathrm{ISM}}={10}^{9.2-(12+\mathrm{log}({\rm{O}}/{\rm{H}}))}$, and the metal loading factors of inflowing gas accretion ζ_in = 0, we calculate the loading factors of outflowing galactic winds ζ_out at each stacked point and linear fit. We get

$$\begin{eqnarray}&&\begin{array}{l}z\sim 2:\mathrm{log}({\zeta }_{\mathrm{out}})=(-0.130\pm 0.072){m}_{10}+(0.408\pm 0.119),\\ z\sim 3:\mathrm{log}({\zeta }_{\mathrm{out}})=(-0.332\pm 0.037){m}_{10}+(0.202\pm 0.035),\\ \mathrm{where},\quad {m}_{10}=\mathrm{log}({M}_{* }/{10}^{10}{M}_{\odot }).\end{array}\end{eqnarray} \tag{ 5 }$$

And we find that $\mathrm{log}({\zeta }_{\mathrm{out}}/\alpha {\mu }_{\mathrm{gas}})$ is only a little bit above zero over the mass range, with ζ_out ≈ 1.01–1.5 × α μ_gas. Thus, our results indicate that the shallower MZR may be attributed to a shallower M_* scaling of the metal loading of the galactic outflows ζ_out at the low-mass end. We generalize their conclusions that outflows ζ_out remain the dominant mechanism other than gas fraction μ_gas that sets the MZR slope, and μ_gas gradually carries more relative importance and rise to nearly the same order as ζ_out for the low-mass regime.

Our MZR slope β ∼ 0.25 is steeper than those reported in Li et al. (2023) at the same redshifts and similar mass range as in Table 3. Although we use the same NIRISS data of the A2744-lensed field, we only match 28 out of 50 galaxies with the on-sky distances (d2d) lower than 1'' to the Abell catalog of Li et al. (2023), and only 18/50 of them are in agreement with our metallicity measurements within 1σ confidence interval. This difference likely arises from the updated calibration files used in our NIRISS data reduction, and from our Bayesian approach in the metallicity inference using multiple line ratios to joint fit other than only [O iii]/[O ii] from Bian et al. (2018). In addition, we include the new JWST/NIRCam imaging data covering the rest-frame optical wavelength ranges for our sample galaxies (Paris et al. 2023), use more complex SFH (DPL), and employ the latest JWST-based lensing model (Bergamini et al. 2023a) for more reliable stellar mass estimates. Another source of difference is their choice of exponentially declining SFH (τ model), which may not be appropriate for our high-redshift star-forming galaxies (Reddy et al. 2012), and might introduce a significant bias in stellar mass M_* estimation (Pacifici et al. 2015; Carnall et al. 2018, 2019).

In agreement with previous work, we also find a tendency for the slope of the MZR to flatten out in the low mass at around M_*/M_⊙ ≲ 10⁹, although not as significant. As for higher redshift z ∼ 3–10, our inferred slopes β are consistent with those by Curti et al. (2023b) and Nakajima et al. (2023), but our intercept Z₈ are ∼0.3 dex higher. At that time, the metal might be enriching and hence the MZR might be building up (Curti et al. 2023b), and it is not until the SFR peaks at Cosmic Noon z ∼ 2–3 that the MZR exhibits a higher intercept.

The MZR measurements are also sensitive to different strong line calibrations, especially for the intercept Z₈ (Kewley & Ellison 2008), as discussed in Section 3.2. In Table 3, we also provide the MZR from stacks using the Sanders et al. (2023) calibration for comparison. Although the measured slopes are significantly steeper than our default B18 MZR, they are still consistent with Heintz et al. (2023) for dwarf galaxies at higher redshift. We fit the stacked result presented by Henry et al. (2021) in the similar mass range, which assumes Curti et al. (2017) calibration. Our slope agrees with theirs β = 0.22 ± 0.03, but the intercept is ∼0.1 dex higher. This agrees with Wang et al. (2022b) and Li et al. (2023), who test that the calibrations of Bian et al. (2018) yielded a steeper MZR than the calibrations of Curti et al. (2017) when analyzing the same data.

Moreover, we compare our results with two simulation works presented separately in Figure 5. Our individual measurements are largely compatible with the result of the simulations IllustrisTNG (Torrey et al. 2019). But several high-metallicity galaxies lift the stacked MZR up high slightly, yielding a steeper slope than they predicted. Our measured slopes are in better agreement with the FIRE simulation results (Ma et al. 2016), which are capable of resolving high-z dwarf galaxies with sufficient spatial resolution.

In addition, all the MZRs discussed above are derived from galaxy populations residing in random fields. There has been continuous discussion about the environmental dependence of MZR shapes at high redshifts (Peng & Maiolino 2014; Bahé et al. 2017; Calabrò et al. 2022; Wang et al. 2023). Here we raise one recent observation of the MZR at z ∼ 2.2 showing a much shallower slope (β = 0.14 ± 0.02), measured using the HST grism spectroscopy of 36 galaxies residing in the core of the massive BOSS1244 protocluster (Wang et al. 2022b). Our work presented here confirms the significant difference between the MZR slopes measured in field and overdense environments, indicating the change in metal removal efficiency as a function of the environment.

Since metallicity estimates heavily rely on line flux measurements, in this section we verify that different methodologies in deriving emission line fluxes from the NIRISS slitless spectroscopy with limited spectral resolution do not result in significant biases on the metallicity derivations.

For grism spectroscopy, it has long been recognized that the morphological broadening effect can change the overall spectral shape and flux levels of galaxies (see, e.g., van Dokkum et al. 2011; Wang et al. 2019, 2020). We thus systematically compare, for the first time, two methods to measure emission line flux from slitless spectroscopy, with and without the consideration of this morphological broadening effect. The 2D forward-modeling analysis of Grizli is depicted in Section 3.1. In this section, we describe the line profile fitting to 1D extracted spectra using LMfit. The morphology of a galaxy has already been taken into account when forward modeling its 2D spectrum by Grizli. The extracted 1D spectra are morphologically broadened along the dispersion direction, and can vary significantly in spectral slope and flux level for the same object due to the different projected 1D morphology (see Figures 8 and 9 of Wang et al. 2019 for examples). Therefore, we regard the 2D line flux as the reference intrinsic value and 1D flux as the measurement not corrected for the morphology. The difference has not yet been fully investigated, and thus demands immediate attention, with the upcoming advent of large slitless spectroscopic surveys, e.g., Euclid, Roman, and the Chinese Space Station Telescope (CSST).

In the top three panels of Figure 6, we show the comparison between line flux measured from 2D or 1D spectra and try to associate it with the half-light radius r₅₀. The flux ratio of 2D to 1D deviates from 1 tangibly, and 2D flux modeled by Grizli are larger in most cases (47/48, 41/48, and 43/48 for [O iii], [O ii], and Hβ, respectively) than 1D flux fitted using LMfit by a median factor of ∼30% (with wide dispersion −0.3–5, where minus factor means 2D flux is lower than 1D flux). This strong offset does not seem to be related to S/N. As expected, we find it does correlate with the half-light radius r₅₀ of the individual galaxies, although not as strong as the Pearson correlation coefficients R shown. The unit of r₅₀ is the pixel, and here 1 pixel corresponds to 003, as illustrated in Section 2. Furthermore, Pearson R decreases as the S/N decreases from the first three brightest lines [O iii], [O ii] to Hβ, convincing us of this weak correlation. Linear fitting is employed in an attempt to describe this phenomenon, although it is based on limited data. This nonzero inconsistency first appears when we use 1D [O iii] flux to normalize our individual spectra for stacking. We rechecked our MZR using 2D [O iii] flux to normalize for stacking, and found the bias of metallicity is lower than 1σ. It indicates that the bias of the two flux measurements may be obscured by the stacking procedure, although we need a larger sample and more tests to verify this assertion. A more significant effect may be seen in the physical quantities directly determined by the line flux value, such as SFR.

Zoom In Zoom Out Reset image size

Since the flux ratio of 2D to 1D exhibits a correlation with the half-light radius r₅₀, we interpret this discrepancy as a morphological broadening effect. The morphological broadening of the spectrum is not due to physical factors such as velocity dispersion or radiative damping, but is simply an observational effect of the extended source (van Dokkum et al. 2011; Wang et al. 2020). For an ideal point source with no physical broadening effect, the emission line will be measured as a δ function. But if we could spatially resolve the galaxy, which is common in slitless spectroscopy, the emission line would be broadened as a result of the superposition of δ functions from individual pixels. Therefore, more parts of the line edge will be drowned in the noise, resulting in lower total line flux modeled by the Gaussian function. And of course, larger sources produce more broadening, yielding lower flux measurements. We, therefore, deem the top three panels of Figure 6 to be the first attempt to quantitatively analyze the impact of the morphological broadening effect. For large sources (r₅₀ > 10), the intrinsic flux could be several times larger than the broadened flux.

Although the 2D measurements are larger than the 1D results, in general, it seems that this bias is the same for all emission lines of the same source. As one can notice in the top three panels of Figure 6, for a given source with the same abscissa r₅₀, the corresponding ordinate values 2D/1D of all three lines are quite close to each other, although our naked eye can only recognize those outliers. And we have tested that these patterns are also independent of their S/Ns. Moreover, we show the line flux ratio in the bottom three panels of Figure 6, and they nearly follow the one-on-one line, with few outliers. That means even if this effect is not taken into account like in the 1D method, the flux ratios do not deviate from the 2D method significantly. Therefore, it indicates that the bias of the morphological broadening effect is systematic. We color-code them with the metallicity or the dust extinction A_v derived in Section 3.2 using 2D Grizli flux ratio. The color patterns demonstrate the physical meaning of these line ratios, i.e., the gas-phase metallicity diagnostics O₃₂ := [O iii]/[O ii], O₃ := [O iii]/Hβ, and the dust extinction indicator Hα/Hβ. The dotted line in the lower right marks the "intrinsic" line ratio in the absence of dust attenuation Hα/Hβ = 2.86. The few sources below it may be due to low S/N and measurement errors (see, e.g., Nelson et al. 2016).

As a consequence, our key result of the metallicity measurement derived from the ratio of two lines in Section 3.2 will not be greatly influenced by the 2D/1D flux measurement method. However the direct line flux (e.g., Hβ) and the derived quantity (e.g., SFR) of a single emission line could be biased, and for a large source, the intrinsic flux could be several times larger than the measured one. The coarse linear fitting here might describe the distinction between 2D/1D forward-modeling flux of emission line to some extent. We interpret this discrepancy as a morphological broadening effect. We recommend carefully checking the way flux is measured to match the scientific requirement and carefully forward modeling the spectrum through the convolution of the morphological broadening effect. The systematic offset, for the first time, may present an important guideline for future work deriving line fluxes with wide-field slitless spectroscopy, especially for large sky surveys to be conducted by, e.g., Euclid, Roman, and CSST, where it is time-consuming for 2D emission line modeling.