The Mass–Metallicity Relation of Dwarf Galaxies at Cosmic Noon from JWST Observations

NIRISS observations of the A2744 field were taken by the GLASS JWST Early Release Science (ERS) program (ERS-1324, Treu et al. 2022). The integration time is 2835 s (≃47 minutes) for pre-imaging and 2 × 5196.5 s ( ≃ 2.9 hr) for GR150 row and column (R+C) grism with three filters (F115W, F150W, F200W). This observation setup provides almost continuous wavelength coverage from 1 to 2.2 μm in the A2744 field.

SMACS 0723 field was targeted by the JWST early release observation program (ERO 2736, Pontoppidan et al. 2022). The JWST/NIRISS grism spectroscopy was obtained in F115W (1–1.3 μm) and F200W (1.7–2.2 μm) filters. The exposure time of the NIRISS grism observations of SMACS 0723 is 2 × 2834.5 s (≃1.6 hr) for each filter, respectively.

The JWST data were reduced and calibrated following the methodology of Wu et al. (2023). The data were reduced using the standard JWST pipeline ⁷ v1.9.4 with calibration reference files "jwst_1041.pmap." First, we modeled the 1/f noise (see Schlawin et al. 2020) and removed it using the tshirt/roeba code. ⁸ Then, we identified "snowball" artifacts from cosmic rays (Rigby et al. 2023) and masked them. We then registered the world coordinate system of mosaicked images using the Pan-STARRS Data Release (DR1) catalog (Flewelling et al. 2020; Gaia DR2 catalog, Gaia Collaboration et al. 2018) for A2744 (SMACS 0723) field. Lastly, the pixel scale of final mosaicked images was resampled to 003 with pixfrac = 0.8. We reduced NIRISS grism data using the Grism Redshift & Line Analysis tool (Grizli; ⁹ Brammer et al. 2022). The 2D grism spectra were drizzled with a pixel scale of 0065.

The JWST/NIRISS grism observations cover the metallicity sensitive lines [O ii] λ λ3727, 3729, Hβ, and [O iii] λ λ4959, 5007, allowing for our gas-phase oxygen abundance estimate. The spectral resolution (R ≡ λ/Δλ) of JWST/NIRISS observations is ≈150, enabling the separation of the [O iii] λ λ4959, 5007 doublet.

We used Grizli to perform target detection, contamination modeling, and spectrum extraction. First, we constructed pre-imaging mosaics, which were used to generate catalogs of sources and associated segmentation maps. Then, we built contamination models for all sources by forward modeling, which eliminates both zeroth-order and higher-order spectra from nearby sources. Finally, we extracted the contamination-removed spectra for detected sources.

We determined galaxy redshifts by spectral template synthesis (SSP) fitting. From the extracted sources, we generated a catalog of galaxies with robust redshift solutions, which required only one significant solution ($\mathrm{log}(p(z)/{\chi }^{2})\gt 1$) for the redshift template fitting. Specifically, at z ≃ 2–3, we expect multiple strong emission lines within the NIRISS wavelength coverage, e.g., [O ii], Hβ, [O iii], and Hα. The redshifts can be well fitted by the resolved [O iii] doublets, and the secondary lines (e.g., [O ii], Hβ, and Hα) allow us to perform a more robust redshift fitting.

Because A2744 cluster is one of the Hubble Space Telescope Frontier Fields (HFF) programs ¹⁰ (e.g., Lotz et al. 2017) and the SMAC 0723 field is part of the Reionization Lensing Cluster Survey (RELICS; Coe & Salmon 2019), ¹¹ each cluster has extensive archival photometric imaging across a range of wavelengths. Specifically, these include HST/Advanced Camera for Surveys (ACS) in F435W, F606W, F814W filters, and the Wide Field Camera 3 (WFC3)/IR in F105W, F125W, F140W, F160W filters. In this study, we use the high-level science products released by Space Telescope Science Institute with a pixel size of 30 mas.

Both A2744 and SMACS 0723 have grism pre-imaging from JWST/NIRISS as references to extract the grism spectra (see Section 2.1). Additionally, SMACS 0723 field has deep JWST/Near Infrared Camera (NIRCam) imaging data in the F090W, F150W, F200W, F277W, F356W, and F444W bands.

For HST imaging in the A2744 field, we used photometric results from the HFF-DeepSpace multiple wavelength catalogs (e.g., Shipley et al. 2018; Nedkova et al. 2021). For HST imaging in the SMACS 0723 field and JWST imaging data in both fields, we use SExtractor (Bertin & Arnouts 1996) to measure Kron-like AUTO photometry of the galaxies on our custom calibrated JWST images and public HST data product. We performed individual measurements and subsequently crossmatched them to obtain a photometric catalog using the detection of NIRISS F200W as the reference. Finally, we built the combined catalog with photometry in F435W, F606W, F814W, F105W, F115W, F125W, F140W, F150W, F160W, and F200W bands for both fields. For the SMACS 0723 field, we also include the F090W, F277W, F356W, and F444W photometry. A 5% photometry relative uncertainty floor is set for JWST bands to involve the systematic zero-point error. Based on this photometric catalog, we perform spectral energy distribution (SED) fitting to derive the physical properties of each galaxy, as explained in Section 3.2.

We selected star-forming galaxies based on the JWST/NIRISS wide field slitless spectroscopy. We used the following criteria to select galaxies from the source catalog as established in Section 2.1.

• 1.  
  With the redshifts determined by NIRISS grism, we select galaxies at 1.8 < z < 2.3 and 2.65 < z < 3.4. For the SMACS 0723 field, the absence of grism observations in the F150W band limited us to only select galaxies at 2.65 < z < 3.4. These redshift ranges ensure that NIRISS spectra completely cover the [O ii], Hβ, and [O iii] emission lines in the A2744 field, and Hβ and [O iii] emission lines in the SMACS 0723 field. Our redshift selection criteria also ensure these emission lines are located in the high response wavelength interval (>90% response at filter center), which avoids the impact of large flux calibration uncertainty. Based on this redshift criterion, we selected 76 galaxies to be the primary sample to apply the following criteria.
• 2.  
  The S/N of the [O iii] line should be greater than 10. This criterion was established by empirical inspection. We did not require the S/N of secondary lines, which could lead to the Eddington bias (Eddington 1913). The galaxies that have secondary line detection (e.g., Hβ or [O ii]) tended to have higher metallicity (e.g., Henry et al. 2021). To minimize this selection effect, only [O iii] detection significance was required instead. 72 galaxies satisfied this criterion.
• 3.  
  To remove active galactic nucleus (AGN) contamination, we excluded galaxies that fall in the "AGN regime" of the mass-excitation (MEx) diagram ([O iii] λ5007/Hβ versus M_⋆) defined by Juneau et al. (2011) and modified by Coil et al. (2015). In all galaxies, we found Hβ to be at least marginally detected (>1.5σ), so there are very few sources with ambiguous [O iii]/Hβ on the MEx diagram. Eight galaxies were classified in the AGN-like category and were removed from our sample. Because we detected no other AGN diagnostic lines (e.g., [S ii] and [N ii]), the MEx diagnostic was the best method available to remove AGN candidates from our sample.
• 4.  
  We performed visual inspections on each individual galaxy and removed heavily contaminated sources. Due to the challenges of extracting the spectrum of multiimages and arclets, we removed these sources following Mahler et al. (2018), Bergamini et al. (2023).

Ultimately, we selected 51 galaxies, where 44 are in A2744, and seven are in SMACS 0723. We show six examples of the JWST/NIRISS color composite stamps and spectra in Figure 1(a). Among all galaxies in the sample, 27 galaxies are at redshift 1.8 < z < 2.3 (z_med = 2.01), and 24 galaxies are at redshift 2.65 < z < 3.4 (z_med = 2.97). The redshifts of 42 galaxies are spectroscopically determined by NIRISS observations for the first time, while the remaining nine galaxies have been previously confirmed by HST/WFC3 or Very Large Telescope/MUSE (e.g., Mahler et al. 2018; Wang et al. 2020; Boyett et al. 2022). Among the nine galaxies, the NIRISS redshifts of two galaxies (ID 1692 and ID 2225) are inconsistent with those in previous studies. For these two galaxies, we rechecked their spectra (shown in Figure 1(a)) and confirm the reliability of the NIRISS redshifts based on the appearance of multiple emission lines with S/N > 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To measure the intensity of galaxy emission lines, we use the Gaussian profile models to fit the lines in continuum-subtracted grism spectra. The line wavelengths are fixed by the galaxy's redshift, which is determined by the SSP fitting as presented in Section 2.1. For the Hα and [O ii] lines, we fit the spectra using a single 1D Gaussian function. For the Hβ + [O iii]λ λ4959, 5007 lines, we fit the spectra using three Gaussian functions. We fix the line ratio of the [O iii] doublet to [O iii]λ4959/[O iii] λ5007 = 1/3 and require the line width of [O iii] doublet to be equal. We further verify the accuracy of our derived flux by comparing our measurements to the modeled flux solution fitted by Grizli and find the two to be consistent. A thorough summary of the best-fit observed emission line fluxes can be found in Appendix.

In the early era of JWST observations, data calibration cannot be completely accurate (e.g., Rigby et al. 2023). To minimize calibration issues, we employed reference calibration files that have been corrected using space-based observations. To quantify the effect of any remaining systematic errors on our data, we compared the NIRISS-extracted spectra with previous HST/WFC3 grism observations. We found that the estimated relative flux calibration error between different filter bands was limited to ≲5%, which is consistent with values reported in previous literature (e.g., Wang et al. 2022a).

The galaxy stellar masses are derived by SED fitting to the photometry shown in Section 2.2. Note that both A2744 and SMACS 0723 clusters significantly magnify the background sources, including galaxies in this work due to the gravitational lensing. In order to calculate the lensing magnification of each galaxy at the appropriate redshift, we adopt the lensing models from Mahler et al. (2018), Golubchik et al. (2022), respectively (see Appendix). To do our analysis, we first correct the magnitudes in all bands by the lensing magnification factor. Note that we do not introduce the systematic error of the lensing models. To simultaneously measure the stellar mass (M_⋆), the stellar dust attenuation (E(B − V)_star), and the best-fit model of the stellar continuum with emission lines, we fit the SED using the Bayesian code BEAGLE (Chevallard & Charlot 2016). For our SED fitting, we fix the galaxy redshift, assume a constant star formation history, and model the dust attenuation curve following Calzetti et al. (2000). We set the dust attenuation in the V band A_V to vary from 0–8 and use the Chabrier (2003) IMF with an upper limit of 300 M_⊙. The metallicity of the interstellar medium (ISM) is assumed to be the same as that of the stellar populations. The nebular line and continuum emission are included in our fitting based on the CLOUDY photoionization code (Ferland et al. 2017). To model the metallicity, we use a flat prior with values ranging from $-2\lt \mathrm{log}(Z/{Z}_{\odot })\lt 0.4$. For our model, we bound the ionization parameter in log-space between $-4\lt \mathrm{log}U\lt -1$.

The derivation of star formation rate (SFR) and gaseous metallicity is based on the attenuation-corrected emission line intensities. We correct the reddening effect of dust attenuation for all measured nebular lines using the E(B − V)_star, which is derived from the SED fitting, assuming the dust reddening for the nebular and stellar components are equal. Indeed, several works (e.g., Shivaei et al. 2020) pointed out that E(B − V)_star could equal to ∼2E(B − V)_gas for galaxies at z = 2; we test using this assumption and found out that the SFR could be increased by 0.1 dex. Furthermore, the relative changing line ratios (e.g., O₃₂) are at an average level of 0.03 dex, which will lead to consistent conclusions in this paper. Therefore, in the rest of the paper, we adopt the assumption of E(B − V)_star = E(B − V)_gas.

Note that the nebular dust attenuation can be also derived from Hα/Hβ Balmer decrement (e.g., Hα/Hβ or Hβ/Hγ) for a few galaxies in our sample. But half of our galaxy sample does not have Hα coverage. The absence of Hα line coverage and lack of significant Hγ line detection make it impossible to measure dust attenuation via the Balmer decrement. To maintain consistency among our entire sample, we adopt the dust attenuation from the SED fitting to perform the dereddening correction for emission lines.

We measure the SFR using the attenuation-corrected Hα luminosity. For galaxies without Hα coverage, we derive the Hα luminosity from the attenuation-corrected Hβ flux by adopting a Balmer decrement of Hα/Hβ = 2.86. In order to estimate the SFR, we assume the Kennicutt (1998) SFR calibration and the Chabrier (2003) IMF. The SFR is modeled as

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

where L_Hα is the Hα luminosity. Throughout this paper, SFR corresponds to the SFR_Hα , which is derived from the attenuation-corrected Hα luminosity. And for those galaxies without Hα coverage, we get Hα luminosity using the attenuation-corrected Hβ assuming a Balmer decrement of Hα/Hβ = 2.86.

Within our sample, it is worth mentioning that the Hα line is blended with [N ii] lines due to the low spectral resolution, which biases our estimate of Hα toward a higher value. However, the contribution of [N ii] decreases with metallicity, and typically, the [N ii]/Hα is less than 0.07 for star-forming galaxies with the stellar mass of M_⋆ < 10⁹ M_⊙ at z = 2–3 (e.g., Bian et al. 2018). Therefore, the [N ii] contamination should be negligible to our SFR measurement, especially when compared with the uncertainty of dust attenuation correction. Besides, Balmer line stellar absorption can bias the intrinsic line luminosity to a lower value. Nevertheless, for the young star-forming galaxies in our sample, the Balmer line stellar absorption is comparatively weak and marginalized.

Using the measurements from SED fitting and SFR calibration, we plot Log(SFR) as a function of the Log(M_⋆), in Figure 2. We find that the relation between SFR and M_⋆ displays a positive correlation across the complete range of stellar masses and redshift. We also find this slope agrees with the more massive galaxies measured in the MOSDEF sample from Sanders et al. (2021) and the extrapolated galaxy star-forming main sequence. Despite a large degree of scatter, the SFR at fixed M_⋆ increases with the redshift (Figure 2). However, when compared with previous studies of star-forming main-sequence galaxies at redshifts of 2 < z < 3, our sample of galaxies has a slightly higher mean SFR.

Zoom In Zoom Out Reset image size

Bian et al. (2018) established empirical relations between the strong metallicity diagnostic line ratios and the direct-T_e metallicity using a sample of local analogs of high-redshift galaxies. These calibrations have been demonstrated as one of the most reliable metallicity calibrations at z > 1 (Sanders et al. 2020).

We derive the gas-phase oxygen abundance by adopting the Bian et al. (2018) calibration and measuring the strong metallicity diagnostic line ratios O₃₂ = log([O iii]λ λ4959, 5007/[O ii]λ λ3727, 3729), and O₃ = log([O iii]λ λ4959, 5007/Hβ). This calibration allows us to minimize the systematic uncertainty of the metallicity measurements when comparing our MZR results with those from Sanders et al. (2021). For the galaxies with both [O iii] and [O ii] coverage (44/51), we use the O₃₂ to measure the gas-phase oxygen abundance as follows:

$$\begin{eqnarray}&&12+\mathrm{log}({\rm{O}}/{\rm{H}})=8.54-0.59\times {{\rm{O}}}_{32}.\end{eqnarray} \tag{ 2 }$$

As shown above, we use the Bian et al. (2018) calibration, which is suitable in the range of 0.3 < O₃₂ < 1.2 to measure the oxygen abundance in our galaxies. For the galaxies (11/48) where O₃₂ is out of this range, we extrapolate this relation linearly to estimate their oxygen abundance. The linearity of the O₃₂ calibration allows us to extrapolate with minimal systematic error. For the same reason, we try not to use the R₂₃ and O₃ to further constrain the metallicity, which will cause more systematic uncertainties, especially for low-mass and low-metallicity galaxies.

For galaxies without [O ii] coverage (7/51), we calculate the metallicity by using line ratio O₃ ≡ [O iii]λ λ 4959,5007/Hβ, and use the following equation:

$$\begin{eqnarray}&&{{\rm{O}}}_{3}=43.9836-21.6211x+3.4277{x}^{2}-0.1747{x}^{3},\end{eqnarray} \tag{ 3 }$$

where $x=12+\mathrm{log}({\rm{O}}/{\rm{H}})$.

For the above O₃–Z relation, there exist two oxygen abundance solutions for a given O₃ value. We adopt the solution within the range of $7.8\lt 12+\mathrm{log}({\rm{O}}/{\rm{H}})\lt 8.4$, which is consistent with the oxygen abundance range derived from O₃₂ indicator in a similar mass range. Our derived metallicities are presented in Appendix.

Additionally, we stack each individual spectrum to increase the S/N and measure the average metallicity in different stellar mass and redshift bins. We perform the spectrum stacking in two redshift bins (1.8 < z < 2.3 and 2.65 < z < 3.4). For our stacking analysis, we only use the galaxies in the A2744 field, because those in the SMACS 0723 field do not have [O ii] coverage. We divide the entire galaxy sample into four (three) stellar mass bins within the redshift range 1.8 < z < 2.3 (2.65 < z < 3.4). Each stellar mass bin contains between four and nine galaxies, which is determined by the bootstrap simulation to minimize the uncertainty of derived metallicities of stacked spectra. Moreover, we ensure uniformity of the median stellar mass of each mass bin when setting the galaxy number in each bin.

For our stacking analysis, we first shift each spectrum into the rest frame according to the grism redshift. We then apply dust-corrections following the E(B − V) values from our SED fitting and the Calzetti et al. (2000) attenuation curve. We then subtract off the continuum for each individual spectrum before normalization and stacking. We normalized each spectrum based on their [O iii] flux and resampled it onto a uniform wavelength grid. Lastly, the spectra in each mass bin are combined by taking the mean of the 3σ clipping of each wavelength grid. In Figure 1(b), we present the final stacked spectra in different stellar mass and redshift bins. We measure the relative emission line flux in the stacked spectra using the same method as for the individual galaxies. We derive gas-phase metallicities using the O₃₂ indicator. The measurement uncertainties are propagated from both the stacking and line flux fitting. The median stellar mass, measured line ratios, median SFR, and the derived metallicity of stacked spectra are presented in Table 1.

Table 1. Properties of Stacked Spectra in Bins of M_⋆ for the z ≃ 2.0 and z ≃ 3.0 Samples

log (M⋆/M⊙)	N	($12+\mathrm{log}({\rm{O}}/{\rm{H}})$)stack	($12+\mathrm{log}({\rm{O}}/{\rm{H}})$)med	O32	O3	SFR [M⊙ yr−1]
(1)	(2)	(3)	(4)	(5)	(6)	(7)
1.8 < zgrism < 2.3
6.64 ± 0.12	8	8.03 ± 0.05	8.00 ± 0.10	0.86 ± 0.08	0.59 ± 0.04	0.61 ± 0.12
7.56 ± 0.09	9	8.11 ± 0.03	8.05 ± 0.06	0.73 ± 0.06	0.79 ± 0.04	1.09 ± 0.25
8.50 ± 0.09	6	8.23 ± 0.03	8.26 ± 0.06	0.52 ± 0.05	0.74 ± 0.04	2.61 ± 0.95
9.25 ± 0.20	4	8.35 ± 0.02	8.38 ± 0.07	0.32 ± 0.03	0.55 ± 0.03	2.83 ± 0.57
2.65 < zgrism < 3.4
7.48 ± 0.09	7	8.00 ± 0.04	7.98 ± 0.08	0.91 ± 0.06	0.75 ± 0.03	5.02 ± 1.77
8.39 ± 0.10	6	8.14 ± 0.03	8.10 ± 0.05	0.67 ± 0.04	0.78 ± 0.03	5.60 ± 2.35
9.62 ± 0.06	4	8.34 ± 0.02	8.40 ± 0.09	0.34 ± 0.03	1.01 ± 0.10	10.64 ± 6.97

Note. Column (1) is the median stellar mass for each bin; column (2) is the number of galaxies in each bin; columns (3) and (4) are gas-phase metallicity measured from the stacked spectra and statistical median of individual galaxies, respectively; columns (5) and (6) are the line ratios O₃₂ = log([O iii]λ λ4959, 5007/[O ii]λ λ3727, 3729), and O₃ = log([O iii]λ λ4959, 5007/Hβ) derived from the stacked spectra; column (7) is the median star formation rate.

Download table as:  ASCIITypeset image

It is important to note that performing the dust-correction either before or after stacking during the stacking procedure could potentially introduce bias to the measured metallicity (e.g., Henry et al. 2021). We performed a test by changing the sequence of dust-corrections. Consequently, we find that the difference in derived metallicity 12+log(O/H) is 0.01–0.03 typically, which is within 1σ error of metallicity ≃0.02–0.04. Therefore, our results are not sensitive to the dust-correction sequence when stacking. We also measure the dust attenuation using the Balmer decrement on the stacked spectra that are not dust-corrected.

We find a clear correlation between gas-phase metallicity and stellar mass for our measurements of both individual galaxies and composite spectra. This is the first-ever effort to derive the MZR down to galaxies of M_⋆ = 10^6.5 M_⊙, at z = 2–3 using grism slitless spectroscopy. Within both the z ≃ 2 and z ≃ 3 sample, we find Spearman correlation coefficients of 0.70 and 0.68, respectively, with the statistical-significant p-value <10⁻⁴. This means a positive correlation with strong evidence to reject the null hypothesis.

Beyond the strong trends of the MZR within each sample, we also look for signs of redshift evolution. Using the stacked spectra for the two samples, we find that the MZRs at z ≃ 2.0 and z ≃ 3.0 show a monotonic evolution toward lower metallicity with increasing redshift. As seen in Figure 3, the scale of this evolution is similar to that of MOSDEF from z ≃ 2.3 to z ≃ 3.3.

Zoom In Zoom Out Reset image size

To quantify the MZR and derive slope and normalization, we fit it using a power law of the following form:

$$\begin{eqnarray}&&12+\mathrm{log}({\rm{O}}/{\rm{H}})=\beta \times \mathrm{log}\left(\displaystyle \frac{{M}_{\star }}{{10}^{8}\,{M}_{\odot }}\right)+Z8,\end{eqnarray} \tag{ 4 }$$

where Z8 is the metallicity at the stellar mass of 10⁸ M_⊙. Using all galaxies in the sample, we get the best-fit relations with 1σ uncertainties:

$$\begin{eqnarray}&&\begin{array}{l}12+\mathrm{log}({\rm{O}}/{\rm{H}})\\ =\,(0.14\pm 0.04)\times \mathrm{log}\left(\displaystyle \frac{{M}_{\star }}{{10}^{8}\,{M}_{\odot }}\right)+(8.14\pm 0.03).\end{array}\end{eqnarray} \tag{ 5 }$$

We also fit the stacked spectra in the z ≃ 2.0 and z ≃ 3.0. The best-fit parameters with 1σ uncertainties are as follows: for z ≃ 2.0, β_z≃2 = 0.13 ± 0.04, and Z8_z≃2 = 8.18 ± 0.02; and for z ≃ 3.0, β_z≃3 = 0.16 ± 0.03, and Z8_z≃3 = 8.08 ± 0.03. These two best-fit relations and the M–Z values of each individual galaxy are shown in Figure 3.

Note that the slope of best-fit MZR for all stacked spectroscopy is 0.14 ± 0.04, consistent with both z ≃ 2 and z ≃ 3 best-fit slopes, implying that there is no significant evolution in the slope between z = 2–3. In contrast to the lack of evolution in the slope, we do find moderate evolution in the MZR normalization with redshift between z = 2–3. Specifically, we find that the difference in the normalization at z ≃ 2 and z ≃ 3 is ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})/{\rm{\Delta }}z=-0.09\pm 0.06$. Beyond redshift evolution, we also probe for differences in the MZR as a function of stellar mass. We find that, at the low-mass end, we measure a much shallower slope than at the high-mass end mass where β ≃ 0.30 (Sanders et al. 2021). This suggests a transition in the slope of the MZR at M_⋆ ≃ 10⁹ M_⊙.

Although our sample is not large enough to allow a statistically robust study of the scatter of the MZR at the low-mass end, we find a tentative increase of the MZR scatter at the stellar mass less than 10⁸ M_⊙. By defining the intrinsic scatter as ${\sigma }_{\mathrm{ins}}=\sqrt{{\sigma }_{\mathrm{obs}}^{2}-{\sigma }_{\mathrm{meas}}^{2}}$, where σ_obs is the observation scatter, and σ_meas is the measurement uncertainties, we measure the intrinsic scatter to be 0.16–0.18 dex at M_⋆ = 10⁸–10⁹ M_⊙. However, we see this value increase to ≈0.23 dex at M_⋆ = 10⁷ M_⊙ bins implying a tentative evolution in the scatter.

Deep JWST/NIRISS spectroscopy of a sample of dwarf galaxies allows us, for the first time, to make a robust measurement of the slope of MZR at z = 2–3 using low-mass galaxies (M_⋆ ≈ 10^6.5–9.5 M_⊙). At z = 2–3, we find a shallower slope in MZR at the low-mass end, in contrast to that from Sanders et al. (2021), which examines a sample of galaxies with M_⋆ = 10⁹–10¹¹ M_⊙. This difference indicates a transition in MZR slope around M_⋆ ≃ 10⁹ M_⊙.

The differences in our measurements of the MZR slope may be indicative of differences in the feedback mechanisms that regulate it (e.g., Wang et al. 2022b; Strom et al. 2022). Typically, stellar feedback drives the gas and metal out of galaxies, which will result in differences in the gas-phase metallicity of the ISM. The feedback of stellar winds can be described by the wind mass loading factor ($\eta \equiv {\dot{M}}_{\mathrm{out}}/\mathrm{SFR}$), which presents the fraction of gas ejected in winds. The wind mass loading factor η is related to the wind velocity (v_w ), with $\eta \propto {v}_{w}^{-2}$ for the energy-driven wind model, $\eta \propto {v}_{w}^{-1}$ for the momentum-driven wind model, and $\eta \equiv \mathrm{const}$ for the constant wind model. We argue that the shallower MZR slope at the low-mass end is caused by galaxies with different masses that could be dominated by different feedback mechanisms.

For example, similar to our analysis, Torrey et al. (2019) used IllustrisTNG simulations and similarly found that the MZR has a shallower slope at the low-mass end. They argue that their results are indicative that feedback in low-mass galaxies may be dominated by constant wind ($\eta =\mathrm{constant}$), while feedback in high-mass galaxies may be dominated by energy-driven wind ($\eta \propto {v}_{w}^{2}$). Similarly, Pillepich et al. (2018) presents the necessity to introduce a constant wind model to match the low-mass end of the galaxy stellar mass function. Specifically, the constant wind velocity reduces the wind mass loading factor. The larger-than-expected wind velocity and lower mass loading factor further result in longer wind recycling times and less direct metal ejection, which manifests as somewhat higher metallicities in the ISM and a shallower MZR slope for low-mass galaxies, in agreement with our results.

Our best fit MZR slope is β = 0.14 ± 0.04, which is consistent with the prediction of the momentum-driven wind model (e.g., Finlator & Davé 2008; Guo et al. 2016). Such a slope suggests that momentum-driven wind is the dominant feedback in dwarf galaxies, in contrast to the constant wind models in Torrey et al. (2019). One possible interpretation is that the measured MZR slope in our sample is the average effect of a transition between the energy-driven wind model and constant wind models (the slope of the MZR in momentum-driven wind models is between the other two models). In this interpretation, the momentum-driven wind model ($\eta \propto {v}_{w}^{-1}$) acts as an intermediate state between the energy-driven wind ($\eta \propto {v}_{w}^{-2}$) and constant wind ($\eta \equiv \mathrm{const}$). Thus, a more general power-law wind model could be written as $\eta \propto {v}_{w}^{-{\alpha }_{\mathrm{fb}}}$ where α_fb is the power index parameter of the feedback. Our results indicate that there could exist a continuous transition of wind model power index from α_fb = 2, to α_fb = 0 in low-mass galaxies with a transition stellar mass of M_⋆ ≃ 10⁹ M_⊙.

The slope of the MZR at the low-mass end remains not well constrained and is under debate. For example, Henry et al. (2021) found that the MZR slope for galaxies with stellar masses M_⋆ > 10⁸ M_⊙ at z ≃ 1.9 could remain steep, similar to that of higher-mass galaxies. However, the low-mass MZR slope we find is shallower. The discrepancy with Henry et al. (2021) could arise from several factors, including the different metallicity calibration methods used. As Sanders et al. (2020) suggested, we employ strong line calibration from Bian et al. (2018) and use O₃₂ as the metallicity indicator. In contrast, Henry et al. (2021) derived metallicities based on Curti et al. (2017) using the joint fitting of multiple line ratios. The choice of strong lines and a calibration method could lead to a systematic discrepancy in MZR slope and normalization. To better compare our results with those from Henry et al. (2021), we recalculated the metallicities from their stacked measurements (overlaid in Figure 3), using our method based on the Bian et al. (2018) calibration. We find that the stacked results at M_⋆ < 10⁹ M_⊙ show a tendency to follow the shallower-sloped MZR compared to the higher-mass points. To better constrain the validity of the slope transition hypothesis, we require further observations of the low-mass MZR and a more robust comparison to simulations.

For fixed stellar mass M_⋆, we find that the MZR normalization moderately evolves from lower O/H to higher O/H from z = 3, to z = 2. Moreover, we also aim to determine any evolution with stellar mass by comparing with Sanders et al. (2021). For the high-mass end (M_⋆ ≃ 10¹⁰ M_⊙), the normalization evolution is ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})/{\rm{\Delta }}z=-0.11\pm 0.02$, while for the result, in this work, we found ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})/{\rm{\Delta }}z\,=-0.09\pm 0.06$ (M_⋆ ≃ 10⁸ M_⊙). It concludes that the evolution of the MZR is independent of galaxy stellar mass, as indicated by an unchanged MZR shape from z = 3, to z = 2 in a wide stellar mass range of 10⁶–10¹⁰ M_⊙.

The evolution of MZR is the consequence of a complicated interplay between star formation efficiency, galactic feedback, gas fraction, gas accretion, and recycling (e.g., Lilly et al. 2013). Thus, studies of MZR evolution can be used to constrain how the above processes evolve with cosmic time. Along these lines, Onodera et al. (2016) explain the evolution of the MZR normalization by suggesting that star formation efficiency evolves weakly with redshift. However, Sanders et al. (2021) suggest an increasing star formation efficiency with redshift to match their mild MZR evolution (${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})/{\rm{\Delta }}z\simeq -0.1$). An increase in star formation efficiency with redshift is consistent with recent studies on the evolution of star formation efficiency (e.g., Liu et al. 2019). In this work, we also find a similar mild MZR evolution (${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})/{\rm{\Delta }}z=-0.09\,\pm 0.06$) at the low-mass end, which suggests that star formation efficiency increases with redshift across a wide stellar mass range.

The gas fraction of galaxies is another crucial factor for regulating the metallicity of galaxies. Observations have found that the gas fraction increases significantly with redshift at fixed M_⋆ and has been quantified by μ_gas ≡ M_gas/M_⋆ ∝ (1 + z)^2.5 (e.g., Tacconi et al. 2018). These large gas fractions can dilute the ISM of galaxies making a lower gas-phase metallicity for higher-redshift galaxies. Theoretically, the evolution of the gas fraction is considered the driving force of MZR evolution in some cosmological simulations (e.g., Ma et al. 2016; Torrey et al. 2019). However, the gas fraction is also predicted to increase with decreasing stellar mass at fixed redshift, which is quantified as ${\mu }_{\mathrm{gas}}\propto {M}_{\star }^{-0.4}$ for z = 2 (e.g., Torrey et al. 2019). In contrast to previous studies, our results indicate that the evolution rate of MZR appears to be the same in the low-mass and high-mass end, implying that the high gas fraction of dwarf galaxies at z = 2–3 may have a weak effect on the evolution of MZR. In the future, the direct measurements of the gas fraction in dwarf galaxies (e.g., using the Square Kilometre Array) are crucial to understanding the effect and contribution of gas reservoirs on galaxy metal enrichment.

Previous works (e.g., Ellison et al. 2008; Mannucci et al. 2010; Andrews & Martini 2013) have suggested that a fundamental relation exists between metallicity and other galaxy properties. This FMR correlates the metallicity with both stellar mass and SFR. For a given stellar mass, the metallicity is enhanced with the decreasing SFR. Several works probing large samples of galaxies from z = 3, to z = 0 have found little redshift evolution in the FMR in galaxies above 10⁸ M_⊙ (e.g., Henry et al. 2013; Sanders et al. 2015, 2018; Curti et al. 2020; Henry et al. 2021). These results emphasize the dominance of stellar mass and SFR on the gas-phase metallicity of galaxies, regardless of the cosmic redshift.

The FMR relation can be explicitly parameterized using a 2D projection with the relation between O/H and is written as

$$\begin{eqnarray}&&{\mu }_{\alpha }\equiv \mathrm{log}\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)-\alpha \times \mathrm{log}\left(\displaystyle \frac{\mathrm{SFR}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right).\end{eqnarray} \tag{ 6 }$$

In the above equation, α is introduced as the parameter to minimize the scatter in O/H. For the dwarf galaxies at z = 2–3 in our sample, the scatter of O/H is minimized at α_best = 0.60 ± 0.13, suggesting the existence of FMR for z = 2–3 dwarf galaxies.

Based on our fitting, we construct an FMR with the linear function form of

$$\begin{eqnarray}&&12+\mathrm{log}({\rm{O}}/{\rm{H}})=(0.20\pm 0.02){\mu }_{{\alpha }_{\mathrm{best}}}+(6.59\pm 0.15),\end{eqnarray} \tag{ 7 }$$

where ${\mu }_{{\alpha }_{\mathrm{best}}}=\mathrm{log}\left({M}_{\star }/{M}_{\odot }\right)-0.60\mathrm{log}\left(\mathrm{SFR}/{M}_{\odot }\,{\mathrm{yr}}^{-1}\right)$. We show the projection of the FMR, the O/H versus ${\mu }_{{\alpha }_{\mathrm{best}}}$ in Figure 4. As seen in Figure 4, we find that the FMR exists within a large stellar mass range (M_⋆ = 10⁶–10¹⁰ M_⊙) at z = 2–3.

Zoom In Zoom Out Reset image size

The existence of the FMR in the low stellar mass end demonstrates that the SFR and stellar mass act as the main physical properties for determining the gas-phase metallicity in dwarf galaxies at z = 2–3. However, it is important to note that the exact measurement of parameter α and interpretation of the FMR can vary depending on the method used to measure it. For example, some works that have measured alpha for the FMR have found quite a range, with a systematic difference depending on whether the FMR was measured using a strong line calibration or direct metallicities in individual galaxies or stacks (e.g., Andrews & Martini 2013; Henry et al. 2021). In the context of these studies, the results presented in this paper suggest that the FMR for dwarf galaxies at z = 2–3 follows a consistent α value with previous measurements (α =0.55–0.7; Andrews & Martini 2013; Sanders et al. 2017; Curti et al. 2020; Sanders et al. 2021). These measurements are over cosmic time (z = 0–3), suggesting a slight evolution of FMR. However, we note that the slope of FMR in the low-mass end is much shallower than previous measurements, and further studies will be needed to fully understand the evolution of the FMR over cosmic time.

In the future, deeper imaging and spectroscopic observations over a larger survey area will allow us to probe the MZR and FMR relations out to even higher redshifts and down to an even lower metallicity. This will allow us to better understand the physical processes driving chemical enrichment at different epochs in the history of the Universe.