Low-metallicity Galaxies from the Dark Energy Survey

We present a new selection of 358 blue compact dwarf galaxies (BCDs) from 5000 square degrees in the Dark Energy Survey, and the spectroscopic follow-up of a subsample of 68 objects. For the subsample of 34 objects with deep spectra, we measure the metallicity via the direct T e method using the auroral [O iii]λ 4363 emission line. These BCDs have an average oxygen abundance of 12+log(O/H) = 7.8, with stellar masses between 107 and 108 M ⊙ and specific star-formation rates between ∼10−9 and 10−7 yr−1. We compare the position of our BCDs with the mass–metallicity (M–Z) and luminosity–metallicity (L–Z) relation derived from the Local Volume Legacy sample. We find the scatter about the M–Z relation is smaller than the scatter about the L–Z relation. We identify a correlation between the offsets from the M–Z and L–Z relation that we suggest is due to the contribution of metal-poor inflows. Finally, we explore the validity of the mass–metallicity–SFR fundamental plane in the mass range probed by our galaxies. We find that BCDs with stellar masses smaller than 108 M ⊙ do not follow the extrapolation of the fundamental plane. This result suggests that mechanisms other than the balance between inflows and outflows may be at play in regulating the position of low-mass galaxies in the M–Z–SFR space.


INTRODUCTION
In the Λ cold dark matter (ΛCDM) scenario for structure formation, structures form hierarchically through mergers of smaller dark matter halos into larger ones.Gas accretes into the halos, fueling star formation Corresponding author: Yu-Heng Lin lin00025@umn.eduand forming the primeval, low-mass galaxies (Cole et al. 2000;Springel et al. 2005).These early, low mass galaxies serve as the building blocks of galaxies.They dominate the number counts at all epochs (McLeod et al. 2021).Nearby low mass galaxies offer a unique perspective to expand our understanding of galaxy formation.Observationally, the evolution of galaxies can be empirically described with the scaling relations between some basic physical properties such as the stellar mass, the star-formation rate (SFR) and the gas Lin et al. metal content.For example, Skillman et al. (1989); Pilyugin (2001); Lee et al. (2006); Guseva et al. (2009); Berg et al. (2012) find that galaxies from low to high metallicity obey a luminosity-metallicity relation (L-Z relation).Later studies focused on the tighter relation between the stellar mass and the oxygen abundance known as the mass-metallicity relation (M-Z relation; Tremonti et al. 2004;Berg et al. 2012;Izotov et al. 2015).
The M-Z relation shows that oxygen abundance increases with stellar mass from the low to intermediate mass systems (≤ 10 9.5 M ⊙ ), and then appears to flatten for more massive spiral galaxies (∼ 10 10.5 M ⊙ ) (Tremonti et al. 2004;Berg et al. 2012;Izotov et al. 2015).Some theoretical models of galaxy formation naturally explain the M-Z relation as the result of a balance between the production of metals in the star-formation episodes, and the dilution caused by inflows of pristine gas, and the removal of metals via outflows (Dayal et al. 2013;Ma et al. 2015;Spitoni et al. 2017).In these models, outflows in low mass galaxies are very efficient at removing metal-enriched gas because of the dwarfs' shallow gravitational potential wells.As the inflow fuels star formation, and star formation drives both the stellar mass growth and outflows, a close relation between the stellar mass, the metallicity, and the star formation rate is expected.Mannucci et al. (2010) have shown that for galaxies more massive than ≳ 10 9 M ⊙ these three properties are indeed correlated and this relation defines a surface in the 3D space.This relation, which is found to hold out to z ∼ 4, is referred to in the literature as the fundamental metallicity relation (FMR, e.g., Curti et al. 2020).At a given mass, the higher the SFR, the lower the metallicity.For galaxies with masses below 10 9 M ⊙ , the FMR is highly uncertain, due to the limited number of well studied low mass galaxies (Christensen et al. 2012;Belli et al. 2013;Telford et al. 2016).
For low mass galaxies, the L-Z relation and M-Z relation are intriguing at the low metallicity end.For a given oxygen abundance, many objects show brighter B band (or g band ) luminosities.In contrast, when compared with the more fundamental M-Z relation, the offsets are substantially reduced (McQuinn et al. 2020).The main reason for the offset from the L-Z relation is the higher SFR of the outliers, where young massive stars contribute a significant fraction of the luminosity (McQuinn et al. 2020).This effect is less pronounced in massive systems, or in the M-Z relation when the stellar masses are determined from data that sample the fainter low mass stars.Inflow of pristine gas or a minor merger can also affect the metallicity measurements, moving galaxies both in the L-Z and in the M-Z plane (Ekta & Chengalur 2010;Pascale et al. 2021).
Despite their importance, identifying low-metallicity dwarf galaxies is challenging, primarily because of their intrinsically-low luminosities and surface brightnesses (see discussion in Almeida et al. 2017).Consequently, these systems are more likely to be detected during a phase of active star formation.During this phase, young massive stars emit a strong continuum and ionize the surrounding interstellar medium (ISM), resulting in overall blue colors, and strong emission line spectra.Studies have successfully identified metal-poor systems using criteria based either on spectroscopic (Izotov et al. 2012;Guseva et al. 2017;Indahl et al. 2021;Hirschauer et al. 2022) or photometric diagnostics (Brown et al. 2007;James et al. 2015;Yang et al. 2017;Hsyu et al. 2018;Senchyna & Stark 2019;Kojima et al. 2020).In this work, we push the detection of blue compact dwarf galaxies (BCDs) to two magnitudes fainter than previous studies, using public images from the Dark Energy Survey (DES).For a sub-sample of 34 galaxies, we present spectroscopy to confirm the redshift and accurately measure nebular oxygen abundances using the direct T e method.
The structure of the paper is as follows.In Section 2, we describe the selection method of the BCDs.We report the spectroscopic follow-up observations and data reduction in Section 3. In Section 4, we describe our measurements.We discuss the results in section 5 and conclude in section 6.Throughout the paper, we assume ΛCDM cosmology with H 0 =67.66 km s −1 Mpc −1 (Planck Collaboration et al. 2020).

SAMPLE SELECTION
The Dark Energy Survey (DES, Abbott et al. 2018) is a broadband imaging survey covering 5,000 deg 2 of the southern sky.We use g, r, i, z catalogs and images released as part of the DES DR1 data release (Abbott et al. 2018), which includes data products derived from wide-area survey observations over 345 distinct nights between August 2013 to February 2016.
We selected objects with signal-to-noise ratio greater than five in the g, r, i bands and three in the z band, respectively, and with EXTENDED_COADD=3.The EXTENDED_COADD flag, explained in detail in Abbott et al. (2018), is used to separate galaxies from point-like stars and QSOs.This signal-to-noise ratio corresponds to a magnitude brighter than 23.7, 23.4, and 22.9, in the g, r, and i filters, respectively.For comparison, the imaging data of the Sloan Digital Sky Survey have a depth of 22.13, 22.70, and 22.20 in the corresponding filters (Padmanabhan et al. 2008).
To select BCD galaxies, we used similar color selection criteria as Hsyu et al. (2018).The criteria were optimized based on the photometric colors of known low metallicity systems such as I Zwicky 18 (Zwicky 1966;Skillman & Kennicutt 1993) and Leo P (Giovanelli et al. 2013;Skillman et al. 2013).The Hsyu et al. (2018) selection identifies blue compact sources with colors consistent with H ii regions.This is because by enforcing g band magnitudes brighter than z band, redder host galaxies are naturally excluded.
The color criteria are optimized to search for galaxies with colors consistent with low redshift (z ∼ 0.01 − 0.02) low-metallicity galaxies such as I Zwicky 18 and Leo P. Thus, we expect to select sources at low redshift (z ∼ 0.02), although by pushing to fainter magnitudes than the SDSS data we are able to search for fainter galaxies, which could be intrinsically faint or at higher redshift.We note that most of our spectra were observed with the Blue (b2k) disperser, which has the shortest observed wavelength of 3795 Å.To cover the [O ii]λ3727 emission line, we only obtain the spectroscopic follow-up for objects with redshift z ≥ 0.02 in this work.
The broadband color selection from the following query resulted in 895 galaxies from the full 5,000 square degrees area.The selection criteria we used are: where z err is the uncertainty of the z band magnitude.The complete SQL query was run on the NOAO datalab interface 1 , and can be found in Appendix A. The photometric errors were taken into consideration to maximize the purity of the selected sample.We further applied a Galactic extinction cut of E(B-V)< 0.066.
The original criteria from Hsyu et al. (2018), developed for SDSS, included an additional color condition based on the u band photometry, which is not available over the DES area.Therefore, we used data from the all-sky GALEX survey and transformed the u band condition from Hsyu et al. (2018) into a condition on the GALEX NUV filter (centered at 2,700 Å; Martin et al. 2005).We cross-matched the Hsyu et al. (2018) objects with the GALEX survey sources using a search aperture of 2.5 ′′ radius.We found that 60 out of the 88 Hsyu et al. (2018) objects also have GALEX data.We performed an orthogonal distance regression (taking 1 https://datalab.noirlab.edu/into consideration the errors in both bands) to find a linear relation between the magnitudes in the u and N U V filters.The best fit linear relation is: (2) Using this relation, we transform the Hsyu et al. ( 2018) u-band-based color criterion as a function of the N U V magnitude as: In Figure 1, we demonstrate this selection method using two stellar population synthesis models generated with Starburst99 (Leitherer et al. 1999).Specifically, we consider a low metallicity model (Z=0.002,12+log(O/H) ≃ 7.83) and a high (Z=0.014,12+log(O/H) ≃ 8.67) metallicity model using the "Geneva track without rotation" (Ekström et al. 2012), both 10 Myr old, redshifted to z = 0.05, and normalized to i magnitude = 20.Fluxes for the brightest recombination lines, including Hα, Hβ, Hγ, are calculated based on Case B line ratios and the intrinsic Hβ line luminosity produced by the model.The nebular emission lines are faint 10 Myr after the star formation.We assume the gas metallicity is the same as the stellar metallicity, and estimate the oxygen line ([O ii]λλ 3726, 3729 doublets and [O iii]λλ 4959, 5007 doublets) fluxes using the Maiolino et al. (2008) strong line calibration.The lower metallicity model (blue line) shows a bluer stellar continuum and brighter emission lines compared to the high metallicty model (red line), resulting in an overall bluer spectral energy distribution.As a result, our technique successfully identifies the low metallicity galaxies.
Out of the 895 galaxies selected on the basis of the optical colors alone, 516 objects had a corresponding GALEX counterpart.In the following we only consider galaxies that have a GALEX detection.We created a Zooniverse Panoptes project2 for the authors to visually inspect all galaxies with a GALEX counterpart.In order to exclude objects where the colors were due to image artifacts, we visually examined the individual band DES images (g, r, i) that entered the color selection.Each object was inspected by at least 3 group members, and we kept in our final clean catalog those 358 BCDs, for which all the experts agreed that the images were clean from artifacts.Figure 2 shows the distribution in the (NUV-g)-(g-r ) color-color space of the final DES BCD galaxy sample as well as the sample of Hsyu et al. (2018).Targets were chosen to be visible at airmass lower than 1.3 during the observations.To optimize the observing efficiency, we prioritized observing bright objects (g ≤ 20).For each object, we first obtained short exposures (200 seconds) optimized for the detection of bright lines and to confirm the galaxy's redshift.We then observed the confirmed galaxies with deeper exposure times aiming at the detection of the auroral [O iii] line required for the direct measurement of the gas-phase metallicity.For each object, we took three exposures to remove the cosmic rays.The total exposure time for each object ranged between 3×600 to 3×1200 seconds.
We used the Blue (b2k) disperser with a 0. ′′ 9 slit, covering the wavelength range between 3795 Å and 6615 Å, at a resolution of R=2400.We also observed 5 galaxies using the Wide (l2k) disperser with a 0. ′′ 9 slit, covering the wavelength range from 3570 Å to 7120 Å, at a resolution of R=1800.The seeing in the 2020 run varied from 0. ′′ 7 to 0. ′′ 9, with sporadic high cloud coverage.The seeing in the 2021 run varied from 0. ′′ 8 to 1. ′′ 4, with medium cloud coverage.For flux calibration, we observed the standard stars EG21 and LTT3864 in the 2020 observing run, and LTT4816 and LTT7379 in 2021 run (Stone & Baldwin 1983).
We used two different acquisition strategies in the two observing runs.In 2020, to properly place the target into the slit, we first centered on a bright star within two arcminutes from the target galaxy, and then the slit was rotated to include both the target and the acquisition star.Because the Blanco telescope does not have an atmospheric dispersion corrector (ADC), this strategy requires a correction in post processing for the wavelength-dependent flux-loss resulting from atmospheric refraction (see Section 3.2 and Appendix B).In 2021, we instead acquired directly on the targets (which therefore had to be slightly brighter than in the 2020 run) and oriented the slit at the parallactic angle to avoid slit-losses.During the 2021 observations, three objects selected from SDSS (labeled in Table 1) were observed as a filling for the first half of the night.The selection criteria of these objects are the same as described in Hsyu et al. (2018).

CTIO Data Reduction
We use the IRAF package (Tody 1986) to reduce and calibrate the spectroscopic data.The reduction and calibration processes include overscan and bias subtraction, flat-fielding, wavelength calibration, background subtraction, atmospheric extinction correction, aperture extraction, and flux calibration.The wavelengths are calibrated with a HgNe lamp.The aperture sizes in the cross-dispersion direction were chosen to optimize the S/N ratio of the continuum, after verifying that the emission lines showed the same spatial distribution as the continuum.
For the 2020 observation run, we correct the wavelength dependent slit-loss following Filippenko (1982, see details in Appendix B).We assume the source has a Gaussian spatial profile with standard deviation equal to half of the seeing full width half maximum (FWHM), and calculate the flux recovery factor ϵ(λ) = I(λ)/I ′ (λ) for both the standard star and targets, where I ′ (λ) is the observed amount of light, and I(λ) is the amount that goes through the slit if there were no refraction.See details in Appendix B.
The atmospheric refraction is most significant at the blue side of the spectra, which affects the flux measurements of [O ii]λ3727 and the derived 12+ log(O/H).The mean and maximum of the flux recovery factor ϵ(3727) are 1.36 and 1.66.Uncorrected, this flux loss would lead to a 0.13 and 0.22 dex lower values in the calculated oxygen abundances.In the 2021 observation run, we observed at the parallactic angle, so we do not apply a flux correction.
Figure 3 Using the short exposures, we confirm that 59 out of the 67 BCD candidates are emission line galaxies, yielding a ∼ 90% success rate for the chosen selection criteria.For the 35 galaxies with deep spectroscopic follow-up, we report the redshifts together with the Galactic extiction-corrected photometry in Table 1, and the 59 confirmed objects in the online material.For all 34 galaxies (14 from the 2020, and 20 from the 2021 run) we were able to detect the auroral [O iii]λ4363 line at a signal-to-noise ratio >5.

ANALYSIS
In this section, we first explain the measurements of the emission line fluxes (Section 4.1), and then the steps involved in the measurement of the galaxies' physical properties, including the gas oxygen abundance (Section 4.3), star formation rate (Section 4.5), and stellar mass (Section 4.4).

Emission line measurements
We estimate the redshift z using the brightest forbidden emission line [O iii] λ5007 Å.For each of the emission lines (i) we project the spectrum to velocity space, and measure the continuum level as the median of the flux density in the velocity range of v ± 500 km s −1 − 1000 km s −1 , after masking other emission lines within this range.The emission line flux is measured by direct integration of the spectral flux density over the continuum-subtracted line profile, typically within | v |≤ 200 km s −1 .The uncertainty of the emission line flux is calculated from the standard deviation of the continuum, adding in quadrature a 2% systematic flux uncertainty resulting from the flux calibration step.Emission line measurements for the sample galaxies are reported in Table 2.
We correct for Galactic extinction using the dust maps in Schlafly & Finkbeiner (2011) and the reddening curve of Cardelli et al. (1989) with R V = 3.1.

Electron temperature and internal dust extinction correction
In addition to Galactic extinction we also account for internal extinction, using the emissivity of the Balmer lines under the Case B assumption.Because the intrinsic Balmer line ratios depend mildly on the electron density, n e , and the electron temperature, T e , we proceed with the following iterative approach.We assume a fixed electron density n e =100 cm −3 throughout the calculation of the dust extinction correction.
We calculate the O ++ electron temperature T e (O iii) using the line ratio of [O iii]λλ4959,5007/[O iii]λ 4363 with the method described in Izotov et al. (2006), and the O + electron temperature T e (O ii) with the empirical relation (Pagel et al. 1992): where t e = T e × 10 −4 .The first estimate of the temperature is calculated without the dust reddening correction.Once we obtain the color excesses E(B−V), we correct the emission line fluxes and repeat the calculation described above until T e (O iii) changes by less than 1000K.We constrain our temperature range to be between 5000K ≤ T e (O iii) ≤ 22000K, since the H ii regions rarely exceed this temperature, and the overestimation of T e (O iii) leads to an underestimation of the oxygen abundance.We correct for the dust reddening and the underlying stellar absorption using the minimum χ 2 approach 3 : where X R (λ), σ 2 X obs are the observed Balmer line ratios and uncertainties over Hβ.X ext (λ) = F λ /Hβ is the product of extinction and intrinsic line ratio under Case B assumption.The extinction at a given wavelength λ is computed as the following: where F 0 λ is the flux without extinction, f (λ) is the reddening function from Cardelli et al. (1989) is the equivalent width of the line, and a H i (≤ 0) is the equivalent width of the underlying stellar absorption at Hβ; k(λ) is the scaling coefficient accounting for the wavelength dependence of the underlying stellar absorption (Aver et al. 2021).We use PyNeb (Luridiana et al. 2015) to calculate the theoretical line ratios under the Case B assumption, as a function of density and temperature.We include the Hα (if present), Hβ, Hγ, and Hδ in our calculation, while Hϵ and H8 are not included due to the blends of [Ne iii] and He i, respectively, and we do not have high S/N detections on H9 in many spectra.
The dust-extinction-corrected emission line fluxes are presented in Table 2, and the T e and derived oxygen abundances are presented in Table 3.The uncertainties of the dust extinction and the T e are computed via a Monte Carlo simulation.We generate 500 realizations by randomizing the Balmer line ratios within ±1σ.
For each realization, we compute the E(B-V) and the T e .We report the 16% and 84% percentiles of the distribution.

Metallicity
We use the total oxygen abundance relative to hydrogen as a proxy of the H ii region metallicity.We also assume that the oxygen total abundance can be written as: The oxygen abundance can be determined directly when the electron density, temperature, and line emissivities are known.This method is referred to as the direct method in the literature.When the gas metallicity is estimated empirically, using only the strongest emission lines, the measurements are said to use the strong line calibration.We summarize the direct method in Section 4.3.1 (see, e.g., Izotov et al. 2006).In Section 4.3.2,we also consider the alternative empirical method to estimate nebular oxygen abundance based on the detection of only strong emission lines calibrated (Maiolino et al. 2008;Jiang et al. 2019).We refer to the latter as the strong-line method.

Direct Method
The intensity of collisionally excited emission lines depends on the ion density, the electron density, and the electron temperature.The direct method measures the oxygen abundance relative to hydrogen based on the extinction-corrected, emission line ratios [O iii] where t 2 = t e (O ii) and t 3 = t e (O iii) are both in 104 K.

Strong Line Calibration
The direct method is an accurate method to determine the oxygen abundance.However, often one would have to use the empirical strong line calibration  Skillman (1989), McGaugh (1991), Pilyugin (2000), and van Zee et al. (2006) have demonstrated that the strong oxygen lines are dependent on both the oxygen abundance and ionization parameters in low metallicity H ii regions.The O 32 ratio is highly sensitive to the ionization parameter, so it can be used in combination with other ratios to break the degeneracy.Jiang et al. (2019) provide a calibration with a combination of R 23 and O 32 based on 835 star forming green-pea galaxies, with metallicities ranging from 7.2 ≤ 12+log(O/H) ≤ 8.6, and 15 galaxies with metallicities lower than 7.6.Since the metallicities of the BCDs fall into this range, we compare our results of the direct method with this strong line calibration.
The metallicity Z = 12+log(O/H) in this method is determined by solving the equation: To derive 12+log(O/H) as a function of R23 and y from equation (10), Jiang et al. (2019) suggest a different formula based on the values of y and log(R23).However, we found the result more reasonable if we only apply the following formula regardless of the value y: (11) We present the results and a comparison in Section 5.4.

Stellar Mass and Luminosity
We estimate the stellar mass of the galaxies using the photometric SED fitting code FAST++ (Kriek et al. 2009) 4 .The input data are the emission line corrected magnitudes in the F U V , N U V , g, r, i, and z bands.We calculate the emission line correction for the g-band magnitude using the available spectra.Specifically, we create a new spectrum where the strongest lines in the wavelength range covered by the filter are masked.We call f ′ (λ) the masked spectrum, and f(λ) the observed spectrum.We then compute the emission line corrected g band magnitude from the observed magnitude, g 0 , as: where P(λ) is the throughput of the filter.For the r band we proceed in a different way, because the spectra do not always cover the Hα wavelength range.Accordingly, we calculate the emission-line-free r band magnitude as follows: where f r0 is the flux density computed from the r band magnitude before emission-line correction, F(Hα) is the Hα emission line flux estimated from the intrinsic Hβ flux and the observed dust attenuation, and W r is the width of the r band filter.Finally, the GALEX F U V filter covers from 1340 Å to 1800 Å.Since all of the objects in our discussion have redshift z < 0.1, the F U V filter does not cover the Lyα emission, and we do not apply emission line correction on F U V .
In the SED fitting, we used the Flexible Stellar Population Synthesis (FSPS) model (Conroy & Gunn 2010), assuming a Chabrier initial mass function (IMF) (Chabrier 2003) and metallicity in the range of Z=0.0008 and Z=0.0031, corresponding to 5% to 20% solar metallicity.We describe the star formation history with a "tau model" (SFR(t) ∝ exp(−t/τ )), where t is the age of the galaxy and τ varies from 6.5 Myr to 10 Myr in steps of 0.5 Myr.Hsyu et al. (2018) estimate the stellar masses of the galaxies in their sample using a color-dependent mass-to-light ratio provided in Bell et al. (2003).However, this relation was derived from objects with masses larger than 10 8.5 M ⊙ , not optimized for low mass objects with starbursts, like those studied in this paper and in Hsyu et al (2018).For consistency, we calculated the stellar masses of BCDs in Hsyu et al. (2018) with the same SED fitting analysis described here.
Finally, to compare with other results in the literature, we need to derive the luminosity of the galaxies in the B−band filter.We used the empirical relation between the B band and the i−band and g − i color, given in Cook et al. (2014): where we used the g emission-line-free magnitude.

Star Formation Rate
We calculate the star formation rate with both the Hβ luminosity and UV continuum luminosity L ν (1500) using the method in Kennicutt (1998).Throughout the paper, we denote the SFR derived from Hβ flux as SFR and SFR derived from UV continuum as SFR(UV).The specific star formation rate (sSFR) is defined as sSFR = SFR/M * .The reported Hβ luminosities, L(Hβ), are calculated using the slit-loss corrected Hβ line flux.We estimate the slit-loss corrected flux by assuming a Gaussian profile with the width similar to the seeing, and calculate the slit loss for the assumed profile observed with a 0.9 ′′ -wide long slit.We estimate the SFR with the extinction-corrected Hβ luminosity adopting an intrinsic Hα/Hβ ratio of 2.86: where L(Hβ) is in units of erg s −1 and the SFR is in M ⊙ yr −1 .ϵ(IMF) is the factor to correct the flattening of the stellar IMF below 1 M ⊙ for a Chabrier IMF, instead of the power law Salpeter IMF (Salpeter 1955) adopted by Kennicutt (1998): where ξ Ch (m) and ξ Sal (m) are the Chabrier IMF and Salpeter IMF, respectively.The SFR of our samples ranges from 0.03 to 3 M ⊙ yr −1 .At these low levels of SFR one would normally have to worry about the sampling of the high mass end of the IMF.However, the presence of strong [O iii] emission lines shows that the most massive stars are still alive in the current bursts.We calculate the UV continuum L ν (1500 Å) from the flux density of GALEX F U V magnitudes, and correct the dust extinction.The SFR(UV) is calculated as (McQuinn et al. 2015a): SFR(UV) = ϵ(IMF) × (2.04 ± 0.81) × 10 −28 L ν (1500 Å).
(17) The SFR(UV) of our samples ranges from 0.02 to 10 M ⊙ yr −1 .The uncertainties of the SFR(UV) are computed by generating 500 realizations by randomizing the E(B-V) and F U V magnitudes within ±1σ, and combined with the uncertainty in equation 17.For each realization, we compute the SFR(UV).We report the 50% percentiles as the SFR(UV), and 16% and 84% percentiles of the distribution.
Lin et al.

DISCUSSION
In this section, we compare the galaxy properties of the selected BCDs in the DES with samples of similar galaxies from the literature.Specifically, we consider the sub-sample of BCDs in Hsyu et al. (2018) for which the metallicity is measured via the direct method, the LVL galaxies in Berg et al. (2012) with their SFR reported in Lee et al. (2009), and the extremely-metal poor (XMP) galaxies with 12+log(O/H) ≤ 7.35 (1/20 Z ⊙ ), defined in McQuinn et al. (2020).The galaxies in these samples are all at low redshift (<0.1) with 12+log(O/H)< 8.25 (0.4 Z ⊙ ) determined with the direct method.

Specific Star Formation Rate
Figure 4 shows the SFR(Hβ) and stellar mass of the DES BCDs (blue stars), those identified in Hsyu et al. (2018, green squares), as well as the LVL galaxies from Lee et al. (2009, black dots) and the XMP galaxies in McQuinn et al. (2020).For reference, we also show the   detected in the GALEX N U V band.For the same N U V magnitudes, higher redshifts imply higher SFRs.
(19) We show how our galaxies compare with the L-Z and M-Z relations in Figure 6 and Figure 7, respectively.
In Figure 6, most of our BCDs are brighter than the B−band luminosity predicted from their metallicity and lie in between the LVL galaxies and the XMP galaxies.This is not surprising considering that the SFRs of the BCDs are biased high.The high SFR indicates that more young, massive, bright stars are currently present in the galaxy.These massive stars contribute a significant luminosity to the B−band, making the BCDs brighter than the L-Z relation.In terms of stellar mass, massive stars are not as impactful as in luminosity.For example, an O-star can be as bright as 10 5 L ⊙ , with only 20 M ⊙ .Therefore, the M-Z relation is not strongly biased by recent star formation and shows a tighter relation among the low mass galaxies.As shown in Figure 7, most of the BCDs are distributed within 1σ from the M-Z relation derived by the LVL galaxies, and many of the XMP galaxies that were outliers in the L-Z relation are also now close to the M-Z trend line.
Figure 7 shows that even in the M-Z relation, some low metallicity outliers remain among the XMP galaxies selected by McQuinn et al. (2020).One possibility to explain the objects that remain outliers in both the L-Z and M-Z relations is that the measured metallicity is not representative of the stellar mass of the underlying galaxy.This lower metallicity is likely attributed to external events.For example, this can happen if the H ii region results from SF in a pocket of pristine gas acquired from the circumgalactic medium (Ekta & Chengalur 2010), or if it is the result of a minor merger with a galaxy/component that has a lower metallicity (Pascale et al. 2021).In these scenarios, the metallicity of the main galaxy would be underestimated, and the object would be an outlier in both L-Z and M-Z relations.
In an attempt to identify galaxies with possible low-metallicity inflows, we show in Figure 8 the luminosity and stellar mass offsets of the galaxies, and we define the offsets as the horizontal distance of the data from the M-Z and L-Z relations defined by the LVL galaxies.Most of the DES BCDs and the Hsyu et al.'s sample have small scatter from the M-Z relation, with ∆M * (Z)<1 σ of the M-Z relaion, and larger scatter from the L-Z relation, with ∆L(Z)∼2 σ of the L-Z relation.For the XMP galaxies there appears to be a correlation between ∆L(Z) and ∆M * (Z), with most of the XMP galaxies located on the upper right side of the plot.This result suggests that the "extreme metal-poor" nature of these objects is a consequence of pristine gas inflow.Although not as significant as for the XMP galaxies (which are extreme by definition), the offsets for the BCDs appear to correlate in a similar way to the XMP galaxies.This suggests that pristine inflows may in fact be very common, as required by galaxy models (Davé et al. 2012;Ma et al. 2015).

Mass-Metallicity-Star Formation Rate Relation
The fundamental metallicity relation (FMR) describes a second order dependence of the M-Z relation on the SFR.Specifically, it is observed that for a given stellar mass, the lower the metallicity, the higher the star formation rate of the galaxy.Many theoretical models reproduce the observed relation as the result of the interplay among various effects: pristine gas brought in with the inflow fueling star-formation, enriched elements being expelled by outflows, and the varying amount of gas available for star formation.The dependence on the SFR is likely a reflection of the dependence on the gas fraction and the efficiency of the baryonic cycle (Bothwell et al. 2016b,a).
Note that the SFR dependence on the fundamental metallicity relation is different from how the SFR affects the L-Z relation.The latter is biased by the SFR since the relation is determined from the luminosity of a single broad band filter.
The FMR in Curti et al. (2020) is derived using SDSS galaxies with stellar mass 8 < log(M * ) < 12, given as agrees well with the M-Z relation found with the LVL galaxies Berg et al. (2012).We compared the FMR derived for the low mass objects with our DES galaxies, but find no significant second-order dependence of the metallicity on the SFR.Objects with 12-log(O/H) < 8.0 do not fall on the FMR surface described in equation ( 20).When calculating the metallicity expected from the FMR relation for the stellar mass and SFR of each galaxy, we find that the metallicity residuals ∆Z = Z FMR − Z dir correlates with the SFR.This is shown in Figure 9, where Z dir is the metallicity measured in Section 4.3.1 using the direct method.The mean of the metallicity residuals ∆Z is 0.45 dex, greater than the uncertainty of ∼0.25 at the low mass end of FMR (see Curti et al. 2020, Figure 11).Intriguingly, the metallicity residuals are actually larger for the lower SFR galaxies, while these lower SFR galaxies have the sSFR closer to the main-sequence galaxies.We fit a linear correlation to the residual ∆Z as a function of the Log(SFR): ∆Z = (−0.14± 0.02) × log(SFR) + 0.23 ± 0.04.( 22  In the updated equation 23 the SFR dependence is weak enough that it can be neglected, and the resulting the M-Z relation is almost identical to the one of Berg et al. (2012).Therefore the dispersion in the M-Z relation of the BCDs cannot be tightened by the star formation rate, and other factors may be driving it.
The FMR plane is the consequence of "main sequence" galaxies reaching equilibrium between gas inflow, outflow, and star formation.The extrapolation of the FMR plane cannot be taken literally, but what is the applicable lowest mass for this equilibrium model?With bursty star formation histories and lower gravity, dwarf galaxies such as our BCDs can be easily driven away from the equilibrium between gas flows and star formation.However, regardless of the burstiness, the LVL galaxies, having similar sSFR to those of the "main sequence" galaxies, do not follow the FMR plane either.Llerena et al. (2023) also find that lower sSFR galaxies have greater metallicity residuals from a z ∼ 3 SFG sample.Detailed mechanisms of the balance between inflows and outflows may be at play in regulating the position of low mass galaxies in the M * -SFR-Z space.For example, Magrini et al. (2012) considers a multi-phase interstellar medium model, and shows that compact regions with dense gas have distinct properties from a diffuse region.The scatter of the M-Z relation may result from galaxies with different gas density.Forbes et al. (2014) shows that a large intrinsic scatter of the mass loading factor (the ratio between the ejected mass from galatic winds and the galaxies' SFR) will weaken the anticorrelation between metallicity and SFR, at a fixed mass.The variation among mass loading factors in low mass galaxies can be enhanced by galactic winds generated from supernova explosions.These galactic winds are metal-enriched (Vader 1986(Vader , 1987)), and may escape (blowout) or even remove the ISM from the galactic halo (blow-away) given their shallow gravitational potential (Mac Low & Ferrara 1999).For example, McQuinn et al. (2015b) finds the low metallicity dwarf galaxy Leo P held on to only 5% of the oxygen atoms, the other 95% was presumably lost to the galaxy halo or the intergalactic medium.The oxygen retained in Leo P is much lower than estimates of 20%-25% of metals retained in more massive galaxies (10 9 M ⊙ < M * < 10 11.5 M ⊙ , Peeples et al. 2014).

Strong Line Calibration
Here we explore the reliability of strong line calibrations in Section 4.3.2 with the BCDs presented in this work.Note that while there are many strong line calibrations (e.g Maiolino et al. 2008;Izotov et al. 2021), we choose Jiang et al. (2019) specifically because it is calibrated on green pea galaxies and is the most directly applicable calibration available for the metallicities of our sample.Figure 10 shows the comparison of the metallicity derived from the direct method (Z dir ) and the strong line calibration (Z SL ) in Jiang et al. (2019).For most objects, the metallicity calculated from the strong line method is slightly lower, within 0.2 dex, compared to the values computed from the direct method.The mean offset is ∆Z = Z SL -Z dir = 0.10 ± 0.16 dex.The second panel of Figure 10 shows ∆Z versus the O 32 ratio.The third panel of Figure 10 shows ∆Z versus T e , and the fourth panel shows ∆Z versus the metallicity offset in the M-Z relation.
As for the photoionization-parameter-sensitive line ratio O 32 , the scatter ∆Z is greater when Log(O 32 ) < 0.5, consistent with the result in Jiang et al. (2019).The observed wavelength of the [O ii]λ 3727 emission line is at the blue end of our spectra, where the uncertainty of the flux calibration is greater.An overestimation on [O ii]λ 3727 would lead to lower Log(O 32 ) and higher metallicity Z dir .A low O 32 line ratio, if measured correctly, represents a lower temperature and a fainter [O iii]λ 4363 line, which increases the uncertainty in the direct method.This is consistent with the tighter correlation between the ∆Z and the electron temperature T e , as T e is determined without the [O ii] emission.The metallicity of the outliers in the M-Z relation is caused by the dilution of the infall pristine gas rather than the whole galaxy.This dilution does not affect the estimated metallicity.As expected, we do not observe a significant correlation between the ∆Z and the metallicity offset from the M-Z relation (∆Z(M * ) = Z dir -Z(M * )).

CONCLUSIONS
We present a new selection of 358 blue compact dwarfs from 5,000 deg 2 in the DES imaging survey, and the spectroscopic confirmation of a sub-sample of 64 galaxies using the CTIO COSMOS 4m telescope.The mean redshift of the confirmed BCDs is 0.04.
For the spectroscopic subsample, we estimate the galaxies' luminosities, stellar masses, and star formation rates from the combination of imaging and spectroscopic data.
For 34 galaxies with deep spectroscopic follow-up we measure the metallicity via the direct T e method using the auroral [O iii]λ 4363 emission line.The metallicity of our objects ranges from 7.37 ≤ 12+log(O/H) ≤ 8.16, with a mean of 12+log(O/H)= 7.8.The SFRs are in the range of 0.03 to 3 M ⊙ yr −1 , with stellar masses between 10 7 to 10 8 M ⊙ .The sSFR of the BCDs are around 10 −9 to 10 −7 yr −1 .
We compare the position of our BCDs with the M-Z and L-Z relation derived from the Local Volume Legacy sample.We find that in the L-Z plane, the DES BCDs fall below the LVL L-Z relation but are above the most extreme metal-poor galaxies.In contrast, in the M-Z plane, the DES BCDs agree with the local M-Z relation, and the scatter around the M-Z relation is smaller than the scatter around the L-Z relation.We identify a correlation between the offsets from the M-Z and L-Z relation that we suggest is due to the contribution of metal-poor inflows.
Finally, we explore the validity of the mass-metallicity-SFR fundamental plane in the mass range probed by our galaxies (< 10 8 M ⊙ ).This is the largest sample of galaxies in this mass range with accurate metallicity that can be used for such a study.Our results show that low mass star forming galaxies do not follow the extrapolation of the fundamental plane.This result suggests that mechanisms other than the balance between inflows and outflows may be at play in regulating the position of low mass galaxies in the M-SFR-Z space.Recent JWST spectroscopic results show that z > 7 galaxies behave in a similar way as our sample of BCDs: they are consistent with the z = 0 M-Z relation, while deviating from the mass-metallicity-SFR fundamental plane (Curti et al. 2022).Understanding this common behavior at the lowest masses will shed light on the metal enrichment of galaxies in the early universe.A. SQL QUERY SELECT coadd object id, ra, dec, snr g, snr r, snr i, snr z, mag auto g, mag auto r, mag auto i, mag auto z, magerr auto g, magerr auto r, magerr auto i, magerr auto z, ebv sfd98 FROM des dr1.galaxiesWHERE snr g > 5 AND snr r > 5 AND snr i > 5 AND snr z > 3 AND (mag auto g -mag auto r) BETWEEN -0.2 and 0.2 AND (mag auto r -mag auto i) BETWEEN -0.7 and -0.1 AND (mag auto i -mag auto z) < 0.1 AND (mag auto i -mag auto z) > -0.4 -(2*magerr auto z) AND mag auto g > 22 AND mag auto r > 22 AND flags g = 0 AND flags r = 0 AND flags i = 0

B. SLIT LOSS CORRECTION
In this section we describe the dispersion correction we apply on the 2020 observation run.Following the equation ( 1) and (2) in Filippenko (1982), we first calculate the refractive index n(λ) T,P , with temperature T ≃ 9 • C and P ≃ 600 mm Hg, according to our observation condition.The atmospheric differential refraction relative to λ = 5000 Å is ∆R(λ) ≡ R(λ) − R(5000 Å) ≃ 20625(n(λ) − n(5000 Å))tanz, (B1) where z is the zenith angle.We assume the target point spread function as a Gaussian profile, with the seeing being the diameter (2σ) of the profile, µ(x, y) = 1 2πσ 2 e −(x 2 +y 2 )/2σ 2 .(B2) The amount of light that goes into the slit with slit width 2a and slit length 2b arc second is

Figure 1 .
Figure 1.The SED and the spectra of two stellar population models at 10 Myr.The lower metallicity model is shown in blue and the higher metallicity model is shown in red.Both spectra are normalized at i magnitude=20.The throughput of the g, r, i, z filters are shown in gray.The lightblue shaded areas indicate our selection criteria.

Figure 2 .
Figure 2. The 358 sources chosen in this paper (gray cross) as well as Hsyu et al. (2018) BCDs (green square).The blue stars represent sources with follow-up spectroscopic confirmation.
(a) shows an example of the b2k spectrum taken in 2020, with the insets zooming-in on the [O ii]λ 3727 doublet and the temperature sensitive [O iii]λ4363 emission line.Figure 3(b) and Figure 3(c) shows the b2k and l2k spectra taken in 2021.

Figure 3 .
Figure 3. Three example spectra of BCDs.(a)The b2k spectrum observed in Jan., 2020.The intensity before (black) and after (blue) the atmospheric refraction correction are shown in the inset zooming-in on the [O ii]λ3727 and [O iii]λ4363 lines.(b) The b2k spectrum observed in June, 2021.(c) The l2k spectrum observed in June, 2021.

Figure 4 .
Figure 4.The mass vs. Log(SFR(Hβ)).Blue stars are the DES BCDs observed in the work, the green squares are the BCDs reported in Hsyu et al. (2018), the black triangles are the LVL galaxies in Berg et al. (2012), and the open squares are the XMP galaxies in McQuinn et al. (2020).The black contours are typical SFG derived from the SDSS DR8 catalog.sSFRs are indicated with the dotted lines.

Figure 5 .
Figure 5. Number counts of g-band from the 358 DES BCD candidates (gray), BCDs selected in Hsyu et al. (2018), and the 35 BCDs presented in this work.

Figure 6 .
Figure 6.Absolute B band Magnitude vs. Metallicity of our spectroscopic sub-sample (blue stars), BCDs selected in Hsyu et al. (2018) (green squares), Local Volume Legacy galaxies selected in Berg et al. (2012) (black triangle dots), and the XMP galaxies (McQuinn et al. 2020) (open squares).The 1σ scatter of the L-Z relation (equation 18) is shown in gray shaded areas.

Figure 7 .
Figure 7. Stellar mass vs. Metallicity of the same galaxy samples as in Figure 6.The 1σ scatter of the M-Z relation (equation 19) is shown in gray shaded areas.Note the tighness of this relationship compared to that in Figure 6.

Figure 8 .
Figure 8.The horizontal offsets of the L-Z relation vs. the M-Z relation.The gray shaded areas are the 1 σ scatter of the L-Z and M-Z relation.

Figure 10 .
Figure 10.First panel: 12 + log(O/H) calculated via the direct (Z dir ) vs. Strong Line Calibration (ZSL) in Jiang et al. (2019).The dashed line is the identical line.The light blue shadows are ±0.1 dex as uncertainty.Second panel: The metallicity difference of the two methods ∆Z = ZSL-Z dir vs. Log(O32).The light blue shadow indicates ±0.16 dex as uncertainty in the strong line calibration as per Jiang et al. (2019).Third panel: ∆Z vs. Te/10 4 K. Fourth panel: ∆Z vs. Metallicity offset in the M-Z relation.

I
(λ) moves the light profile away from the center of the slit,∆x ≃ ∆Rsin(θ); ∆y ≃ ∆Rcos(θ), (B4)where θ is the difference between the slit position angle and the parallactic angle.With the refraction, the amount of light that goes into the slit is ′ (λ) is the observed flux density, and I(λ) is the flux density that result to intrinsic emission line ratios.We define the flux recovery factor ϵ(λ) as ϵ(λ) ≡ I(λ)/I ′ (λ), and multiple the factors along the spectra.

Table 1 .
The coordinates, the measured spectroscopic redshift, the magnitudes in FUV, NUV, g, r, i, z bands of the BCDs presented in this work.The magnitudes are corrected for Galactic extinction.The full sample with spectroscopic confirmed redshift are included in the online material.a: The last three objects are selected from SDSS.

Table 2 .
The corrected lines fluxes of the observed BCDs, normalized to Hβ.The Hβ flux is reported in the unit of 10 −15 erg s −1 cm −2 .The 3σ upper limit is reported for low S/N emission lines.

Table
The derived properties for the observed BCDs.For object J2120-5054, the [O iii]λ4363 has S/N < 5 ,therefore the Te[O iii] are not reported.