A Large Sample of Extremely Metal-poor Galaxies at $z<1$ Identified from the DESI Early Data

The DESI experiment is a new generation ground-based cosmological survey, aiming to explore the expansion history and structure growth rate of the universe with unprecedented precision and to understand the nature of dark energy (Levi et al., 2013; DESI Collaboration et al., 2016a, b, 2023a, 2023b). It will accurately measure the redshifts of about 40 million galaxies and quasars over a 14,000 square degree footprint. The instrument is a multi-object spectroscopic system installed at the prime focus of the Mayall 4m telescope at Kitt Peak, Arizona (DESI Collaboration et al., 2022; Silber et al., 2023; Miller et al., 2023). It features a 3.2-deg field-of-view and 5000 fiber robots that can position optical fibers to 5000 objects and collect their spectra simultaneously. The fibers (diameter of 1.5″) are fed to ten 3-arm spectrographs, which provide a continuous wavelength coverage of 3600–9800Å. The spectral resolution varies from $R\sim$ 2000 in the blue end to 5000 in the red, which is high enough to resolve the [O II] doublet.

The DESI spectroscopic targets include approximately 50 million extragalactic and Galactic objects, mainly selected from the DESI Legacy Imaging Surveys (Zou et al., 2017; Dey et al., 2019; Schlegel et al., 2023). The imaging data are about 2 mag deeper than the SDSS, covering a sky area of 20,000 deg² in optical $grz$ bands and including the latest six-year near-infrared observations from the Wide-field Infrared Survey Explorer (WISE) satellite. The primary targets to be observed in dark time are LRGs at $0.4<z<1.1$, ELGs at $0.6<z<1.6$ and QSOs at $z<3.5$ (Raichoor et al., 2020; Ruiz-Macias et al., 2020; Yèche et al., 2020; Zhou et al., 2020; Chaussidon et al., 2023; Zhou et al., 2023; Raichoor et al., 2023b). During bright time, DESI will carry the Bright Galaxy Survey (BGS) of galaxies at $z<0.6$ and Milky Way survey of stars (Allende Prieto et al., 2020; Ruiz-Macias et al., 2020; Cooper et al., 2022; Hahn et al., 2023). There are also some secondary target programs to utilize the spare fibers and DESI capabilities, such as the LOW-Z Secondary Target Survey targeting low-redshift ($z<0.03$) dwarf galaxies (Darragh-Ford et al., 2022).

The DESI project is expected to be a five year survey. Before the main survey started on May 14, 2021, there were several stages of survey validation for about six months (DESI Collaboration et al., 2023a), which aimed to refine the target selections and fiber assignment, improve the data reduction pipeline and observing efficiency, and validate the instrumental performances (Alexander et al., 2023; Bailey et al., 2023; DESI Collaboration et al., 2023a; Kirkby et al., 2023; Raichoor et al., 2023a; Schlafly et al., 2023; Lan et al., 2023), etc. An automatic spectroscopic data pipeline has been developed to reduce raw observations into wavelength- and flux-calibrated spectra, and to perform redshift measurements for all observed targets (Guy et al., 2023). The spectroscopic classifications and redshifts are determined using the Redrock pipeline¹¹1https://github.com/desihub/redrock (Bailey et al., 2023; Brodzeller et al., 2023), which fits a suite of templates of stars, galaxies, and quasars to the DESI spectra. There are several internal data releases generated by these data reduction and analysis pipelines, such as Fuji data assembly for the survey validation (732 observed tiles) and Guadalupe data assembly for the first two-month of the main surveys (653 observed tiles). Although the Early Data Release (EDR) only includes the Fuji data DESI Collaboration et al. (2023b), we consider both Fuji and Guadalupe data assemblies as the DESI early data in this paper. Figure 1 shows the distribution of the observed DESI tiles in the full DESI footprint.

The underlying stellar continuum is reproduced by using a full spectral fitting code of STARLIGHT Cid Fernandes et al. (2005), where the combination of 45 single stellar populations (SSPs) from Bruzual & Charlot (2003) and the Chabrier (2003) initial mass function (IMF) are adopted. The SSPs include 15 ages ranging from 1 Myr to 13 Gyr and three different metallicities (0.01, 0.02, 0.05). During the fitting, the main optical emission lines are masked out.

The modelled stellar continuum is subtracted from the observed spectrum and the fluxes of emission lines are then measured in the residual spectrum by Gaussian profile fitting. If not explicitly specified, all the spectral line measurements are calculated in rest frame. We are interested in the following emission lines: [O II]${\lambda\lambda}$3727,3729, H$\gamma$, [O III]$\lambda$4363, H$\beta$, [O III]$\lambda\lambda$4959,5007, [N II]${\lambda\lambda}$6548,6583, H$\alpha$, and [S II]${\lambda\lambda}$6716,6731. Due to their proximity, the following line groups are defined: [O II]${\lambda\lambda}$3727,3729; H$\gamma$ and [O III]$\lambda$4363; [O III]$\lambda\lambda$4959,5007; [N II]${\lambda\lambda}$6548,6583 and H$\alpha$; [S II]${\lambda\lambda}$6716,6731. Emission line fluxes in each group are fitted with Gaussian functions simultaneously in order to mitigate the contamination of neighboring lines. The steps of emission-line measurements are described below:

1. 1.

   The local background is estimated by a linear fit to the spectrum within a 70Å wavelength window around each line or line group, excluding regions affected by emission lines (15Å width). The linear fit is performed with a robust 3$\sigma$ clipping algorithm.

2. 2.

   The fitted local linear background is subtracted from the spectrum and the background error is estimated from the RMS of the linear fit.

3. 3.

   The background-subtracted emission line profiles (within a 15Å window) are modeled with Gaussian functions, where the mean, standard deviation and amplitude are free parameters to be fitted. The parameter uncertainties are provided by Gaussian fitting with the error spectrum as the weight.

4. 4.

   In addition to the total flux, Gaussian width, and line centers, we also calculate the rest-frame equivalent width for each emission line, whose error is estimated through error propagation.

Figure 2 shows a randomly selected DESI galaxy spectrum in rest-frame wavelength and the Gaussian fits to the emission-line profiles.

There are a total of 2,856,369 unique galaxies with successful measurements of spectral lines in the Fuji and Gaudalupe data sets. These objects are selected from the redshift catalogs with SPECTYPE $=$ “GALAXY” and ZWARN $=0$, which means that they are spectroscopically confirmed galaxies with reliable redshift measurements. Table 1 presents the sample selection steps and criteria. In order to eliminate spurious emission line detections due to cosmic rays and large random noise, we only use emission lines with Gaussian FWHM $>$ 1.2Å, which is empirically set to be slightly larger than the pixel scale of 0.8 Å in observed wavelength. To obtain credible metallicity measurements as shown in Section 3, we require that the [O III]$\lambda$4363, [O II]$\lambda$3726, [O II]$\lambda$3729, [O III]$\lambda$4959, and [O III]$\lambda$5007 flux measurements have high-enough SNRs as shown in Table 1. Due to the lower limit of SNR([O III]$\lambda$4363), the fitted FWHMs of this line might be abnormally large in some cases. In addition, the DESI data pipeline might mistakenly classify some broad-line AGNs as “GALAXY”. So we set an upper FWHM limit of 3 Å for [O III]$\lambda$4363 and H$\beta$, although it might exclude some galaxies with galactic outflow that would be worth further investigation in the future. The inclusion of [O III]$\lambda$5007 limits our galaxies to the redshift range of $z\lesssim 0.96$. We obtain a total of 1,633 [O III]$\lambda$4363-selected galaxies that satisfy all the above SNR and FWHM cuts.

The gas-phase extinction ($E(B-V)$) is calculated using the Balmer decrement, that is the change of the flux ratio of H$\alpha$/H$\beta$ or H$\beta$/H$\gamma$ relative to its intrinsic value. The intrinsic ratio (H$\alpha$/H$\beta$)₀ is 2.86 and (H$\beta$/H$\gamma$)₀ is 2.137, under the Case B recombination for $T_{e}=10^{4}$ K and electron density of $N_{e}=100$ cm^-3 (Hummer & Storey, 1987). Here, we impose 3$\sigma$ cuts on H$\beta$ and either H$\alpha$ or H$\gamma$ emission lines. We preferentially use the H$\alpha$/H$\beta$ to calculate the extinction if the H$\alpha$ line is not redshifted out of the wavelength coverage ($z\lesssim 0.49$) and meets the SNR requirement. As a common practice, the $E(B-V)$ is set to zero when the observed line ratios are lower than the intrinsic values (about 32% of the galaxies). It could be due to the uncertainties of the line ratio or statistical deviations from the theoretical values. The median reddening is about 0.08 mag for the [O III]$\lambda$4363-selected galaxies. All the emission-line fluxes are corrected with the extinction using the attenuation law of Cardelli et al. (1989).

Only metal-poor star-forming galaxies are the focus of this paper, so AGNs should be also excluded. Their optical spectra are dominated by non-thermal emission, making it complicated to study the stellar population. Some diagnostic diagrams with emission line ratios are usually used to discriminate different ionizing sources. Typically, the BPT diagrams (Baldwin et al., 1981) of [O III]/H$\beta$ vs. [N II]/H$\alpha$ and [O III]/H$\beta$ vs. [S II]/H$\alpha$ are the standard tools to separate star-forming galaxies from AGNs. Figure 3 shows the method of using the BPT diagrams to select star-forming galaxies. The data points in this figure are all the [O III]$\lambda$4363-detected galaxies. In the diagram of [O III]$\lambda$5007/H$\beta$ vs. [N II]$\lambda$6583/H$\alpha$ in Figure 3a, the empirical starburst limit as shown in solid line is adopted from Kauffmann et al. (2003):

where a modeling error of 0.1 dex (Kauffmann et al., 2003) is added to the original formula. The diagnostic diagram of [O III]/H$\beta$ vs. [S II]/H$\alpha$ in Figure 3b provides a supplementary classification, since the [N II]$\lambda$6583 line is not always detected and often blended with H$\alpha$. The separation curve is proposed by Kewley et al. (2001):

where we also add an additional error of 0.1 dex (Kewley et al., 2001). The original separation line as shown in dashed curve in Figure 3b is too close to the star-forming sequence. It would classify normal star-forming galaxies as AGNs. The above two BPT diagrams can only be used at redshifts $z<0.49$. At higher redshift, we attempt to use the diagnostic diagram of [O III]$\lambda$5007/H$\beta$ vs. [O II]/H$\beta$ first proposed by Lamareille et al. (2004) and revised in Lamareille (2010). Figure 3c shows the galaxy distribution on this diagram. The dashed line in this plot shows the updated demarcation curve of Lamareille (2010) to separate starbursts and AGNs. The solid line shows the separation limit plus an uncertainty of 0.15 dex. From this figure, we can see that a significant fraction of galaxies in the star-forming regions of the former two BPT diagrams lying above the curve, so the Lamareille (2010) separation is not optimal especially for galaxies with small [O II]/H$\beta$ and large [O III]$\lambda$5007/H$\beta$. Hence, we adopt the demarcation curve with the addition of 0.15-dex uncertainty (Lamareille, 2010), which is expressed as:

This new proposed separation limit is used to screen AGNs at higher redshifts in our sample via the the diagram of [O III]$\lambda$5007/H$\beta$ vs. [O II]/H$\beta$. We should note that the additional errors are added to both ordinate and abscissa in Figure 3 in order to set looser selection cuts of star-forming galaxies and make discrimination curves more suitable to our data. By using the above three BPT diagnostic diagrams, we get 1,623 star-forming galaxies. All the following analyses in this paper are based on this sample.

The gas-phase metallicity can be estimated using the calibrations developed from either theoretical models or empirical methods. The direct method to measure the metallicity is the so-called $T_{\rm e}$ method that relies on the determination of the electron temperature ($T_{\rm e}$) from auroral emission lines. $T_{\rm e}$ can be measured by using the line intensity ratio of two different ionization species of the same element, such as [O III]$\lambda\lambda$4959,5007/[O III]$\lambda$4363 (Aller, 1984). Although the direct $T_{e}$ method is subject to some caveats, such as metallicity underestimation caused by possible temperature fluctuation or gradients in high-metallicity star forming regions (Kewley & Ellison, 2008; Méndez-Delgado et al., 2023), it is currently the most reliable and accurate method to measure the metallicity.

In the direct $T_{e}$ method, the two-zone model is assumed, where the low and high ionization zones are traced by low and high ionization emission lines, respectively. The electron temperature in the high ionization zone $T_{e}$([O III]) can be estimated using the line ratio of [O III]$\lambda$4363 and [O III]$\lambda\lambda$4959,5997, given the electron density ($N_{\rm e}$) that can be calculated with the [O II] or [S II] doublets. The electron temperature in the low ionization zone $T_{e}([\mbox{{O~{}II}}])$ can be traced by [O II]$\lambda$7320,7330 lines or [N II]$\lambda$5755, which are too weak to be detected in the DESI spectra. An iterative method from Nicholls et al. (2014) is adopted to determine $T_{e}$([O II]) and the total oxygen abundance:

where $Z$ is the oxygen abundance. The iteration process starts with the O⁺⁺ abundance and the total oxygen abundance is the summation of the O⁺ and O⁺⁺ abundances determined by the two-zone electron temperatures. Usually, the calculation converges within several iterations.

The [S II] doublet (6716 and 6731 Å) can be detected for galaxies at $z<0.45$, while at higher redshift we can use the [O II]$\lambda\lambda$3726,3729 lines to calculate the electron density. We require that the flux SNRs of these lines are greater than 4. If the SNRs are not high enough or the line ratios are out of the valid ranges to determine $N_{\rm e}$, we assume a constant $N_{\rm e}$ $=100$ cm^-3 (Andrews & Martini, 2013; Ly et al., 2014; Gao et al., 2017). The electron temperature is insensitive to $N_{\rm e}$ in the relatively low density regime and assuming $N_{\rm e}$ $=100$ cm^-3 is consistent with the measurements from the [O II] or [S II] doublet. We use the Python package PYNEB²²2http://www.iac.es/protecto/PyNeb/ developed by Luridiana et al. (2015) to calculate the electron density, temperature, and oxygen abundance. As mentioned in Gao et al. (2017), the atomic recombination and collision strength data adopted in PYNEB are outdated, with the consequence that the electron temperature would be overestimated and thus the metallicity underestimated. We replace them with the updated data for O⁺, O⁺⁺, and S⁺ from Froese Fischer & Tachiev (2004); Kisielius et al. (2009); Tayal & Zatsarinny (2010); Storey et al. (2014). If the [S II] or [O II] doublets have enough SNR, the module getCrossTemDen in PYNEB is used to calculate $N_{\rm e}$ and $T_{\rm e}$([O III]) simultaneously. If $N_{\rm e}$ is set to 100 cm^-3, the PYNEB module getTemDen is used to calculate $T_{\rm e}$([O III]). The abundance of O⁺/H⁺ is determined by the PYNEB module getIonAbundance with [O III]$\lambda\lambda$4959,5007/H$\beta$ and $T_{\rm e}$([O III]). The O⁺⁺/H⁺ abundance is obtained by the same PYNEB module with [O II]$\lambda\lambda$3726,3929/H$\beta$ and $T_{\rm e}$([O II]) estimated by Equation (4). The total oxygen abundance 12+log(O/H) is calculated by adding the two abundances. The error estimation is implemented through Monte Carlo simulations. We add random errors to the fluxes of the measured emission lines according to their esitmated uncertainties. This process is repeated 1000 times. Each combination of the fluctuated fluxes is fed to PYNEB to calculate $N_{\rm e}$, $T_{\rm e}$, and 12+log(O/H). The final value of each parameter is calculated as the median value and the error is estimated from the 68% confidence interval of the parameter distribution.

Figure 4 shows the distributions of $N_{\rm e}$, $T_{\rm e}$([O III]) and 12+log(O/H) and their error distributions. There are a total of 702 and 692 galaxies with $N_{\rm e}$ calculated by [S II] and [O II], respectively. The remaining 229 galaxies have the assumed value of $\mbox{$N_{\rm e}$}=100$ cm^-3. About 90% of the galaxies have an electron density in the range of 76–629 cm^-3. The assumed $\mbox{$N_{\rm e}$}=100$ cm^-3 in this paper is in the valid range. There are about 1,608 galaxies with effective $T_{\rm e}$ measurement. The electron temperature of 90% of the galaxies ranges from 11,272 to 19,559 K and the median value is 14,421 K. The median error for $N_{\rm e}$ and $T_{\rm e}$([O III]) is 150 cm^-3 and 1,326 K. Although the relative error of $N_{\rm e}$ is large, it has little effect on the determination of the oxygen abundance due to its insensitivity to the electron density. The oxygen abundance of 90% of galaxies in our [O III]$\lambda$4363-selected galaxies ranges from 7.55 to 8.26 dex and the median value is 7.96 dex. The median metallicity error is 0.11 dex.

The stellar mass of each galaxy is estimated through stellar population synthesis fitting using the multi-wavelength photometric data, as well as the DESI optical spectrum. The photometric data come from the data release DR9 of the DESI legacy imaging surveys (hereafter, LS-DR9), which are mainly used for the target selections of the DESI spectroscopic survey (Zou et al., 2017; Dey et al., 2019). The imaging data in LS-DR9 provide both optical and infrared photometry, including 3 three optical bands ($g,r$, and $z$) and 4 infrared bands ($W1$, $W2$, $W3$, and $W4$) from WISE. Due to much shallower photometric depths and lower spatial resolutions of $W3$ and $W4$ bands, we only use the $W1$ and $W2$ to construct the photometric SED. The limiting magnitudes in $g$, $r$, $z$, $W1$, and $W2$ are 24.7, 23.9, 23.0, 20.7, and 20.0 mag, respectively (Zou et al., 2022). Although the WISE imaging is relatively shallow and has a low resolution of about 6″ (by contrast, 0.9″ for $z$ band), LS-DR9 provides consistent force-modeling photometry for both optical and infrared bands. Considering the importance of the near-infrared bands in estimating the stellar mass, we require that the two WISE bands have effective photometry with SNR $>2$.

The DESI optical spectra, on the other hand, cover the wavelength range of 3650–9800 Å, providing complementary information about the stellar population. However, the continuum SNR especially for ELGs is low, so that we can hardly constrain the stellar mass by only using full spectral fitting techniques. Instead, we obtain artificial broad-band spectrophotometry by convolving the optical spectra with self-defined 10 contiguous broad bands, which start from 3,650 Å and have bandwidths of about 615Å. Each self-defined filter is a boxcar filter with a flat transmission in a window of 615 Å width. Before convolution, we first scale each galaxy spectrum by multiplying by a scaling factor. The scaling factor is calculated as an aperture correction by comparing the total photometric fluxes in the $r$ band with the spectrophotometric flux in this filter. Here we assume that the light collected by the fiber is proportional to the total light of the galaxy and the color of each galaxy is the same inside and outside of the fiber. The 10 customized broad-band spectrophotometry as well as the 5 broad-band photometry, are used for determining the stellar mass. Both 5-band photometry and spectroscopic based photometry are used together to infer stellar masses. We should also note that in some cases the DESI targets are not selected from the DESI imaging surveys, in which case the spectra are scaled to the the Gaia $G$-band photometry.

We use the version 2022.1 of the Code Investigating GALaxy Emission (CIGALE; Boquien et al., 2019; Yang et al., 2020, 2022) to model the spectral energy distributions (SEDs) of galaxies. The CIGALE code can efficiently reconstruct the galaxy SED ranging from the far-ultraviolet to the infrared with highly flexible modeling of composite stellar population, gas emission, and dust attenuation and radiation. This modeling returns important physical properties of the galaxies, including stellar mass, star formation rate (SFR), dust extinction, and stellar age. In CIGALE, we adopt the simple stellar population (SSP) models of Bruzual & Charlot (2003), a Chabrier (2003) initial mass function (IMF), and a delayed star formation history (SFH). We select 21 stellar population ages (1 Myr – 13 Gyr), five stellar metallicities (0.0004, 0.004, 0.008, 0.02, 0.05), and 14 star formation timescales $\tau$ for the SFH ranging from 1 Myr to 80 Gyr. The Charlot & Fall (2000) extinction law and Draine et al. (2014) dust IR emission models are adopted to account for the dust attenuation of the interstellar medium and possible emission excess in WISE W1 and W2 bands, respectively. Due to strong emission lines in the spectra caused by the ionization of high-mass star radiation, the nebular emission is taken into account. However, no AGN emission is considered in the SED modeling, as we have already excluded those galaxies with AGN. The redshift is fixed to the spectroscopic redshift when constructing the models. Figure 5 shows two examples of the SED fitting with CIGALE to present two typical galaxies with old and young stellar populations, respectively.

To assess the quality of our stellar mass measurements, we cross-match the entire DESI galaxy sample with the COSMOS2020 catalog compiled by Weaver et al. (2022) for the Cosmic Evolution Survey (COSMOS). They have performed multi-wavelength photometry (35 bands) and calculated photometric redshifts for galaxies down to the depth of $i\sim 27$ over a 2 deg² footprint. The Farmer catalog and stellar mass derived with the LePhare SED fitting tool is used in this paper. There are more than 6300 galaxies observed by DESI located in the COSMOS field. We limit the galaxies for comparison with the photometric redshift accuracy ($|\Delta z|=|(z_{\mathrm{phot}}-z_{\mathrm{spec}})/(1+z_{\mathrm{spec}})|$, where $z_{\mathrm{phot}}$ and $z_{\mathrm{spec}}$ are photometric and spectroscopic redshifts, respectively) less than 0.02 and the logarithmic stellar mass uncertainty less 0.1 dex (artificially set to relatively smaller values for better comparison), generating a set of about 3600 galaxies. Figure 6 shows the comparison of the stellar mass calculated using the DESI and COSMOS data. It can be seen that the stellar mass measurements of these two surveys are consistent. There is a median bias of about 0.08 dex for $\Delta\log{M_{*}}=\log{M_{*}}[\mathrm{DESI}]-\log{M_{*}}[\mathrm{COSMOS}]$ . It is mainly caused by different configurations in the SED fitting codes for constructing galaxy spectral models (e.g. different IMFs, model spectral libraries, star formation histories). The dispersion $\sigma_{\Delta\log{M_{*}}}$ is about 0.16 dex, which is comparable to the estimated uncertainty with the DESI data (median of 0.18 dex). We also compare the stellar mass derived with only 5-band photometry and find that it has a similar bias but slightly larger dispersion of $\sigma_{\Delta\log{M_{*}}}\sim 0.17$ dex. Thus, the photometry plus spectrophotometry can indeed improve the stellar mass estimation. In particular, the spectrophotmetry is indispensable for those galaxies without DESI LS photometry.

The star formation rate is calculated from the de-reddened H$\alpha$ or H$\beta$ luminosity using the Kennicutt relation (Kennicutt, 1998):

where $F(\mbox{H$\alpha$})$ is the extinction-corrected H$\alpha$ flux, $d_{L}$ is the luminosity distance, and $L(\mbox{H$\alpha$})$ is the H$\alpha$ luminosity. We preferentially use the H$\alpha$ luminosity to estimate the SFR. At $z>0.49$ where the H$\alpha$ line is redshifted out of the wavelength coverage, the SFR from the dereddened H$\beta$ luminosity is calculated assuming the intrinsic flux ratio of $(\mbox{H$\alpha$}/\mbox{H$\beta$})_{0}=2.86$. Note that the emission-line luminosity is aperture-corrected using the same scaling factor as the one when we estimate the stellar mass through SED fitting. The SFR error is mainly from the calibration uncertainty of Equation (5) (about 30%). Note that the uncertainty from the aperture correction is not included.

From the oxygen abundance distribution in Figure 4, we can see that 902 [O III]$\lambda$4363-selected galaxies are metal-poor galaxies, which are defined to have 12+log(O/H) $<7.99$ ($Z<0.2\mbox{$Z_{\odot}$}$). No galaxy has an oxygen abundance larger than the solar metallicity. There are 223 XMPGs with 12+log(O/H) $<7.69$ ($Z<0.1\mbox{$Z_{\odot}$}$). These XMPGs account for 13.9% of all the [O III]$\lambda$4363-detected galaxies. They are about 0.01% of the total DESI observed galaxies, indicating their rarity. We define those galaxies with $Z+E(Z)<0.1\mbox{$Z_{\odot}$}$ as confirmed XMPGs and those with $Z<0.1\mbox{$Z_{\odot}$}$ but $Z+E(Z)\geq 0.1\mbox{$Z_{\odot}$}$ as XMPG candidates, where $E(Z)$ is the abundance uncertainty. There are 95 confirmed XMPGs, 13 and 2 of which have metallicity less than 1/20 and 1/30 $Z_{\odot}$, respectively. The remaining 134 galaxies are XMPG candidates. Only six of our XMPGs are previously identified, if a cross-matching radius of 2 arcsecs is considered.

Figure 7 shows the distributions of basic properties of all our XMPGs. About 74% of the XMPGs are located at relatively low redshift of $z<0.3$. There are 58 $z>0.3$ XMPGs with the highest redshift up to 0.96. This is by far the largest sample of high-redshift XMPGs. The magnitude histogram in Figure 7b presents a bimodal distribution, which corresponds to the two types of DESI spectroscopic targets (brighter for BGS and fainter for ELG). Our XMPGs are almost dust-free, and the median $E(B-V)$ is about 0.01 mag. The median $N_{e}$ and $T_{e}$ are 182 cm^-3 and 18,927 K, respectively. The median [O III]$\lambda$5007 equivalent width (EW([O III]$\lambda$5007)) is about about 528 Å. The median SFR, stellar mass ($M_{*}$), and specific star formation rate (sSFR $=$ SFR/$M_{*}$) are 0.65 $M_{\odot}$ yr^-1, 4.65$\times 10^{7}\mbox{$M_{\odot}$}$, and $1.14\times 10^{-8}$ yr^-1. About 89% of our XMPGs are dwarf galaxies with $M_{*}<10^{9}\mbox{$M_{\odot}$}$.

For comparison, we also present the parameter distributions of all [O III]$\lambda$4363-detected sample in Figure 7. We notice that the XMPGs and the parent sample have similar distributions of redshift, magnitude, extinction, electron density, and EW([O III]$\lambda$5007). The XMPGs have lower metallicities and hence high $T_{e}$. From the SFR and mass distributions, we find that the XMPGs tend to present relatively lower SFRs than the parent sample, mostly due to their lower mass. Further investigations show that the sSFR distributions of these two samples are more similar. The median sSFR of XMPGs is around 1.14 $\times 10^{-8}$ yr^-1, which is comparable to that of other [O III]$\lambda$4363-detected galaxies (1.34 $\times 10^{-8}$ yr^-1), but about 8 times higher than that of normal star-forming galaxies (the control sample selected in Section 4.4), whose median sSFR $\sim 1.36\times 10^{-9}$ yr^-1. Metal-poor galaxies are undergoing intense star formation. Table 2 and 3 in the appendix list the confirmed XMPGs and XMPG candidates identified in this work, respectively.

The two most metal-poor galaxies in our sample, with an oxygen abundance $Z/\mbox{$Z_{\odot}$}<1/30$, are: DESIJ150535.89+314639.4 at $z=0.054$ and DESIJ092331.28+645111.3 at $z=0.005$. Figure 8 shows their color images composed with the $g$, $r$, and $z$-band imaging. The DESI photometric reduction pipeline implements a forward modeling for each object with a point spread function (PSF) profile or a galaxy profile (Dey et al., 2019). The intrinsic shape parameters used in this paper are from the Legacy Survey DR9 (Schlegel et al., 2023). The model and residual images in Figure 8 show that the modeling works very well. The significant morphological discrepancy of these two XMPGs as shown in Figure 8 indicates possible different origins to their low metallicities. Figure 9 shows the observed spectra and photometric SEDs of these two XMPGs. Strong [O III]$\lambda$5007 and relatively weak [O II] emission lines in the spectra indicate that they have the high-ionization environments. Strong H$\alpha$ emission lines suggest intense star formation.

DESIJ150535.89+314639.4 This is the most metal-poor galaxy (DESI TARGETID $=$ 39628517750604871) identified by us. The oxygen abundance is 12+log(O/H) $=7.16\pm 0.05$ ($Z\sim 2.9\%\mbox{$Z_{\odot}$}$). The lowest metallicity reported so far is about 1.6%$Z_{\odot}$ (Kojima et al., 2020). DESIJ150535.89+314639.4 is located at the redshift of $z=0.054$. The $g$-band magnitude is 20.45 mag and the absolute magnitude is $M_{g}=-16.43$ mag. The optical colors are $(g-r)=-0.081$ and $(r-z)=-0.37$. The stellar mass and SFR are $1.51\times 10^{7}\mbox{$M_{\odot}$}$ and 0.22 $M_{\odot}$ yr^-1, respectively, giving a sSFR of $1.43\times 10^{-8}$ yr^-1. From Figure 8, we can see that DESIJ150535.89+314639.4 is a blue compact dwarf (BCD) galaxy with the half-light radius $r_{\mathrm{eff}}=640$ pc. The mean surface brightness $\mu_{g}$ within the half-light radius is about 21.22 mag arcsec^-2. It falls in the color cuts defined by Hsyu et al. (2018) to select BCD candidates and conforms to the BCD observational criteria listed in Morales-Luis et al. (2011).

DESIJ092331.28+645111.3 This is the second most metal-poor galaxy (DESI TARGETID $=$ 39633440282250119). The metallicity is 12+log(O/H) $=7.20\pm 0.09$ ($Z\sim 3.2\%\mbox{$Z_{\odot}$}$). It is located at $z=0.0054$ and $(\alpha,\delta)=(140.88032\arcdeg,64.85315\arcdeg)$. We can see from Figure 8 that DESIJ092331.28+645111.3 has a complex morphology, exhibiting a tadpole-like or cometary shape. Such a head-tail morphology is common among XMPGs as shown in Morales-Luis et al. (2011); Sánchez Almeida et al. (2016a); Kojima et al. (2020). The DESI fiber pointed to the head (i.e. the north-west part). The $g$-band apparent and absolute magnitudes of the head are $g=20.46$ and $M_{g}=-11.38$. The optical colors are $(g-r)=0.15$ and $(r-z)=-0.53$. The stellar mass and SFR are $1.38\times 10^{5}\mbox{$M_{\odot}$}$ and 0.0034 $M_{\odot}$ yr^-1, respectively, giving sSFR $=2.42\times 10^{-8}$ yr^-1. If we take the south-east part into account and assume the same mass-luminosity ratio (M/L) as the north-west one, the total stellar mass would be 1.67 times larger (may be larger due to possibly higher M/L for the tail). The sSFR is slightly larger than DESIJ150535.89+314639.4, suggesting that its mass assembly is faster although lower mass. Nevertheless, this extremely metal-poor head shows a blue color and compact size ($r_{\mathrm{eff}}=172$ pc). It connects to an extended fainter tail with a lower surface brightness. It was found that the XMPG heads and tails are dynamically connected and the tails have older stellar populations and are more metal-rich than the heads (Sánchez Almeida et al., 2015).

DESIJ150535.89+314639.4 and DESIJ150535.89+314639.4 present two different morphologies, suggesting that they might have distinct physical origins. It is also possible that DESIJ150535.89+314639.4 might be in the stage of bursty star formation, so that the low surface brightness tail is outshone by the bright starburst clump. Very faint structures can be seen around this galaxy. To check these, we need higher-resolution and deeper images. We will systematically investigate the morphologies and stellar population properties of our XMPG sample in the future.

The number of XMPGs identified here with DESI is almost double what is currently known. There are 338 XMPGs collected from the literature that are cross-matched with the DESI photometric catalogs (the photometric errors in $g$, $r$, and $z$ bands are limited to be less than 0.2 mag). Figure 10 compares the distributions of the DESI (this work) and literature XMPGs on the color-color diagram of $g-r$ vs. $r-z$. The literature XMPGs are collected by Sui et al. (in preparation). The magnitudes have been corrected for the Galactic extinction. In the right panel of this figure, we can see that most XMPGs in the literature seem to lie in a low-redshift sequence. Only a few XMPGs are identified at redshift of $z>0.3$. Our local XMPGs at $z<0.1$ also follow this sequence. A large sample of relatively higher-redshift XMPGs are identified to the right and the top of the sequence. This is mainly caused by the strong H$\alpha$ and [O III]$\lambda$5007 lines redshifting into or out of the $r$ and $z$ band, which is similar to the effect of green pea galaxies (Cardamone et al., 2009). These galaxies and their positions in the color space will help to identify XMPG candidates at $z<1$ from future large-scale imaging surveys.

Figure 11 displays the distribution of our XMPGs on the BPT diagram of [O III]$\lambda$5007/H$\beta$ vs. [O II]/H$\beta$. Other [O III]$\lambda$4363-detected galaxies are also overplotted in this figure. Almost all [O III]$\lambda$4363-detected galaxies except one galaxy are star-forming galaxies lying below the demarcation line. This galaxy is selected as a star-forming galaxy via the BPT diagrams of [O III]/H$\beta$ vs. [N II]/H$\alpha$ and [O III]/H$\beta$ vs. [S II]/H$\alpha$. Compared to other [O III]$\lambda$5007-detected galaxies, our XMPGs present slightly smaller [O III]/H$\beta$, but still have high ionization conditions due to large [O III]/[O II] ratios. We also check the loci of these galaxies in the Mass-Excitation (MEx) diagnostic diagram (Juneau et al., 2011, 2014) and find that only 2.5% of our XMPGs and 5.2% of [O III]$\lambda$4363-detected galaxies might host an AGNs. Different diagnostic diagrams may have slightly different AGN predictions. In addition, there are parameter measurement uncertainties in these diagrams. So there could be some inconsistent classifications using different diagnostic diagrams.

To characterize the environment of our XMPGs, we adopt the statistics adopted in the analyses of the so-called blueberry galaxies (Yang et al., 2017), which is the distance ($D$) of the nearest galaxy. We construct a control sample from the DESI spectroscopic catalog for comparison. The comparison galaxy samples are selected through closest matching in both spatial distance and stellar mass with the total spectroscopically-confirmed star-forming galaxies. This makes sure that the control sample has both the same redshift and mass distributions as those of our XMPGs. Figure 12 shows the cumulative distributions of $D$ for these two samples. From this figure, we can see that there is an obvious difference between the two distributions. The Kolmogorov-Smirnov test gives a p-value close to zero, proving the significance. We conclude that our XMPGs inhabit relatively low-density environments. Nevertheless, we do not rule out that XMPGs might be mergers or closely interact with neighbour satellites, because the statistics used in this paper is on a large scale.

As investigated by Cardamone et al. (2009) and Yang et al. (2017), both green pea and blueberry galaxies are strong [O III] and compact emitters at different redshifts. These galaxies are relatively metal-poor and are found to be situated in low-density environments. More than half of our XMPGs have effective radii less than 3 kpc and equivalent width of [O III]$\lambda$5007 larger than 200 Å. These features are used to select green pea-like galaxies. It could be foreseeable that our XMPGs are also located in low-density regions. In addition, Filho et al. (2015) and Sánchez Almeida et al. (2016a) explored the large-scale environment of XMPGs using the $N$-body cosmological numerical simulations. They found that XMPGs tend to reside in low-density environments similar to blue compact dwarf galaxies. About 75% of their XMPGs have a strong tendency to reside in voids and sheets.

It has been well established that the metallicity of star-forming galaxies correlates with the stellar mass (Tremonti et al., 2004), forming the so-called mass-metallicity relation (MZR). The MZR reveals that more massive galaxies appear to be more chemically enriched and low-mass galaxies for a steeper relation in this mass-metallicity plane. It was also found that the MZR evolves with redshift. The relation still holds at high redshifts out to 3 and beyond, while high-redshift galaxies are observed to be less enriched (Erb et al., 2006; Maiolino et al., 2008; Moustakas et al., 2011; Zahid et al., 2013; Ly et al., 2016). Analytic models and numerical simulations of galactic chemical evolution have been developed to interpret the MZR (Davé et al., 2011; Lilly et al., 2013; Torrey et al., 2018). These theoretical models disclose that the existence of MZR is the consequence of the secular interplay among different physical processes including stellar formation, gas outflow, and gas recycling and accretion.

Figure 13 presents the MZR of our XMPGs (colored points) and all [O III]$\lambda$4363-detected galaxies (grey points). It can be seen that our XMPGs at higher redshift tend to be the galaxies with higher stellar mass, mainly due to the observational selection effect. In much MZR research, the metallicity determinations are based on strong-line calibrations. However, such determinations generate nearly 1 dex differences depending on the choice of calibration method (Kewley & Ellison, 2008; Moustakas et al., 2010). For this reason, the MZRs based on the $T_{e}$ method are chosen for consistent comparisons. Andrews & Martini (2013) (hereafter AM13) obtained a MZR at $z\sim 0.08$ via the direct method for normal star-forming galaxies using stacked SDSS spectra. This MZR, spanning a wide stellar mass range of $\log{M_{*}}=$7.4–10.5, is plotted in Figure 13. We also display the MZR at $0.5<z<1$ derived by Ly et al. (2016), who used the direct $T_{e}$ method to measure the metallicity for several tens of galaxies at $z<1$ with Keck/MMT spectra. From this figure, we can see that [O III]$\lambda$4363-detected galaxies tend to lie below the local and Ly et al. (2016) MZRs, as expected since the [O III]$\lambda$4363 line is more easily detected in low-metallicity galaxies.

Our XMPGs cover a wide range of stellar masses, from $M_{*}=\sim 10^{5}$ to $\sim 10^{10}\mbox{$M_{\odot}$}$. About 44% of our XMPGs at $z<0.25$ are located below the AM13 MZR, and 100% of our XMPGs with $z>0.5$ are located below the Ly et al. (2016) MZR, taking the MZR uncertainties (about 0.2 dex) into account. The deviations of our XMPGs from the MZRs of normal star-forming galaxies indicate that they might have distinct evolutionary pathways or reside in different environments. In addition, we also present the MZR evolution from the simulations of Ma et al. (2016). Their MZRs at $z=0,1,3$, and 6 are shown in Figure 13. These simulations succeed in reproducing many observed properties (including MZR) of normal star forming galaxies at $0<z<6$. It can be clearly seen that a considerable number of our XMPGs are located within or out of the MZR at $z=$3–6. These XMPGs are probably the analogs of high-$z$ young galaxies, which are excellent samples for studying the evolution stage of the early universe.

In this analysis of the MZR relation, it is impossible to completely avoid sample selection effect. Our XMPG galaxy sample comes from the DESI early data, which includes the observations from the survey validation (SV) and first two months of the main surveys. The sample selection is extremely complicated because of different observing strategies during different SV stages and the main survey (DESI Collaboration et al., 2023a), different selection methods of galaxy targets, the incompleteness from the fiber assignment and redshift measurement, and the bias induced by the XMPG detection method. We use the [O III]$\lambda$4363 line to identify XMPGs in this paper. This emission line is very weak and tends to be detectable in metal-poor galaxies, so it bias the galaxies with low metallicity. Nevertheless, it is obvious that the XMPGs deviate from the normal MZR, illustrating their peculiarity. It is also possible that the lower MZR of our samples relative to normal galaxies is a physical trend. Juneau et al. (2014) mentioned that there might be a relation between the emission-line luminosity and the metallicity trends. The low-redshift galaxies with high emission-line luminosity (like our metal-poor galaxy samples) follow a systematically lower MZR, similar to those galaxies at high redshift. It has been also observed that the MZR has a secondary parameter dependency on the SFR (Ellison et al., 2008; Mannucci et al., 2010). There is a more fundamental mutual relationship between $M_{*}$, metallicity and SFR, that is so-called fundamental metallicity relation (FMR). It is found that there is no significant redshift evolution of FMR up to $z\sim 3$ (Mannucci et al., 2010). The existence of FMR means that galaxies with higher sSFR galaxies present lower metallicities.