# The Direct-Method Oxygen Abundance of Typical Dwarf Galaxies at Cosmic High-
Noon111The data presented herein were obtained at the W. M. Keck Observatory,
which is operated as a scientific partnership among the California Institute
of Technology, the University of California, and the National Aeronautics and
Space Administration. The Observatory was made possible by the generous
financial support of the W. M. Keck Foundation.
Based on observations made with the NASA/ESA Hubble Space Telescope, obtained
from the Data Archive at the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
under NASA contract NAS5-26555. These observations are associated with
programs #9289, #11710, #11802, #12201, #12931, #13389, #14209.

Timothy Gburek Department of Physics & Astronomy, University of California,
Riverside, CA 92521, USA  Brian Siana Department of Physics & Astronomy,
University of California, Riverside, CA 92521, USA  Anahita Alavi IPAC,
California Institute of Technology, 1200 E. California Boulevard, Pasadena, CA
91125, USA  Najmeh Emami Minnesota Institute for Astrophysics, University of
Minnesota, Minneapolis, MN, 55455, USA  Johan Richard Univ Lyon, Univ Lyon1,
Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230,
Saint-Genis-Laval, France  William R. Freeman Department of Physics &
Astronomy, University of California, Riverside, CA 92521, USA  Daniel P. Stark
Steward Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ
85721, USA  Christopher Snapp-Kolas  Department of Physics & Astronomy,
University of California, Riverside, CA 92521, USA

###### Abstract

We present a Keck/MOSFIRE, rest-optical, composite spectrum of 16 typical,
gravitationally-lensed, star-forming, dwarf galaxies at $1.7\lesssim z\lesssim
2.6$ ($z_{\rm{mean}}=2.30$), all chosen independent of emission-line strength.
These galaxies have a median stellar mass of
log($M_{\ast}$/$\rm{M_{\odot}}$)med = 8.29 and a median star formation rate of
$\rm{SFR_{H\alpha}^{med}=2.25\ M_{\odot}\ yr^{-1}}$. We measure the faint,
electron-temperature-sensitive, [O III] $\lambda$4363 emission line at
$2.5\sigma$ ($4.1\sigma$) significance when considering a bootstrapped
(statistical-only) uncertainty spectrum. This yields a direct-method oxygen
abundance of $12+\log(\rm{O/H})_{\rm{direct}}=7.87^{+0.24}_{-0.22}$
($0.15^{+0.11}_{-0.06}\ \rm{Z_{\odot}}$). We investigate the applicability at
high-$z$ of locally-calibrated, oxygen-based, strong-line metallicity
relations, finding that the local reference calibrations of Bian et al. (2018)
best reproduce ($\lesssim 0.15$ dex) our composite metallicity at fixed
strong-line ratio. At fixed $M_{\ast}$, our composite is well-represented by
the $z\sim 2.3$ direct-method stellar mass$\,-\,$gas-phase metallicity
relation (MZR) of Sanders et al. (2020). When comparing to predicted MZRs from
the IllustrisTNG and FIRE simulations, having recalculated our stellar masses
with more realistic non-parametric star formation histories
(log($M_{\ast}$/$\rm{M_{\odot}}$)med = 8.92), we find excellent agreement with
the FIRE MZR. Our composite is consistent with no metallicity evolution, at
fixed $M_{\ast}$ and SFR, of the locally-defined fundamental metallicity
relation. We measure the doublet ratio [O II] $\lambda$3729/[O II] $\lambda
3726=1.56\pm 0.32$ ($1.51\pm 0.12$) and a corresponding electron density of
$n_{e}=1^{+152}_{-0}\ \rm{cm^{-3}}$ ($n_{e}=1^{+49}_{-0}\ \rm{cm^{-3}}$) when
considering the bootstrapped (statistical-only) error spectrum. This result
suggests that lower-mass galaxies have lower densities than higher-mass
galaxies at at $z\sim 2$.

galaxies: abundances - galaxies: dwarf - galaxies: evolution - galaxies: high-
redshift - galaxies: ISM

††facilities: Keck:I (MOSFIRE), HST (WFC3, ACS)

##  1 Introduction

The gas-phase metallicity, or gas-phase oxygen abundance ($12+\rm{log(O/H)}$)
of the interstellar medium (ISM) of galaxies, is a cornerstone in the study of
galaxy formation and evolution. The metallicity traces the stellar mass
buildup of galaxies through the enrichment over time of the ISM by heavy
elements produced via stellar nucleosynthesis. Galaxies, however, are not
closed boxes and have inflows of metal-poor gas$\,-\,$the fuel for star
formation$\,-\,$ from the circumgalactic medium (CGM) and intergalactic medium
(IGM) as well as outflows of metal-laden gas triggered by feedback from
supernovae and/or active galactic nuclei (AGN). This modulation of the
enrichment of the ISM, via gas flows, shows metallicity to also be an
important physical property in the study of the cycle of baryons into, out of,
and within (baryon recycling) galaxies.

The combination of these processes is reflected in the scaling relation of the
gas-phase metallicity ($Z$) with the stellar mass ($M_{\ast}$) of star-forming
galaxies, known more succinctly as the mass$-$metallicity relation or MZR.
This relation demonstrates that stellar mass and metallicity are positively
and tightly correlated, whereby in the local Universe, below a characteristic
mass of $M_{\ast}\approx 10^{10.0-10.5}\ \rm{M_{\odot}}$, the MZR is described
by a power-law. Above this characteristic mass, or “turnover” mass, the MZR
flattens and asymptotically approaches an upper-limit oxygen abundance.
Locally, the MZR has been shown to exist over five decades in $M_{\ast}$ from
$10^{6}\lesssim M_{\ast}/\rm{M_{\odot}}\lesssim 10^{11}$ (e.g., Tremonti et
al., 2004; Lee et al., 2006; Kewley & Ellison, 2008; Andrews & Martini, 2013;
Maiolino & Mannucci, 2019; Curti et al., 2020; Sanders et al., 2021). The MZR
has also been shown in numerous studies to exist at high-$z$ out to $z>3$,
though with an evolution such that galaxies at fixed $M_{\ast}$ have lower
metallicities at higher redshifts (e.g., Erb et al., 2006; Maiolino et al.,
2008; Henry et al., 2013a, b; Zahid et al., 2013, 2014a, 2014b; Steidel et
al., 2014; Maiolino & Mannucci, 2019; Sanders et al., 2015, 2020, 2021; Strom
et al., 2022). Additionally, at higher redshifts, the turnover mass is found
to be larger than seen locally (Zahid et al., 2013, 2014a, 2014b), and at
$z>2$, it is unknown whether the turnover mass exists at all. At $z>2$ and
$M_{\ast}\gtrsim 10^{9}\ \rm{M_{\odot}}$, the MZR has been described by a
single power-law (e.g., Sanders et al., 2021).

In constraining the shape, scatter, and evolution in the MZR, insight is
gained into the physics of how star formation processes and baryon flows are
connected and how galaxy growth is structured and regulated (e.g., Finlator &
Davé, 2008; Davé et al., 2012; Ma et al., 2016; Torrey et al., 2014, 2019).
For instance, the slope of the low-mass end of the MZR can relate galactic
metal retention to how efficient outflows (which remove gas and metals from
the ISM/galaxy) and stellar feedback are in regulating star formation and
stellar mass growth (e.g., Torrey et al., 2014). Correlated scatter in the MZR
can inform of secondary dependencies of the metallicity to properties such as
gas-mass and SFR, giving further insight into current conditions of a galaxy
as well as elucidating more fundamental relationships between mass,
metallicity, and other properties (e.g., Ma et al., 2016; Torrey et al.,
2019).

Through the empirical study of the scatter in the MZR, Mannucci et al. (2010)
and Lara-López et al. (2010) found that metallicities of galaxies in the Sloan
Digital Sky Survey (SDSS) do in fact have a secondary dependence on SFR, a
dependence reaffirmed in more recent work, albeit to varying degrees of the
strength of that dependence (e.g., Yates et al., 2012; Andrews & Martini,
2013; Maiolino & Mannucci, 2019; Curti et al., 2020; Sanders et al., 2021).
This $M_{\ast}-\rm{SFR}-\rm{O/H}$ relation is referred to as the fundamental
metallicity relation (FMR) and displays a reduced scatter in metallicity of
$\sim 0.05$ dex (compared to $\sim 0.1$ dex in the MZR; Tremonti et al.,
2004). In effect, the FMR is a 3D surface that posits that metallicity is
anti-correlated with SFR such that, at fixed $M_{\ast}$, galaxies with above-
average (below-average) SFRs will have below-average (above-average) O/H. In
addition to a reduction in the intrinsic scatter of metallicity, Mannucci et
al. (2010) also suggested that the FMR is redshift-invariant out to $z\sim
2.5$. If true, this naturally explains the evolution in the normalization of
the MZR to be observations at various redshifts of different regions of the
locally-defined FMR; this is physically motivated by the anti-correlation of
O/H with SFR and the observed increase of SFR with redshift at fixed
$M_{\ast}$ (e.g., Speagle et al., 2014; Whitaker et al., 2014; Sanders et al.,
2021). As galaxy samples at high-$z$ have increased in size, evidence has
grown that the FMR is indeed redshift-invariant, at least to within $\sim 0.1$
dex, out to $z\sim 2.5$ and even possibly $z\sim 3.3$ (e.g., Henry et al.,
2013a, b; Cresci et al., 2019; Maiolino & Mannucci, 2019; Sanders et al.,
2018, 2020, 2021). However, this evidence is largely based on metallicities
indirectly-calculated via prescriptions calibrated in the local Universe, and
it is unknown how accurate these methods are at high-$z$.

In order to properly assess the evolution of the MZR and FMR, metallicities
must be estimated accurately at low- and high-$z$. This requires an accurate
understanding of the nebular physical conditions of star-forming galaxies at
different redshifts. Fortunately, a procedure that addresses both of these
requirements exists and is applicable at various redshifts, this procedure
being the “direct” method of oxygen abundance determination. This method
relies on first estimating the electron temperature ($T_{e}$) and electron
density ($n_{e}$) of the ionized nebular gas as these properties are
responsible for the strength of the collisionally-excited oxygen emission
lines needed for this procedure ([O II] $\lambda\lambda$3726, 3729 and [O III]
$\lambda\lambda$4959, 5007). These properties are then considered together
with flux ratios of the collisionally-excited lines to hydrogen Balmer
recombination lines in order to estimate the total oxygen abundance (e.g.,
Izotov et al., 2006; Osterbrock & Ferland, 2006). Unfortunately, this direct
method relies on weak auroral emission lines to calculate $T_{e}$, which is
determined from the flux ratio of strong emission lines (e.g., [O III]
$\lambda\lambda$4959, 5007) to auroral emission lines (e.g., O III]
$\lambda\lambda$1661, 1666 or [O III] $\lambda$4363) of the same ionic
species. While use of [O III] $\lambda$4363 is common in this methodology as
it lies in the rest-optical with [O III] $\lambda\lambda$4959, 5007, this line
is $\sim 30-100\times$ fainter than [O III] $\lambda$5007 (e.g., Jones et al.,
2015, Figure 1), typically decreasing in strength with increasing galactic
metallicity. As such, large, representative samples of [O III]
$\lambda$4363-emitters (and thus direct metallicities) have been difficult to
acquire with current facilities and instrumentation, especially at high O/H
(and $M_{\ast}$ by the MZR) and at $z>1$ where only a handful of [O III]
$\lambda$4363 detections exist, mostly thanks to gravitational lensing
(Brammer et al., 2012a; Christensen et al., 2012; Stark et al., 2013; James et
al., 2014; Patrício et al., 2018; Gburek et al., 2019; Sanders et al., 2016a,
2020). Moreover, due to the faintness of the auroral lines, both in the UV and
optical, the currently-detected auroral-line-emitters at $z>1$ are clearly
biased and are more representative of extreme emission line galaxies (EELGs)
than of “typical” star-forming galaxies seen at these redshifts (Sanders et
al., 2020). These high-$z$ auroral-line-emitters tend to fall well above the
mean $M_{\ast}-SFR$ relation defined by typical galaxies at a given redshift
and have flux ratios indicative of higher ionization parameters and lower
metallicities than average.

To overcome the current limitations of the direct-metallicity method, and
therein study more representative samples of galaxies across a wider dynamic
range of metallicities and redshifts, indirect “strong-line” methods of
determining oxygen abundance were developed (Jensen et al., 1976; Alloin et
al., 1979; Pagel et al., 1979). These methods allow for metallicity estimation
when [O III] $\lambda$4363 cannot be detected. Instead, strong-line methods
rely on locally-calibrated empirically (e.g., Pettini & Pagel, 2004; Jones et
al., 2015; Bian et al., 2018; Curti et al., 2020) or theoretically-determined
(e.g., McGaugh, 1991; Kewley & Dopita, 2002; Dopita et al., 2013) relations
between metallicity and flux ratios of strong, rest-optical, nebular emission
lines. However, while these strong-line methods have proven very useful in
understanding the enrichment of local galaxies, they have several drawbacks of
their own. For example, depending on the strong-line index and calibration
used, metallicity estimates can vary by up to 0.7 dex (Kewley & Ellison,
2008). In part, this is due to how the strong-line methods are calibrated.
Calibrations based on photoionization models tend to produce higher
metallicity estimates than empirical, $T_{e}$-based calibrations (Curti et
al., 2020, Figure 3). Empirical calibrations can also suffer from sample
selection effects whereby individually-detected [O III] $\lambda$4363-emitters
yield metallicities of more extreme star-forming regions whereas metallicities
from galaxy samples stacked in order to detect [O III] $\lambda$4363 may be
more representative of “typical” galaxies (Curti et al., 2017; Sanders et al.,
2020) that fall on the $M_{\ast}-SFR$ relation.

When considering high-$z$ galaxies, it is unknown if these locally-calibrated
strong-line relations, reflective of H II-region conditions in the local
Universe, are applicable for estimating metallicity. Excitation diagrams have
shown that star-forming region conditions likely evolve with redshift; this is
most notably seen in the [O III] $\lambda$5007/H$\beta$ vs. [N II]
$\lambda$6583/H$\alpha$ Baldwin$\,-\,$Phillips$\,-\,$Terlevich (N2-BPT;
Baldwin et al., 1981) diagram, where the locus of star-forming, high-$z$
galaxies is offset from the locus of local, star-forming, SDSS galaxies (e.g.,
Steidel et al., 2014; Shapley et al., 2015; Strom et al., 2017, 2018; Kashino
et al., 2017, 2019; Runco et al., 2022). While it is a current matter of
debate as to what is driving this evolution in the locus and thus the H II
region physical conditions (see Kewley et al. (2013) for an analysis of
several possibilities such as the ionization parameter, electron density,
hardness of the ionizing spectrum, and N/O abundance ratio), it is clear that
caution must be taken when applying strong-line metallicity methods at high
redshift. Calibrations are needed that are derived from objects with analogous
physical conditions to typical, star-forming, high-$z$ galaxies.

In this paper, we analyze a $\langle z\rangle=2.3$ composite spectrum of 16
gravitationally-lensed, typical, star-forming dwarf galaxies selected
independent of emission-line strength. In particular, we study the direct-
method metallicity from this composite, derived from a detection of the
$T_{e}$-sensitive [O III] $\lambda$4363 auroral-line. The paper is organized
as follows: In Section 2, we discuss our observations, data reduction, and
sample selection. In Section 3, we discuss our spectral-fitting and stacking
methodologies, introduce our composite spectrum, and calculate physical
properties of our stacking sample and composite. In Section 4, we present our
analysis and discussion in regard to how representative our sample is of
typical, $z\sim 2.3$, star-forming dwarf galaxies, the applicability of
locally-calibrated strong-line metallicity diagnostics at high-$z$, the slope
and normalization of the $z\sim 2.3$ MZR, and the redshift evolution of the
FMR. In Section 5, we summarize our results. Finally, in Appendix A, we
briefly describe the reasoning and methods behind our refitting of the Bian et
al. (2018) strong-line metallicity relations. Throughout this paper,
uncertainties reflect our bootstrapped error spectrum for the composite unless
stated otherwise. We assume a $\Lambda$CDM cosmology with $H_{0}$ = 70 km $\rm
s^{-1}$ $\rm Mpc^{-1}$, $\Omega_{\Lambda}$ = 0.7, and $\Omega_{m}$ = 0.3.

##  2 Observations, Data Reduction, and Sample Selection

The focus of this paper is the careful analysis of a stack$\,-\,$from which
the oxygen abundance is directly measured$\,-\,$of 16 gravitationally-lensed,
star-forming, dwarf galaxies at the peak of cosmic star formation. These
galaxies at $1.7<z<2.6$ have stellar masses of $\log(M_{\ast}/M_{\odot})<9.0$
and probe typical dwarf galaxies in this epoch, complimenting the recent large
statistical studies of more massive galaxies at these redshifts, such as the
Keck Baryonic Structure Survey (KBSS-MOSFIRE; Steidel et al., 2014) and the
MOSFIRE Deep Evolution Field survey (MOSDEF; Kriek et al., 2015). In this
section, we detail the photometric and spectroscopic observations and data
reduction of these galaxies and the larger parent surveys from which the
galaxies are drawn. We also discuss the selection strategy of these 16 objects
chosen for stacking.

###  2.1 Photometric Data and Reduction

The galaxy stacking sample is drawn from a spectroscopic follow-up survey of
the photometric Hubble Space Telescope (HST) survey of Alavi et al. (2014,
2016), which was conducted to study faint, low-mass, star-forming galaxies
gravitationally-lensed by the foreground galaxy clusters Abell 1689, MACS
J0717.5+3745, and MACS J1149.5+2223, among others (hereafter A1689, MACS
J0717, and MACS J1149, respectively). This HST survey compliments the Hubble
Frontier Fields (HFF; Lotz et al., 2017) survey of lensing clusters by both
adding deep near-ultraviolet (UV) images of the HFF clusters (of which MACS
J0717 and MACS J1149 are members) to the deep HFF optical and near-infrared
(IR) datasets as well as by adding or including deep near-UV to near-IR
photometry of another lensing cluster, A1689.

For galaxies lensed by A1689, near-UV images were taken over two programs in
the F225W, F275W, and F336W bandpasses with the Wide Field Camera 3
(WFC3)/UVIS channel on the HST. As part of Program ID 12201 (PI: B. Siana),
F275W was observed for 30 orbits, and F336W was observed for 4 orbits. As part
of Program ID 12931 (PI: B. Siana), F336W was observed for an additional 14
orbits (18 orbits total), and F225W was observed for 10 orbits. In the
optical, we used existing HST photometry, taken with the Advanced Camera for
Surveys (ACS)/WFC channel, in the F475W, F625W, F775W, and F850LP bandpasses
(PID: 9289; PI: H. Ford) as well as in the F814W bandpass (PID: 11710; PI: J.
Blakeslee). A summary of the number of orbits for each near-UV and optical
filter, as well as the $5\sigma$ depths for a $0\farcs 2$ radius aperture, can
be found in Alavi et al. (2016, Table 1). In the near-IR, existing images
taken over 1-2 orbits with the F125W and F160W filters and the HST WFC3/IR
channel (PID: 11802; PI: H. Ford) were used. We note that the near-IR
footprint for A1689 is smaller than the near-UV and optical footprints,
covering 10 of the 13 stacking sample galaxies (see sample selection in
Section 2.3) lensed by A1689.

Galaxies behind the lensing clusters MACS J0717 and MACS J1149 were observed
with the WFC3/UVIS channel for 8 orbits in both the F275W and F336W bandpasses
as part of the Alavi et al. (2016) HST survey under Program ID 13389 (PI: B.
Siana). In the optical and near-IR, these clusters were observed with HST
Director’s Discretionary time as part of the Hubble Frontier Fields survey
(Lotz et al., 2017). As with all clusters in this survey (6 total), MACS J0717
and MACS J1149 were observed for 70 orbits with each ACS/WFC and WFC3/IR (140
orbits total). These HFF clusters are observed in the F435W, F606W, and F814W
filters with ACS/WFC and the F105W, F125W, F140W, and F160W filters with
WFC3/IR (PID: 13498 for MACS J0717; PID: 13504 for MACS J1149; PI: J. Lotz).
As for the optical and near-UV photometry of A1689, the depths and orbits
(both from the HFF survey and other projects) for each filter are listed for
MACS J0717 in Alavi et al. (2016, Table 1). This information can be found for
MACS J1149 via the Mikulski Archive for Space Telescopes (MAST) website for
the HFF survey.222https://archive.stsci.edu/prepds/frontier/macs1149.html

The data reduction, calibration, and photometric measurements for MACS J0717
and MACS J1149 are detailed in Alavi et al. (2016), as is the UV data
reduction and calibration for A1689. The reduction and calibration of the
optical data from A1689, as well as the photometric measurements for this
cluster, are discussed in Alavi et al. (2014). The near-IR photometry of A1689
was reduced in the same way as the UV and optical data with the exception that
a larger pixel scale of $0\farcs 08$ was used in the final drizzled images. As
described in Alavi et al. (2014, 2016), our main photometric catalog for A1689
is built on the UV and optical images with a pixel scale of $0\farcs 04$. For
the areas of A1689 with near-IR coverage, the multi-band photometry (from UV
to near-IR) was remeasured on images with larger pixel scales and that are
PSF-matched to the F160W data. The estimations of photometric redshifts, which
were used to select the spectroscopic follow-up survey sample detailed in
Section 2.2, are described in Alavi et al. (2016).

####  2.1.1 Lens Models

When working with objects gravitationally-lensed by foreground galaxy
clusters, accurate lens models are imperative for correcting observed
photometry and spectroscopy for the lensing magnification. This correction is
necessary for the determination of an object’s intrinsic properties (e.g.,
stellar mass, SFR, etc.). Alavi et al. (2016) detail the lens models
considered and used for the HFF clusters and A1689, all of which, while
constructed with different assumptions and methodologies, are constrained by
the location and redshift of known multiply-imaged systems. As stated in Alavi
et al. (2016), for the HFF clusters we use the lens models derived by the
Clusters As TelescopeS (CATS)
collaboration,333https://archive.stsci.edu/prepds/frontier/lensmodels/
specifically the models of Limousin et al. (2016) and Jauzac et al. (2016) for
MACS J0717 and MACS J1149, respectively. For A1689, we use the lens model of
Limousin et al. (2007). These parametric models are all derived via mass
reconstruction done with the
LENSTOOL444https://projets.lam.fr/projects/lenstool/wiki software (Jullo et
al., 2007).

###  2.2 Spectroscopic Data and Reduction

As a follow-up to the photometric HST survey of Alavi et al. (2014, 2016), a
spectroscopic survey was conducted between 2014 January and 2017 March to
obtain near-IR (rest-optical) spectroscopy of select galaxies with the Multi-
Object Spectrometer For InfraRed Exploration (MOSFIRE; McLean et al., 2010,
2012) on the 10 m Keck I telescope. Galaxies for this survey were selected to
have high magnifications, observed optical magnitudes (F606W or F625W) less
than 26.0 (AB), and photometric redshifts in three redshift ranges,
$1.37\leqslant z\leqslant 1.70$, $2.09\leqslant z\leqslant 2.61$, and
$2.95\leqslant z\leqslant 3.80$, so that the strong, rest-optical, nebular
emission lines of the galaxies lie in the near-IR atmospheric transmission
windows. Early selection of galaxies in MACS J0717 and MACS J1149 used
photometric redshifts from the CLASH survey (Postman et al., 2012). In all,
151 sources were observed across 9 masks. For galaxies that fall into the two
lowest redshift ranges, the strong, nebular emission lines targeted are [O II]
$\lambda\lambda$3726, 3729, H$\beta$, [O III] $\lambda\lambda$4959, 5007,
H$\alpha$, and [N II] $\lambda\lambda$6548, 6583. To this end, observations of
galaxies in the lowest redshift range were conducted using the Y-, J-, and
H-band filters, whereas the J-, H-, and K-band filters were used for the two
highest redshift ranges. We note that while we targeted the strong, nebular
emission lines in the highest redshift range as well, H$\alpha$ and the [N II]
doublet were not observed as they fall outside of the K-band’s wavelength
coverage.

Observations used an ABBA dither pattern with a $2\farcs 5$ dither spacing.
The individual exposure time for J-band and H-band data was 120 s and was 180
s for Y-band and K-band data. In total, across the 9 masks, the J-band was
observed between 48 m and 112 m, the H-band between 56 m and 112 m, and the
K-band between 60 m and 120 m, for average total exposure times of 81 m, 85 m,
and 82 m, respectively. Data in the Y-band were taken for one mask in A1689
for a total of 96 m. In each mask, we used $0\farcs 7$-wide slits, yielding
spectral resolutions of $R=3388,\ 3318,\ 3660,\ \rm{and}\ 3610$ for the Y-,
J-, H-, and K-bands,
respectively.555https://www2.keck.hawaii.edu/inst/mosfire/grating.html Our
typical FWHM seeing for a given mask/filter combination was $0\farcs 71$.

The spectroscopic data obtained with MOSFIRE were reduced with the MOSFIRE
Data Reduction Pipeline666https://keck-
datareductionpipelines.github.io/MosfireDRP/ (DRP). This DRP returns a 2D
flat-fielded, wavelength-calibrated, background-subtracted, and rectified
spectrum for each slit in a given mask. This output spectrum is a composite of
the multiple spectra taken at each nod position. For Y-, J-, and H-band
spectra, wavelength calibration is performed using the night-sky lines,
whereas a combination of night-sky lines and a neon arc lamp is used for
K-band spectra owing to the faintness of the sky lines and the dominance of
thermal noise at the red end of the band. Once the 2D spectra were produced,
the 1D spectra were extracted using the custom IDL software
BMEP777https://github.com/billfreeman44/bmep from Freeman et al. (2019). Each
spectrum is flux-calibrated with two stars. A standard star of spectral type
B9 V to A2 V is first used to apply a wavelength-dependent calibration. It is
ensured that this standard star was observed at an air mass similar to that of
the mask under consideration. Following this step, an absolute flux
calibration is conducted using a star that was included in the corresponding
mask.

###  2.3 Sample Selection for Dwarf Galaxy Stack

The galaxies that comprise the stack mentioned in the opening of this section
are drawn from the photometric and spectroscopic surveys detailed above. These
galaxies are required to have a robust spectroscopic redshift and
spectroscopic coverage of the strong, rest-optical, nebular emission lines: [O
II] $\lambda\lambda$3726, 3729, H$\beta$, [O III] $\lambda$4959, H$\alpha$,
and [N II] $\lambda\lambda$6548, 6583. Additionally, these galaxies must have
spectroscopic coverage of H$\gamma$ and the faint [O III] $\lambda$4363
auroral emission line. The auroral line is essential for determining gas-phase
metallicity directly as it is a component of the emission-line ratio used to
estimate electron temperature ($T_{e}$; see Section 3.5 for more details).
These redshift and coverage requirements yield a sample of 18 galaxies and 24
total spectra when accounting for multiply-imaged systems, of which we have
four in our sample. A final cut is made on stellar mass (see Section 3.2 on
mass estimation) to ensure that our sample lies in the dwarf galaxy regime
($\log(M_{\ast}/M_{\odot})<9.0$). With this cut, two galaxies are removed from
our sample, yielding a final count of 16 galaxies (22 total spectra) ranging
in redshift from $z=1.70$ to $z=2.59$ ($z_{\rm mean}=2.30$).

We note here that H$\gamma$ coverage is included as a requirement so as to
provide another Balmer decrement with which to estimate the dust extinction
from the stack. Due to the close proximity of H$\gamma$ (4340 $\rm{\AA}$) to
[O III] $\lambda$4363, this inclusion does not affect our sample size. We also
note here that we do not require spectroscopic coverage of the [O III]
$\lambda$5007 line of the [O III] $\lambda\lambda$4959, 5007 doublet so as to
maximize our galaxy count by including those sources for which [O III]
$\lambda$5007 falls just redward of a given filter. Instead, when necessary,
we make use of the $T_{e}$-insensitive intrinsic intensity ratio of the
doublet: [O III] $\lambda$5007/[O III] $\lambda$4959 = 2.98 (Storey & Zeippen,
2000). Lastly, while we do not select galaxies based on the strength of any
given emission line, we do note that each spectrum has a signal-to-noise ratio
(S/N) for [O III] $\lambda$5007 of $\rm{S/N}>5$, ensuring accurate
normalization of each spectrum (by [O III] $\lambda$5007) during the stacking
process (see Section 3.3). A summary of our sample, and some of the galaxies’
physical properties, are listed in Table 1.

Table 1: Summary and Properties of $\langle z\rangle=2.3$ Dwarf Galaxy Stacking Sample Galaxy | Spec. ID | $z\tablenotemark{a}$ | R.A.bbRight Ascension: hh:mm:ss.sss; Declination: dd:mm:ss.sss; Equinox: J2000 | Dec.bbRight Ascension: hh:mm:ss.sss; Declination: dd:mm:ss.sss; Equinox: J2000 | $\log(\frac{M_{\ast}}{\rm{M_{\odot}}})_{\rm{fid}}\tablenotemark{c}$ | $\log(\frac{M_{\ast}}{\rm{M_{\odot}}})\tablenotemark{d}$ | $\frac{\rm{SFR}}{\rm{M_{\odot}\ yr^{-1}}}\tablenotemark{e}$ | $\log(\frac{L_{\rm{H\alpha}}}{\rm{erg\ s^{-1}}})\tablenotemark{f}$
---|---|---|---|---|---|---|---|---
1 | A1689-1037 | 1.70089 | 13:11:35.197 | -01:20:25.040 | $7.71^{+0.20}_{-0.37}$ | $8.12^{+0.22}_{-0.23}$ | $0.151\pm 0.004$ | $40.511^{+0.013}_{-0.013}$
2 | A1689-1197 | 1.70261 | 13:11:29.689 | -01:20:08.769 | $8.45^{+0.05}_{-0.05}$ | $8.87^{+0.13}_{-0.10}$ | $-\tablenotemark{h}$ | $-\tablenotemark{h}$
| A1689-370 | 1.70257 | 13:11:32.406 | -01:21:16.027 | $8.31^{+0.05}_{-0.05}$ | $8.71^{+0.12}_{-0.11}$ | $1.779\pm 0.020$ | $41.583^{+0.005}_{-0.005}$
| CompositeggThe redshift, stellar masses, and H$\alpha$ luminosity are weighted-averages. The SFR is calculated from this luminosity. | 1.70259 | $-$ | $-$ | $8.36^{+0.03}_{-0.04}$ | $8.79^{+0.08}_{-0.08}$ | $1.779\pm 0.020$ | $41.583^{+0.005}_{-0.005}$
3 | A1689-280 | 1.70316 | 13:11:31.886 | -01:21:26.014 | $7.91^{+0.04}_{-0.05}$ | $8.70^{+0.23}_{-0.21}$ | $0.692\pm 0.011$ | $41.173^{+0.007}_{-0.007}$
4 | A1689-257 | 1.70355 | 13:11:26.426 | -01:21:31.277 | $7.81^{+0.08}_{-0.10}$ | $8.47^{+0.18}_{-0.14}$ | $1.809\pm 0.027$ | $41.590^{+0.006}_{-0.007}$
5 | A1689-1751 | 2.38159 | 13:11:31.333 | -01:19:18.559 | $8.82^{+0.07}_{-0.09}$ | $9.23^{+0.17}_{-0.12}$ | $3.071\pm 0.145$ | $41.820^{+0.020}_{-0.021}$
6 | A1689-232 | 2.38709 | 13:11:32.794 | -01:21:27.893 | $7.06^{+1.10}_{-7.06}$ | $8.70^{+0.20}_{-0.19}$ | $2.368\pm 0.185$ | $41.707^{+0.033}_{-0.035}$
7 | M0717-3958 | 2.39329 | 07:17:27.442 | +37:45:25.475 | $8.69^{+0.11}_{-0.16}$ | $9.42^{+0.09}_{-0.12}$ | $2.072\pm 0.139$ | $41.649^{+0.028}_{-0.030}$
| M0717-4517 | 2.39330 | 07:17:27.050 | +37:45:09.695 | $8.99^{+0.07}_{-0.09}$ | $9.46^{+0.09}_{-0.09}$ | $2.350\pm 0.128$ | $41.704^{+0.023}_{-0.024}$
| CompositeggThe redshift, stellar masses, and H$\alpha$ luminosity are weighted-averages. The SFR is calculated from this luminosity. | 2.39329 | $-$ | $-$ | $8.84^{+0.07}_{-0.08}$ | $9.44^{+0.07}_{-0.07}$ | $2.223\pm 0.094$ | $41.680^{+0.018}_{-0.019}$
8 | A1689-1059 | 2.41141 | 13:11:25.228 | -01:20:19.309 | $8.40^{+0.03}_{-0.03}$ | $9.03^{+0.17}_{-0.21}$ | $5.041\pm 0.464$ | $42.035^{+0.038}_{-0.042}$
9 | A1689-1467 | 2.51903 | 13:11:26.118 | -01:19:42.837 | $8.21^{+0.10}_{-0.13}$ | $8.76^{+0.18}_{-0.17}$ | $2.287\pm 0.441$ | $41.692^{+0.077}_{-0.093}$
10 | A1689-1216 | 2.54082 | 13:11:31.981 | -01:20:07.173 | $8.70^{+0.06}_{-0.07}$ | $9.07^{+0.15}_{-0.15}$ | $0.822\pm 0.144$ | $41.248^{+0.070}_{-0.084}$
| A1689-1292 | 2.54064 | 13:11:26.528 | -01:19:55.146 | $8.90^{+0.06}_{-0.07}$ | $9.07^{+0.11}_{-0.13}$ | $0.667\pm 0.097$ | $41.157^{+0.059}_{-0.068}$
| A1689-537 | 2.54046 | 13:11:29.795 | -01:21:05.969 | $8.60^{+0.08}_{-0.10}$ | $9.05^{+0.17}_{-0.19}$ | $1.725\pm 0.397$ | $41.570^{+0.090}_{-0.114}$
| CompositeggThe redshift, stellar masses, and H$\alpha$ luminosity are weighted-averages. The SFR is calculated from this luminosity. | 2.54065 | $-$ | $-$ | $8.70^{+0.04}_{-0.04}$ | $9.07^{+0.08}_{-0.08}$ | $0.755\pm 0.079$ | $41.211^{+0.043}_{-0.048}$
11 | A1689-470 | 2.54112 | 13:11:26.213 | -01:21:09.695 | $7.28^{+0.31}_{-7.28}$ | $8.49^{+0.22}_{-0.21}$ | $1.235\pm 0.160$ | $41.425^{+0.053}_{-0.060}$
12 | A1689-1451 | 2.54201 | 13:11:28.682 | -01:19:42.849 | $8.22^{+0.08}_{-0.09}$ | $9.10^{+0.17}_{-0.16}$ | $0.756\pm 0.128$ | $41.211^{+0.068}_{-0.080}$
13 | A1689-722 | 2.54247 | 13:11:33.915 | -01:20:52.526 | $8.78^{+0.20}_{-0.36}$ | $8.81^{+0.22}_{-0.24}$ | $3.771\pm 0.291$ | $41.910^{+0.032}_{-0.035}$
14 | M0717-1531 | 2.55185 | 07:17:32.547 | +37:45:02.348 | $9.08^{+0.11}_{-0.14}$ | $9.94^{+0.25}_{-0.06}$ | $5.558\pm 0.441$ | $42.078^{+0.033}_{-0.036}$
| M0717-3187 | 2.55159 | 07:17:35.089 | +37:45:48.120 | $8.90^{+0.10}_{-0.14}$ | $9.67^{+0.17}_{-0.15}$ | $5.965\pm 0.268$ | $42.109^{+0.019}_{-0.020}$
| M0717-5970 | 2.55167 | 07:17:30.613 | +37:44:22.798 | $8.82^{+0.17}_{-0.29}$ | $9.66^{+0.08}_{-0.10}$ | $3.969\pm 0.283$ | $41.932^{+0.030}_{-0.032}$
| CompositeggThe redshift, stellar masses, and H$\alpha$ luminosity are weighted-averages. The SFR is calculated from this luminosity. | 2.55168 | $-$ | $-$ | $8.93^{+0.07}_{-0.09}$ | $9.72^{+0.07}_{-0.07}$ | $5.109\pm 0.178$ | $42.041^{+0.015}_{-0.015}$
15 | A1689-217 | 2.59181 | 13:11:27.623 | -01:21:35.622 | $8.23^{+0.04}_{-0.04}$ | $9.22^{+0.17}_{-0.18}$ | $9.194\pm 0.313$ | $42.297^{+0.015}_{-0.015}$
16 | M1149-2185 | 2.59366 | 11:49:40.162 | +22:25:07.571 | $8.83^{+0.05}_{-0.06}$ | $9.33^{+0.04}_{-0.04}$ | $33.602\pm 3.037$ | $42.859^{+0.038}_{-0.041}$

##  3 Measurements and Stacking Methodology

In this section, we detail our methodologies for fitting the spectroscopy and
photometry of the dwarf galaxies in our stacking sample. We also discuss how
various physical properties are estimated either for the individual galaxies
or for the “sample-average” dwarf galaxy, represented by a composite spectrum
of these dwarfs. We begin by discussing the measurements made for individual
galaxies and then proceed to the construction and analysis of the composite
spectrum.

###  3.1 Fitting the Individual Emission-Line Spectra

Each emission-line spectrum in our stacking sample (22 total), corresponding
to either the single image of a galaxy or one of a multiply-imaged galaxy, is
fit using the Markov Chain Monte Carlo Ensemble sampler
emcee888https://emcee.readthedocs.io/en/v2.2.1/ (Foreman-Mackey et al., 2013).
Lines are fit to each spectrum’s continuum, and single-Gaussian profiles are
fit to the emission lines. To minimize the impact on the spectral-fitting from
pixels contaminated by sky lines, we removed, prior to fitting, any pixels
with a corresponding error spectrum value $>3\times$ the median error value
over the range of the fit.

When fitting, each spectroscopic band (Y, J, H, K) was considered separately.
For each spectrum, the slope and intercept of the continuum were free
parameters. In the band containing H$\alpha$ and the [N II] doublet (the H\-
or the K-band), the free parameters also included the redshift of the
spectrum, the width of the emission lines (each line having the same width),
and the amplitudes of the H$\alpha$ and [N II] $\lambda$6583 lines, with the
amplitude of [N II] $\lambda$6583 constrained such that [N II]
$\lambda$6583/[N II] $\lambda$6548 = 2.95 (Acker et al., 1989).

In the band containing H$\gamma$, [O III] $\lambda$4363, H$\beta$, and [O III]
$\lambda$4959 (the J\- or H-band), two fits were conducted due to the large
wavelength separation between [O III] $\lambda$4363 and H$\beta$. The portion
of the spectrum containing H$\beta$ and [O III] $\lambda$4959 (and [O III]
$\lambda$5007 if covered) was fit first, having the free parameters of line-
width for the filter, redshift, and emission-line amplitudes. If [O III]
$\lambda$5007 is within the spectrum’s wavelength coverage, its amplitude was
fit with the constraint that [O III] $\lambda$5007/[O III] $\lambda$4959 =
2.98 (Storey & Zeippen, 2000). Otherwise, the amplitude of [O III]
$\lambda$4959 was fit, and the line’s flux was multiplied by the
aforementioned intensity ratio in order to estimate the [O III] $\lambda$5007
flux. With a best-fit width and redshift in-hand from the first fit to the
filter, the fainter H$\gamma$ and [O III] $\lambda$4363 lines were then fit
with these two parameters fixed.

Finally, in the band containing the [O II] $\lambda\lambda$3726, 3729 doublet
(the Y\- or J-band), the redshift and line-width were fixed to the values fit
to the complete spectrum’s highest S/N line in order to avoid complications
resulting from the doublet lines’ small wavelength separation. In addition to
the continuum parameters, only the [O II] lines’ amplitudes were free
parameters in these fits.

Ultimately, the final redshift given to the full spectrum is the weighted
average of the redshifts fit to the J (H)- and H (K)-bands.

####  3.1.1 Slit-Loss Correction

When measuring the emission-line fluxes from spectra observed through slit
masks, care must be taken to account for loss of flux outside of the slits in
order to recover the full integrated flux values. To this end, our emission-
line fluxes were slit-loss-corrected on a galaxy-by-galaxy basis using the
methodology of Emami et al. (2020).

###  3.2 SED-Fitting and Stellar Mass Estimation

To determine the stellar masses of the galaxies in our stacking sample, we fit
spectral energy distributions (SEDs) to our HST near-UV to near-IR photometry
(we note that three galaxies lensed by A1689 lack near-IR photometry; see
Section 2.1). At high redshift, observations suggest that high equivalent
width emission lines are fairly common, particularly in lower-mass galaxies
like those in our sample (Reddy et al., 2018). Therefore, prior to SED-
fitting, we subtracted off any contribution to the photometry from the slit-
loss-corrected, nebular emission lines. We also added an additional $3\%$ flux
error, in quadrature, to all bands in order to account for systematic errors
in the photometry (Alavi et al., 2016). To this emission-line-corrected
photometry, we then fit Bruzual & Charlot (2003) stellar population synthesis
models using the SED-fitting code
FAST999https://w.astro.berkeley.edu/~mariska/FAST.html (Kriek et al., 2009).
We assume constant star formation histories (SFH), a Chabrier (2003) initial
mass function (IMF), stellar metallicities of 0.2 $\rm{Z_{\odot}}$ or 0.4
$\rm{Z_{\odot}}$, and a Calzetti et al. (2000) dust attenuation curve. The
redshifts of the galaxies are fixed to their fit spectroscopic values.

Uncertainties on the properties estimated by FAST (e.g., stellar mass, SFR,
$A_{V}$, etc.) are derived using a Monte Carlo approach where the photometry
being fit is perturbed within its uncertainties and is then refit, this
process being repeated 300 times. From these 300 realizations of the SED,
$68\%$ confidence intervals are determined for each estimated property. In
Table 1, we list the best-fit stellar mass, and its uncertainty, of each
galaxy in our stacking sample. The stellar mass associated with our full-
sample composite, detailed below, is taken to be the median of these
individual masses, log($M_{\ast}$/$\rm{M_{\odot}}$)med = 8.29. We note that
the best-fit SEDs, stellar masses, and all other affected properties are de-
magnified based on the lensing models discussed in Section 2.1.1.

###  3.3 The Composite Spectrum

While the individual galaxies in our sample display several nebular emission
lines at high-S/N (e.g., [O III], H$\alpha$), the galaxies are still
inherently faint even with high magnification via gravitational-lensing. As a
result, many other useful, fainter lines are undetected or marginally-detected
in our individual spectra. Such lines can include H$\beta$ and H$\gamma$ for
estimating extinction from dust, the [O II] doublet for calculating electron
density, and especially [O III] $\lambda$4363 for estimating electron
temperature and metallicity directly. A composite spectrum, or stacked
spectrum, of all of our sample galaxies offers a solution to this problem by
including in our study both galaxies for which we have individual line
detections and galaxies for which we only have upper limits. This composite
gives the advantages of both increasing the S/N of faint spectral features and
displaying the average spectrum and properties of dwarf galaxies like those in
our sample. Additionally, we use composites of the individual spectra of
galaxies multiply-imaged by lensing in order to increase the effective
exposure times and S/N of those galaxies’ spectroscopy.

Our methodology for creating composite spectra is similar when stacking
multiple images of sources (4 multiply-imaged galaxies; see Table 1) or all of
the galaxies in our sample (16 total). We first create the composites for our
multiply-imaged galaxies as these composites represent their corresponding
galaxies in the full-sample stack. We begin by shifting the slit-loss-
corrected, observed spectra to the rest-frame and converting the flux
densities into luminosity densities. Each spectrum is then normalized by its
slit-loss-corrected, [O III] $\lambda$5007 emission-line luminosity. This
normalization serves two purposes. It de-magnifies each spectrum implicitly by
dividing the magnified spectrum by its magnified [O III] $\lambda$5007
luminosity. It also, in the case of stacking our full sample, prevents our
composite electron temperature (see Section 3.5) from being biased by the
brightest [O III] $\lambda$5007 source (Sanders et al., 2020). We note here
that, prior to normalizing, the spectra and [O III] $\lambda$5007 luminosities
are not corrected for dust extinction both due to the faintness of H$\beta$
and H$\gamma$ and the sky line contamination of these lines in several of our
individual sources. Following normalization, each spectroscopic band’s science
spectrum and corresponding error spectrum are resampled with the Python tool,
SpectRes101010https://spectres.readthedocs.io/en/latest/,111111https://github.com/ACCarnall/SpectRes
(Carnall, 2017), onto a common wavelength grid with a rest-frame dispersion
$-$ for the full-sample stack $-$ of 0.38 $\rm{\AA\ pix^{-1}}$ in the band (Y
or J) containing [O II], 0.47 $\rm{\AA\ pix^{-1}}$ in the band (J or H)
containing H$\gamma$ through [O III] $\lambda$4959, and 0.63 $\rm{\AA\
pix^{-1}}$ in the band (H or K) containing H$\alpha$ and [N II]. (Hereafter,
these bands will be referred to as the YJ-band, JH-band, and HK-band,
respectively.) These rest-frame dispersions are computed by shifting the
MOSFIRE J-, H-, and K-band observed-frame
dispersions${}^{\rm{\ref{footnote:grating}}}$ to the median redshift of the
stacking sample, $z_{\rm{med}}\approx 2.465$. Once resampled, the spectra from
the full galaxy sample are combined at each wavelength element by taking the
median value of all luminosity densities at that point. When stacking the
spectra of a multiply-imaged galaxy, the average at each pixel is taken
instead, with luminosity densities weighted by their associated uncertainty
values. Finally, the composite spectrum of the full stacking sample is
multiplied by the median [O III] $\lambda$5007 luminosity of the sample,
whereas the composite for a multiply-imaged galaxy is multiplied by the
lowest-luminosity [O III] $\lambda$5007 measurement (a proxy for the least
magnified measurement). We note that, in this work, only ratios of emission
lines are used from the composite of the full stacking sample since individual
luminosity measurements rely on a normalization dependent on our stacking
methodology.

Figure 1: The rest-optical composite spectrum of the 16 typical, star-forming,
$\langle z\rangle=2.30$ dwarf galaxies in our stacking sample. The stacked
spectrum is shown as the blue unfilled histogram, and its best-fit model is
displayed in red. Here we show two offset uncertainty spectra (see Section
3.3), one estimated via bootstrapping in order to account for cosmic variance
(the light blue-gray-filled histogram), and the other a statistical-only error
spectrum (the dark blue-gray-filled histogram). We note that the bootstrapped
uncertainty spectrum was used when fitting the model shown here. In the inset
box corresponding to the green-shaded spectral region, we show a zoom-in of
the weak [O III] $\lambda$4363 auroral line, measured at 2.5$\sigma$
(4.1$\sigma$) significance when considering the bootstrapped (statistical-
only) uncertainty spectrum. We note that, as a result of our stacking
procedure (see Section 3.3), the luminosity density values ($L_{\lambda}$) of
the composite have a constant but arbitrary magnification. While this makes
individual emission-line luminosities unreliable, it does not affect strong-
line ratios.

The uncertainty spectrum of the stack of all of our sample galaxies is derived
via a Monte Carlo approach with bootstrapping. We first create a bootstrapped
sample of number count $N_{\rm{boot}}=16$, the number of galaxies in our full
stack, by randomly drawing galaxies for the sample with replacement. For each
galaxy in this bootstrapped sample, its science spectrum (already shifted and
converted to the rest-frame and luminosity densities, respectively) is
perturbed according to its error spectrum. This perturbed spectrum is then
normalized by its corresponding [O III] $\lambda$5007 emission-line
luminosity, which has been perturbed based on its own uncertainty value. The
normalized, perturbed spectra are then resampled and stacked according to the
procedure detailed above. This process is repeated 500 times in order to
create an array of composite luminosity densities at each wavelength element.
Our uncertainty spectrum is then comprised of the standard deviations of the
values in each of these arrays. By constructing our error spectrum via
bootstrapping, our uncertainties represent both our measurement errors and
sample variance. The composite spectrum of our full stacking sample, as well
as its bootstrapped uncertainty spectrum (the light-gray-shaded region), are
shown in Figure 1. In Figure 1, we also show the statistical-only uncertainty
spectrum (the dark-gray-shaded region), which was created as described above,
but without bootstrap-resampling each iteration. We note that, unless
specified otherwise, the stated uncertainties on measurements derived from the
composite spectrum reflect the use of the bootstrapped uncertainty spectrum
and its consideration of sample variance.

####  3.3.1 Fitting the Composite Emission-Line Spectrum

Prior to their inclusion in the full-sample stack, the composite spectra of
the multiply-imaged galaxies are fit in the same manner as the spectra of
individual galaxy images (see Section 3.1). For the full-sample composite, we
adopt a slightly different fitting methodology. While spectral-fitting is
still facilitated with emcee, and the continuum is still fit with a line, the
emission lines are fit with a combination of two Gaussian profiles due to the
non-Gaussian shape of the high-S/N lines (Steidel et al., 2016). The resultant
profiles fit to the emission lines consist of a “broad” Gaussian component and
a “narrow” Gaussian component, both of which are centered on the rest-frame
wavelength of the emission line in question. In order to maintain consistent
line profiles for all emission lines in a given spectral band (YJ, JH, or HK),
the fitting of these profiles, and the spectrum, is done in two rounds.

In the first round, the two Gaussian components are fit to the brighter,
higher-S/N lines in the composite: H$\beta$ and [O III] $\lambda$4959 in the
JH-band and, separately, H$\alpha$ in the HK-band. For each spectral band, the
broad Gaussian has its $1\sigma$-width fixed at $100\ \rm{km\ s^{-1}}$ while
the narrow-Gaussian $1\sigma$-width is left as a free parameter.121212We
tested the validity of fixing the broad-Gaussian $1\sigma$-width by comparing
luminosities fit to H$\beta$, [O III] $\lambda$4959, and H$\alpha$ with either
a fixed broad-width ($\sigma_{\rm{b}}=100\ \rm{km\ s^{-1}}$) or the broad-
width left as a free parameter. When comparing the luminosities assuming a
fixed-width to the weighted-average luminosities of three runs where the
broad-width was left free, we found an average percentage difference of $\sim
0.4\%$ and $\sim 1.0\%$ when assuming the bootstrapped or statistical-only
uncertainty spectrum, respectively. The weighted-average luminosities fell
well within the uncertainties of the fixed-width luminosities. We therefore
find our choice of fixed broad-Gaussian $1\sigma$-width to be robust.
Additionally, the free parameters for each fit include the slope and intercept
of the continuum, the amplitude(s) of the narrow-Gaussian profile(s), and a
narrow-to-broad amplitude ratio for the spectral band.

In the second round, each spectral band is fit in its entirety. The broad
Gaussian keeps its $1\sigma$-width fixed to $100\ \rm{km\ s^{-1}}$ in all
bands. In the JH\- and HK-bands, the narrow $1\sigma$-widths and amplitude
ratios fit in the first round are held fixed and applied to all lines in the
corresponding bands: H$\gamma$, [O III] $\lambda$4363, H$\beta$, [O III]
$\lambda$4959 in the JH-band and H$\alpha$ and [N II] $\lambda\lambda$6548,
6583 in the HK-band. Free parameters in these bands are the narrow-Gaussian
amplitudes of each line and the linear continuum parameters. (The [N II]
$\lambda$6583 narrow-Gaussian amplitude is constrained in the same manner as
in the individual spectra; see Section 3.1.) In the YJ-band containing [O II]
$\lambda\lambda$3726, 3729, the amplitude ratio from the JH-band is adopted,
but the narrow $1\sigma$-width is left as a free parameter. Like with the
other bands, the narrow-Gaussian amplitudes and linear continuum parameters
are also fit.

The spectral model from these two rounds of fitting can be seen as the red
line in Figure 1. The total luminosity of each emission line, representing the
addition of the broad- and narrow-Gaussian component luminosities, is given in
Table 2 relative to the total H$\beta$ luminosity. Of particular interest for
this study is the $2.5\sigma$ ($4.1\sigma$) detection of [O III] $\lambda$4363
in our composite with the bootstrapped (statistical-only) error spectrum,
which will be used in Sections 3.5 and 3.6 to estimate the electron
temperature and gas-phase metallicity directly.

Table 2: Emission-Line Luminosities of the Composite Spectrum

Line | $\lambda_{\rm rest}\tablenotemark{a}$ | $L_{\rm{meas}}\tablenotemark{b}$ | $L_{\rm{corr}}\tablenotemark{b,c}$ | $A_{\lambda}\tablenotemark{d}$
---|---|---|---|---
$[$O II$]$ | 3726.032 | $0.61\pm 0.12$ | $0.70\pm 0.13$ | 0.67
$[$O II$]$ | 3728.815 | $0.95\pm 0.16$ | $1.09\pm 0.19$ | 0.67
H$\gamma$ | 4340.459 | $0.40\pm 0.08$ | $0.44\pm 0.09$ | 0.59
$[$O III$]$ | 4363.209 | $0.11\pm 0.05$ | $0.11\pm 0.05$ | 0.58
H$\beta$ | 4861.321 | $1.00\pm 0.16$ | $1.00\pm 0.16$ | 0.51
$[$O III$]$ | 4958.910 | $1.95\pm 0.30$ | $1.89\pm 0.29$ | 0.49
H$\alpha$ | 6562.794 | $3.14\pm 0.48$ | $2.70\pm 0.42$ | 0.36
$[$N II$]$ | 6583.448 | $0.13\pm 0.07$ | $0.11\pm 0.06$ | 0.35

ccfootnotetext: All luminosities are dust-corrected using the corresponding
extinction magnitude ($A_{\lambda}$) given in the last column. Balmer lines
are corrected for stellar absorption. See Section 3.4. The listed
uncertainties do not include systematic errors associated with the dust
correction, though these errors are propagated throughout all of our
calculations.

Note. — The luminosity of [O III] $\lambda$5007 can be calculated via the
intrinsic ratio [O III] $\lambda$5007/[O III] $\lambda$4959 = 2.98.

###  3.4 Dust Extinction and SFRs of the Sample

In order to estimate a galaxy’s intrinsic emission-line luminosities, from
which its galactic properties and interstellar medium (ISM) conditions are
derived, a wavelength-dependent correction to the observed luminosities must
be made to account for extinction from nebular dust. This correction is
typically quantified via Balmer recombination-line ratios of observed hydrogen
emission-line luminosities. Ideally, dust extinction would have been
compensated for on a galaxy-by-galaxy basis prior to stacking our sample.
Unfortunately, many of our individual galaxy spectra have Balmer lines that
are too faint or too impacted by sky lines for this approach to be used.
Instead, the nebular dust extinction “typical” of star-forming, dwarf galaxies
like those in our sample is estimated via our full-sample composite spectrum.

Prior to using Balmer decrements to estimate dust extinction, corrections to
Balmer emission-line luminosities must be made to account for Balmer
absorption in the atmospheres of (primarily A-type) stars. Using the Balmer
emission lines from our individual spectra and the model continua from our
SED-fitting, we estimated the sample-median stellar absorption corrections for
H$\gamma$, H$\beta$, and H$\alpha$ and applied them to our composite Balmer
emission-line luminosities, increasing the measured values by $\sim 5.0\%$,
$\sim 1.7\%$, and $\sim 0.4\%$, respectively.

To calculate the nebular dust extinction with our composite spectrum, we
assumed Case B intrinsic Balmer ratios of H$\alpha$/H$\beta$ = 2.79 and
H$\alpha$/H$\gamma$ = 5.90 for an electron temperature and electron density of
$T_{e}$([O III]) = 15,000 K and $n_{e}$ = 100 $\rm{cm^{-3}}$,
respectively131313We note that while an electron temperature is assumed when
selecting Balmer ratio values, the dependence of those ratios on $T_{e}$ is
weak for typical temperatures in H II regions. To confirm that our assumption
of $T_{e}$ = 15,000 K is valid, we dust-corrected the composite spectrum then
calculated $T_{e}$ (see Section 3.5) assuming intrinsic Balmer ratios
corresponding to $T_{e}$ = 10,000 K, 12,500 K, 15,000 K, and 20,000 K (Dopita
& Sutherland, 2003). With each variation, we consistently calculated from our
dust-corrected composite spectrum a $T_{e}\sim\ $15,000 K. (Dopita &
Sutherland, 2003). Further assuming the extinction curve of Cardelli et al.
(1989) with $R_{V}$ = 3.1, we find a “typical” color excess of
$E(\bv)_{\rm{gas}}=A_{V}/R_{V}=0.14^{+0.11}_{-0.09}$. This result allows us to
correct our composite emission-line luminosities for extinction due to dust,
yielding the typical intrinsic values of these lines in star-forming, dwarf
galaxies at high redshift. The dust-corrected (and stellar-absorption-
corrected in the case of the Balmer lines) emission-line luminosities of the
composite spectrum, relative to H$\beta$, are listed in Table 2.

When considering SFRs, similar to how our stellar masses are being reported,
we calculate the SFR for each individual galaxy and report the composite SFR
as the median value of the sample. These SFRs are calculated with slit-loss-
corrected, dust-corrected, H$\alpha$ luminosities ($L_{\rm{H\alpha}}$), de-
magnified according to the lensing models in Section 2.1.1. The dust-
extinction correction of each $L_{\rm{H}\alpha}$ value is conducted with the
“typical” extinction estimate for the sample found via the composite spectrum.
The $L_{\rm{H}\alpha}$ values are not corrected for stellar absorption, which
on average would result in an increase of $<1\%$. The $L_{\rm{H\alpha}}$
values are converted to SFRs using Equation 1 below:

$\mathrm{SFR}\ (\mathrm{M_{\odot}\ yr^{-1}})=4.645\times 10^{-42}\
L_{\mathrm{H}\alpha}\ (\mathrm{ergs\ s^{-1}})$ (1)

This equation is of the same form as the relation in Kennicutt (1998) for
calculating SFRs from recombination lines. However, the conversion factor here
has been recalculated assuming a metallicity of 0.2 $\rm{Z_{\odot}}$ and a
Chabrier (2003) IMF. Our estimates for the SFRs of the individual galaxies are
given in Table 1, and the SFR associated with the composite spectrum is taken
to be the sample-median value of $\rm{SFR_{med}}=2.25\
\rm{M_{\odot}\,yr^{-1}}$.

###  3.5 Electron Temperature and Electron Density

The “direct” calculation of metallicity relies on collisionally-excited oxygen
emission lines and the nebular properties of electron temperature ($T_{e}$)
and electron density ($n_{e}$), which are responsible for the strength of the
collisionally-excited lines. Electron temperature is calculated in two
ionization zones of the star-forming, H II regions. In the $\rm{O^{++}}$ zone,
$T_{e}$([O III]) is calculated using the electron-temperature-sensitive
emission-line ratio [O III] $\lambda\lambda$ 4959, 5007/[O III] $\lambda$4363
and the IRAF routine NEBULAR.TEMDEN (Shaw & Dufour, 1994). We note that while
$T_{e}$([O III]) does have a dependence on electron density, $n_{e}$, below
$n_{e}\approx 10^{3}\ \rm{cm^{-3}}$, $T_{e}$([O III]) is insensitive to
$n_{e}$ (Osterbrock & Ferland, 2006; Izotov et al., 2006) and can be
calculated assuming the typical $z\sim 2$ H II region electron density of a
few hundred per cubic centimeter (Sanders et al., 2016b). For our calculation,
we assume $n_{e}=150\ \rm{cm^{-3}}$ and obtain an electron temperature in the
$\rm{O^{++}}$ region of $T_{e}$([O III]) = 15,800 $\pm$ 3,100 K. We note in
regard to the assumed $n_{e}$ that Gburek et al. (2019) studied a galaxy in
our present stacking sample, A1689-217, that had a similar electron
temperature of $T_{e}$([O III]) = 14,300 K. When calculating this temperature,
they found that assuming any $n_{e}<10^{3}\ \rm{cm^{-3}}$ changed their result
by $<0.5\%$, suggesting our current assumption is robust.

Ideally, the electron temperature in the $\rm{O^{+}}$ ionization region is
calculated using measurements of the [O II] $\lambda\lambda$7320, 7330 auroral
emission-line doublet. Unfortunately, for the galaxies in our stacking sample,
we do not have spectroscopic coverage of these lines. Instead, we calculate
$T_{e}$([O II]) via the $T_{e}$([O III])$-T_{e}$([O II]) relation of Campbell
et al. (1986), reprinted here in Equation 2:

$T_{e}\text{([{O~{}II}])}=0.7\,T_{e}\text{([{O~{}III}])}+3000\ \text{K}$ (2)

Use of this equation gives us an electron temperature in the $\rm{O^{+}}$
region of $T_{e}$([O II]) = 14,100 $\pm$ 2,200 K.

The electron density, $n_{e}$, can be derived with the doublet ratio [O II]
$\lambda$3729/[O II] $\lambda$3726, the $T_{e}$([O II]) electron temperature,
and the IRAF routine NEBULAR.TEMDEN. For our composite spectrum, we calculate
[O II] $\lambda$3729/[O II] $\lambda$3726 = 1.56 $\pm$ 0.32 (1.51 $\pm$ 0.12
when using the statistical-only uncertainty spectrum and associated fits).
This corresponds to an electron density of $n_{e}=1^{+152}_{-0}\ \rm{cm^{-3}}$
($n_{e}=1^{+49}_{-0}\ \rm{cm^{-3}}$), where the “best-fit” $n_{e}$ value is
set to the low-density limit of $n_{e}=1\ \rm{cm^{-3}}$ as a result of the
best-fit [O II] ratio exceeding the maximum theoretical bound of [O II]
$\lambda$3729/[O II] $\lambda$3726 $\lesssim 1.5$ (Osterbrock & Ferland, 2006;
Sanders et al., 2016b).

We note that we have significantly detected the component emission lines of
the [O II] $\lambda\lambda$3726, 3729 doublet and resolved their individual
peaks (S/N(3726, 3729) = (6.4, 7.8); see Figure 1). The electron density
associated with our dwarf galaxy sample and with the best-fit ratio of these
lines is significantly lower than the densities found in more massive galaxies
at $1.5\lesssim z\lesssim 2.5$, which typically lie in the range of
$n_{e}\approx 100-300\ \rm{cm^{-3}}$ (Steidel et al., 2014; Sanders et al.,
2016b; Kashino et al., 2017; Kaasinen et al., 2017; Davies et al., 2021).

Table 3: Properties of the Dwarf Galaxy Composite

Property |  Value
---|---
$z_{\rm{mean}}\tablenotemark{a}$ |  2.30
$\log(M_{\ast}/\rm{M_{\odot}})^{\rm{med}}_{\rm{fiducial}}\tablenotemark{a,b}$ |  8.29
$\log(M_{\ast}/\rm{M_{\odot}})_{\rm{med}}\tablenotemark{a,c}$ |  8.92
$\rm{SFR}^{\rm{med}}_{\rm{H}\alpha}\ (\rm{M_{\odot}\ yr^{-1}})\tablenotemark{a}$ |  2.25
$E(\bv)_{\rm{gas}}$ |  $0.14^{+0.11}_{-0.09}$
$n_{e}^{\rm{boot}}\ (\rm{cm^{-3}})\tablenotemark{d}$ |  $1^{+152}_{-0}$
$n_{e}^{\rm{stat}}\ (\rm{cm^{-3}})\tablenotemark{d}$ |  $1^{+49}_{-0}$
$T_{e}$([O II]) (K) |  14,$100\pm 2$,200
$T_{e}$([O III]) (K) |  15,$800\pm 3$,100
$12+\log(\rm{O^{+}/H^{+}})$ |  $7.27^{+0.28}_{-0.19}$
$12+\log(\rm{O^{++}/H^{+}})$ |  $7.71^{+0.27}_{-0.20}$
$12+\log(\rm{O/H})_{\rm{direct}}$ |  $7.87^{+0.24}_{-0.22}$
$Z\ (\rm{Z_{\odot}})$ |  $0.15^{+0.11}_{-0.06}$

bbfootnotetext: Our fiducial median stellar mass assuming constant SFHs.

ccfootnotetext: The median stellar mass assuming non-parametric SFHs. See
Section 4.3.2.

Note. — All uncertainties here (except for $n_{e}^{\rm{stat}}$) derive from
the composite bootstrapped error spectrum. See Figure 1.

###  3.6 Direct Oxygen Abundance

We directly calculate the oxygen abundance, or gas-phase metallicity, of our
composite spectrum using the ionic abundance equations of Izotov et al.
(2006). These equations utilize the values of $T_{e}$([O II]), $T_{e}$([O
III]), and $n_{e}$ given in the preceding section as well as the dust-
corrected emission-line ratios of [O II] $\lambda\lambda$3726, 3729/H$\beta$
and [O III] $\lambda\lambda$4959, 5007/H$\beta$. We assume that the total
oxygen abundance is the summation of the ionic abundances in the H II region
$\rm{O^{+}}$ and $\rm{O^{++}}$ ionization zones as seen in Equation 3. Any
higher ionization states of oxygen are deemed to have a negligible
contribution to the metallicity.

$\mathrm{\frac{O}{H}\ \approx\ \frac{O^{+}}{H^{+}}+\frac{O^{++}}{H^{+}}}$ (3)

From our composite spectrum, we report a typical gas-phase metallicity for
high-redshift, star-forming, dwarf galaxies of
$12+\log(\rm{O/H})=7.87^{+0.24}_{-0.22}$ ($0.15^{+0.11}_{-0.06}\
\rm{Z_{\odot}}$; Asplund et al., 2021). This metallicity estimate, as well as
the calculations from our composite spectrum of the other physical properties
detailed in Section 3, are summarized in Table 3.

##  4 Results and Discussion

This section of the manuscript will take the measurements derived in the
previous section from our dwarf galaxy sample and composite spectrum and
analyze them in the context of strong-line abundance diagnostics and global
galaxy scaling relations. Prior to this, however, it is crucial to look at our
sample and stack compared to the broader galaxy population at $1.7\lesssim
z\lesssim 2.6$ in order to assess how representative our sample is of typical,
dwarf galaxies at this epoch.

###  4.1 How Representative is our Sample?

Here we will consider two main diagnostics, the [N II] Baldwin-Phillips-
Terlevich diagram (N2-BPT; Baldwin et al., 1981) and the $M_{\ast}-\rm{SFR}$
relation, or “star-forming main sequence.”

####  4.1.1 N2-BPT Diagnostic Diagram

In Figure 2, we show the location of our $\langle z\rangle=2.30$ stack of
star-forming, dwarf galaxies on the [O III] $\lambda$5007/H$\beta$ vs. [N II]
$\lambda$6583/H$\alpha$ BPT diagnostic diagram. This diagram is a useful tool
for distinguishing between star-forming galaxies (SFGs) and AGN through
optical strong-line ratios (Kewley et al., 2001; Kauffmann et al., 2003;
Kewley et al., 2013, K01, K03, K13, respectively), all without needing to
apply a dust correction to the line fluxes. The SFG locus of the BPT is also a
probe of changing physical conditions in star-forming regions with redshift
(e.g., Kewley et al., 2013, Figure 2). This has been an active area of
research in numerous high-redshift statistical studies (e.g., Steidel et al.,
2014; Shapley et al., 2015, S14, S15, respectively) which have shown that
high-$z$ SFGs cluster around a locus offset toward higher [O III]
$\lambda$5007/H$\beta$ and/or [N II] $\lambda$6583/H$\alpha$ when compared to
the star-forming locus of $z\sim 0$ SDSS galaxies. We show this in Figure 2,
where the SFG locus of $z\sim 0$ SDSS galaxies is given by the solid red line
(Kewley et al., 2013), and the offset SFG loci of the $z\sim 2.3$ KBSS-MOSFIRE
and MOSDEF surveys are displayed by the purple (Steidel et al., 2014) and
magenta (Shapley et al., 2015) lines, respectively. Recent work by Runco et
al. (2022) has shown that these high-$z$ SFG loci actually converge with
consistent emission-line-fitting applied to each sample and reaffirm the
existence of an offset in the BPT between local and high-$z$ SFGs.

We show, via the black star in Figure 2, that our stack of $z\sim 2.3$ star-
forming galaxies is also offset from the SDSS SFG locus, lying in parameter
space consistent with KBSS-MOSFIRE and MOSDEF. While we cannot place each
individual galaxy in our sample on this plot due to skyline contamination,
particularly of H$\beta$, we note that our stack lies below the “maximum
starburst” demarcation (dotted black line) of Kewley et al. (2001), the
empirical demarcation (dashed brown line) between SFGs and AGN of Kauffmann et
al. (2003), and the theoretical, $z=2.30$, upper-limit SFG locus (red dotted-
dashed line) of Kewley et al. (2013). We also note that the large uncertainty
of our stack in log([N II] / H$\alpha$) is the result of [N II] $\lambda$6583
only being detected in the composite spectrum with 2$\sigma$ significance.

Figure 2: The [O III] $\lambda$5007/H$\beta$ vs. [N II]
$\lambda$6583/H$\alpha$ BPT diagnostic diagram. Our $\langle z\rangle=2.30$
dwarf galaxy composite is shown as the black star. The $z\sim 2.3$
$M_{\ast}$-binned stacks of star-forming (SF) MOSDEF galaxies from Sanders et
al. (2021, S21) are shown as gray squares. The median $M_{\ast}$ of each of
these stacks, and of our composite (in bold), is listed in the upper right-
hand corner of the plot by order of increasing [N II]/H$\alpha$ of the stacks,
highlighting that our stack is an extension to lower $M_{\ast}$ (and O/H via
the MZR) of the MOSDEF survey. The $z\sim 2.3$ SF sequences of the MOSDEF
(Shapley et al., 2015, S15) and KBSS-MOSFIRE (Steidel et al., 2014, S14)
surveys are shown by the magenta and purple lines, respectively. Like these
sequences and the MOSDEF stacks, our composite also lies offset from the
$z\sim 0$ mean SF sequence given by the red line and parameterized by Kewley
et al. (2013, K13). We plot the $z=2.30$ SF sequence upper-limit from K13 as
the red dot-dashed line. The demarcation between SF galaxies and AGN of
Kauffmann et al. (2003, K03) is given by the dashed brown line, and the
“maximum starburst” curve of Kewley et al. (2001, K01) is given by the dotted
black line.

While the [N II] $\lambda$6583 measurement in our composite is fairly
uncertain, the location of our stack along the $x$-axis of the N2-BPT is
interesting when compared to the $M_{\ast}$-binned stacks of SFGs (the blue-
gray squares) from the MOSDEF survey and Sanders et al. (2021, S21). This is
because of the monotonic relationship that exists between log([N II]
$\lambda$6583/H$\alpha$) and metallicity; as this strong-line ratio increases,
metallicity increases (Pettini & Pagel, 2004; Maiolino et al., 2008; Curti et
al., 2017; Bian et al., 2018; Sanders et al., 2021). By the mass-metallicity
relation (e.g., Tremonti et al., 2004; Sanders et al., 2021), as this ratio
increases, the stellar mass of galaxies should then also increase on average.
We see this with the MOSDEF stacks, where in the upper right-hand section of
the plot, we list the median stellar masses of the stacks in order of
increasing [N II]/H$\alpha$. Here in this list we have also included, in bold,
the median stellar mass of our dwarf galaxy stack, which, based on its
positioning in the BPT, predictably has the lowest listed median stellar mass.
In this, we show that our sample is a complementary extension in stellar mass
to the MOSDEF (and KBSS-MOSFIRE) survey, extending its mass range into the
dwarf galaxy regime ($M_{\ast}<10^{9}\ \rm{M_{\odot}}$).

####  4.1.2 The Star-Forming Main Sequence

One of the primary goals of this study is to analyze our dwarf galaxy
composite relative to the $z\sim 2.3$ mass-metallicity relation (MZR) of star-
forming galaxies (see Section 4.3), which is a scaling relation between galaxy
stellar mass ($M_{\ast}$) and gas-phase oxygen abundance, or metallicity
(O/H). However, in order to properly contextualize our findings in relation to
the broader $z\sim 2.3$ dwarf galaxy population, we must consider the SFRs
associated with the stack and the sample that comprises it. This is due to the
existence of the fundamental metallicity relation (FMR) between $M_{\ast}$,
SFR, and O/H which has been demonstrated locally at $z\sim 0$ (e.g., Mannucci
et al., 2010, 2011; Lara-López et al., 2010; Andrews & Martini, 2013; Curti et
al., 2020) and at high redshift out to $z\sim 3.3$ (e.g., Henry et al., 2013a,
b; Cresci et al., 2019; Sanders et al., 2018, 2021). The FMR demonstrates
that, at a fixed $M_{\ast}$, a galaxy with an above-average (below-average)
SFR will typically have a below-average (above-average) metallicity.
Therefore, if our stacking sample (and thus composite spectrum) is biased in
SFR, it will not have an average metallicity representative of typical dwarf
galaxies at $z\sim 2.3$. This would be problematic when comparing the
metallicity of our composite to the low-mass end of the MZR.

To investigate whether our stacking sample has a bias in SFR, we plot our
sample and its median values against the $M_{\ast}-\rm{SFR}$ star-forming main
sequence (SFMS) in Figure 3. Galaxies that lie on this relation, which is
redshift-dependent, are considered to be representative of the typical galaxy
at that corresponding stellar mass and redshift. In Figure 3, we compare to
the $z\sim 2.3$ SFMS parameterizations of Sanders et al. (2021) and Whitaker
et al. (2014).

Figure 3: $M_{\ast}$ vs. dust-corrected $L_{\rm{H}\alpha}$ on the left-hand
axis and SFR on the right-hand axis. Our stacking sample is shown by the
colored circles, with blue circles representing our four $z\sim 1.7$ galaxies
and red circles our twelve $2.4\lesssim z\lesssim 2.6$ galaxies. The purple
star lies at the median $M_{\ast}$ and $L_{\rm{H}\alpha}$ (or SFR) of our
$z\sim 2.3$ sample. The black squares show the $z\sim 2.3$ $M_{\ast}$-binned
stacks of MOSDEF galaxies from Sanders et al. (2021). We compare our sample
median values against the $M_{\ast}-\rm{SFR}$ relations, or “star-forming main
sequences” (SFMS), of Sanders et al. (2021, black line) and Whitaker et al.
(2014, orange and green lines) to determine how representative our sample is
of typical, $z\sim 2.3$ dwarf galaxies. The Sanders et al. (2021) relation is
a power-law derived from H$\alpha$ SFRs. The orange and green relations of
Whitaker et al. (2014) are parameterized by a broken power-law and second-
order polynomial, respectively; both relations derive from UV+FIR SFRs. Dashed
portions of the SFMSs are extrapolations.

In comparing to the $z\sim 2.3$ SFMS of Sanders et al. (2021), we do so with
stellar masses and SFRs calculated in a manner highly consistent with the
methodologies adopted in Sanders et al. The stellar masses of our stacking
sample and the MOSDEF galaxies used to calibrate the SFMS both rely on
emission-line-corrected photometry and the SED-fitting (with FAST; Kriek et
al., 2009) assumptions of constant star formation histories, the Calzetti et
al. (2000) attenuation curve, and the Chabrier (2003) IMF. The SFRs in both
studies are calculated from dust-corrected (via the Cardelli et al. 1989
extinction curve) H$\alpha$ luminosities ($L_{\rm{H}\alpha}$). We note that
the conversion factor between $L_{\rm{H}\alpha}$ and SFR (see, for example,
Equation 1) is dependent on assumptions such as the stellar metallicity and
can vary between different authors. Therefore, we plot dust-corrected
$L_{\rm{H}\alpha}$ on the left-hand axis of Figure 3 so that our stacking
sample (red and blue circles) and the Sanders et al. (2021) MOSDEF stacks
(black squares) can be directly compared without the additional SFR
conversion. Ultimately, however, we find that the SFR conversion factor used
by both studies is very similar, and we continue our analysis of how
representative our stacking sample is via SFR, given on the right-hand axis of
Figure 3.

The $z\sim 2.3$ SFMS of Sanders et al. (2021) is parameterized as a power-law
over the stellar mass range 9.0 $<$ log($M_{\ast}$/$\rm{M_{\odot}}$) $<$ 11.0.
In Figure 3, we plot their best-fit relation141414The fitting of the SFMS was
done with the four lowest-mass bins. Additionally, there is evidence that the
lowest-mass ($M_{\ast}<10^{9.5}\ \rm{M_{\odot}}$) bin may be biased high in
SFR (Shivaei et al., 2015; Sanders et al., 2021). over this range as a solid
black line and extrapolate into the dwarf galaxy regime as seen by the dashed
black line. In plotting our stacking sample in the $M_{\ast}-\rm{SFR}$
parameter space, we differentiate the $z\sim 1.7$ galaxies from the
$2.4\lesssim z\lesssim 2.6$ galaxies by blue and red points, respectively. We
note that three of the four $z\sim 1.7$ galaxies lie along or above the
extrapolation of the $z\sim 2.3$ SFMS. These objects are biased high in SFR
relative to typical $z\sim 1.7$ galaxies that have a lower SFR at fixed
$M_{\ast}$ due to the redshift evolution of the SFMS (Speagle et al., 2014;
Whitaker et al., 2014; Sanders et al., 2021). The other group of $z>2.38$
galaxies are found to scatter around the Sanders et al. extrapolation.
Considered together, at the median mass of the complete $\langle z\rangle=2.3$
stacking sample, log($M_{\ast}$/$\rm{M_{\odot}}$)med = 8.29, the median SFR of
the sample, log(SFR/$\rm{M_{\odot}\,yr^{-1}}$)med = 0.35, lies
$\Delta\rm{log(SFR)}\sim 0.19$ dex above the extrapolation of the $z\sim 2.3$
SFMS of Sanders et al. (2021). This median point is shown as the purple star.

With the relative offset in SFR of our stacking sample in hand, we estimate
the bias in O/H, resulting from the FMR, of our composite spectrum.
Considering the strength of the SFR dependence of direct-method O/H at fixed
$M_{\ast}$ from Sanders et al. (2020),

$\Delta\rm{log(O/H)}\approx-0.29\times\Delta\rm{log(SFR\,/\,M_{\odot}\,yr^{-1})}$
(4)

our sample stack is biased by $\Delta\rm{log(O/H)}\sim-0.06$ dex, a value half
the statistical uncertainty of our composite direct-method metallicity
estimate ($\sigma_{\rm{stat}}\approx 0.12$ dex). We therefore conclude that,
when comparing to the SFMS of Sanders et al. (2021), our stacking sample of
dwarf galaxies does not have a major bias in SFR or O/H and, on average, is
representative of typical dwarf galaxies at $z\sim 2.3$ with $M_{\ast}\gtrsim
10^{8}\ \rm{M_{\odot}}$.

In Figure 3, we also plot the $2.0<z<2.5$ SFMS parameterizations of Whitaker
et al. (2014), which were fit to $M_{\ast}$-binned stacks above a mass-
completeness limit of $10^{9.2}\ \rm{M_{\odot}}$. Whereas Sanders et al.
(2021) fit the SFMS with a power-law, Whitaker et al. (2014) fit the SFMS with
both a second-order polynomial (green line) and a broken power-law (orange
line) for which a separate slope was fit above and below a characteristic mass
of log($M_{\ast}$/$\rm{M_{\odot}}$) = 10.2. We note that in our recreation of
the second-order polynomial fit, we use the more precise polynomial
coefficients given in the erratum (Whitaker et al., 2020) to Whitaker et al.
(2014, Table 1), in order to more accurately portray the curve. Similar to
Sanders et al. (2021) and the MOSDEF survey, the sample of Whitaker et al.
(2014) is composed of star-forming galaxies from the CANDELS fields (Grogin et
al., 2011; Koekemoer et al., 2011), though has a larger galaxy count and
different sample selection. Like with the MOSDEF galaxies and our stacking
sample, the stellar masses are determined with FAST assuming a Calzetti et al.
(2000) attenuation curve and Chabrier (2003) IMF; however, the star-formation
histories are taken to be exponentially declining. Unlike in Sanders et al.
(2021) and our stacking sample though, the SFRs are estimated from the
combination of rest-frame ultraviolet (UV) light and light re-radiated by dust
in the far-infrared (FIR).

Figure 3 shows that, at least qualitatively, the parameterizations from
Whitaker et al. (2014) generally agree with the power-law fit ($\beta=0.75$)
and MOSDEF stacks of Sanders et al. (2021), though begin to diverge as a
result of steeper slopes near unity below log($M_{\ast}$/$\rm{M_{\odot}}$)
$\lesssim 10$. This divergence, in part due to the choices of parametric form
made by each author, is particularly pronounced in the dwarf galaxy mass
regime where our sample lies. Unfortunately, there does not currently exist a
high-redshift statistical sample in this regime that would allow us to
confidently use a given parameterization. Therefore, we determine the most
applicable SFMS based on consistency with our study in physical property
estimation as well as the selection approach of our stacking sample. As
detailed above, our methodologies for estimating $M_{\ast}$ and SFR are most
analogous to those of Sanders et al. (2021), mitigating systematic
uncertainties between estimations of physical properties calculated with
different techniques. Additionally, when selecting our stacking sample, we did
not require a detection of [O III] $\lambda$5007, H$\alpha$, or [O III]
$\lambda$4363 which typically bias a sample toward higher SFR at fixed
$M_{\ast}$. Our sample selection instead suggests a more representative sample
like we see when comparing to the SFMS of Sanders et al. (2021). We therefore
conclude that the SFMS of Sanders et al. (2021) is the most applicable
comparison and that our stacking sample is representative of typical, star-
forming, dwarf galaxies at $z\sim 2.3$.

Figure 4: Strong emission-line ratios as a function of direct-method
metallicity. From left-to-right and top-to-bottom we consider the oxygen-based
strong-line ratios of $\rm{O_{32}}$, $\rm{O_{3}}$, $\rm{O_{2}}$ and
$\rm{R_{23}}$. Our $z\sim 2.3$ composite is displayed via the black star in
each panel. For reference, we also include the $z\sim 2.4$ KBSS-LM1 composite
of 30 star-forming galaxies also selected independent of emission-line
strength (gray diamond; Steidel et al., 2016), as well as the median values of
the $z\sim 2.2$ auroral-line sample of Sanders et al. (2020, gray square). We
compare these points to the locally-calibrated, strong-line metallicity
relations of Maiolino et al. (2008, orange line), Jones et al. (2015, green
line), Curti et al. (2017, cyan line), Bian et al. (2018, both the local
reference$\,-\,$purple line$\,-\,$and high-$z$ analog$\,-\,$red
line$\,-\,$relations), and Sanders et al. (2021, DIG-corrected; brown line).
We note the metallicity-insensitivity of the $\rm{O_{3}}$ and $\rm{R_{23}}$
indices at the oxygen abundances considered here. We find that our composite
and KBSS-LM1 have their metallicities best reproduced by the local reference
relations of Bian et al., in contrast to the auroral-line sample of Sanders et
al. that favors the high-$z$ analog relations.

###  4.2 Strong-Line Metallicity Calibrations at High-z

A major outstanding issue and active area of research in high-$z$ astronomy is
how to accurately calculate the gas-phase metallicities of the star-forming
galaxies in the various large, statistical, spectroscopic surveys at $z>1$ and
$M_{\ast}\gtrsim 10^{9}\,\rm{M_{\odot}}$ (e.g., 3D-HST; Brammer et al. 2012b,
KBSS-MOSFIRE; Steidel et al. 2014, MOSDEF; Kriek et al. 2015, FMOS-COSMOS;
Kashino et al. 2019). The cause of this problem is two-fold. For one, auroral
lines such as O III] $\lambda\lambda$1661, 1666 or [O III] $\lambda$4363,
needed for direct, $T_{e}$-based metallicity estimation, are exceedingly
faint, especially at high-redshift and with increasing galaxy stellar mass.
Additionally, it is not fully understood how changing physical conditions with
redshift in star-forming regions affect locally-calibrated, strong-line ratio
metallicity diagnostics for indirect metallicity estimation. In other words,
the accuracy and applicability of these strong-line metallicity calibrations
at high-redshift is an open question which several studies have tried to
address (Jones et al., 2015; Sanders et al., 2016a; Patrício et al., 2018;
Gburek et al., 2019; Sanders et al., 2020).

In Figure 4, we revisit this issue with our $z\sim 2.3$ dwarf galaxy
composite. Here we plot our stack (the black stars) in the parameter space of
$T_{e}$-based oxygen abundance versus various commonly-used, oxygen-based,
strong emission-line ratios. The dust-corrected emission-line ratios and
direct metallicity ($12+\log(\rm{O/H})=7.87^{+0.24}_{-0.22}$) of our stack are
measured (see Section 3) from the composite spectrum in Figure 1. Our stack in
these plots is compared to several locally-calibrated strong-line metallicity
diagnostics from the literature (Maiolino et al., 2008; Jones et al., 2015;
Curti et al., 2017; Bian et al., 2018; Sanders et al., 2021) in an attempt to
shed light on which strong-line ratios are serviceable at the typical
metallicity of our stacking sample as well as which calibrations most
favorably reproduce our composite metallicity at fixed strong-line ratio.
While we cannot comment on the shape or slope of the various calibrations, we
can get a sense of the appropriate normalization of these relations when
considering high-$z$ dwarf galaxies. In Figure 4, we consider four strong-line
ratios: $\rm{O_{32}}$ = log([O III] $\lambda$5007/[O II] $\lambda\lambda$3726,
3729), $\rm{O_{3}}$ = log([O III] $\lambda$5007/H$\beta$), $\rm{O_{2}}$ =
log([O II] $\lambda\lambda$3726, 3729/H$\beta$), and $\rm{R_{23}}$ = log(([O
III] $\lambda\lambda$4959, 5007 + [O II] $\lambda\lambda$3726,
3729)/H$\beta$). We do not consider strong-line metallicity calibrations based
on [N II] in this work due to our low detection significance ($\sim 2\sigma$)
of [N II] $\lambda$6583 in the composite spectrum as well as concerns in the
literature (e.g., Masters et al., 2014, 2016) of elevated N/O abundance ratios
at high-redshift.

Of immediate note when considering the location of our stack relative to the
$\rm{O_{3}}$\- and $\rm{R_{23}}$-based calibrations in the right-hand panels
of Figure 4 is that our stack lies at or near the apex of these relations in
what is called the “turnover” between the high and low metallicity branches of
these calibrations. In these turnover regimes, the strong-line ratio is
insensitive to the metallicity of a galaxy, giving these relations little
value as useful metallicity indicators over the oxygen abundance range spanned
by the turnover region (somewhere roughly between
$7.7\lesssim\rm{12+log(O/H)}\lesssim 8.3$ depending on the strong-line
index$\,-\,\rm{O_{3}}$ or $\rm{R_{23}}\,-\,$ and calibration used). We
therefore do not recommend the use of these strong-line indices for $z\sim 2$
dwarf galaxies similar to our those in our stacking sample. These results and
conclusion are not particularly surprising as [O III] $\lambda$4363-emitter
studies (Sanders et al., 2016a; Gburek et al., 2019; Sanders et al., 2020) and
photoionization modeling (e.g., Steidel et al., 2014) have shown that it is
quite common for $z\sim 2$ star-forming galaxies, over a couple orders of
magnitude in stellar mass, to have metallicities that lie within these
insensitive turnover regions.

A more interesting result is revealed when looking at the dwarf galaxy stack
relative to the $\rm{O_{32}}$\- and $\rm{O_{2}}$-based calibrations in the
left-hand panels of Figure 4. In particular, we focus on the comparisons with
the strong-line metallicity relations of Bian et al. (2018), who used stacked
SDSS spectra to create $T_{e}$-based, empirical metallicity calibrations over
the metallicity range $7.8<\rm{12+log(O/H)}<8.4$. These calibrations were
created from two distinct SDSS samples, a reference sample of galaxies lying
within $\pm 0.05$ dex of the $z\sim 0$ star-forming sequence of the N2-BPT
diagram (parameterized by Kewley et al. 2013) and a high-redshift analog
sample of SDSS galaxies lying within $\pm 0.04$ dex of the offset $z\sim 2.3$
star-forming sequence of the BPT diagram (Steidel et al., 2014). With these
selection criteria based on nebular emission-line ratios, the calibrations of
Bian et al. (2018) should represent the conditions of star-forming regions in
low- and high-redshift galaxies, respectively. We note that the Bian et al.
calibrations presented in Figure 4 were re-fit for this study for the reasons,
and via the methods, described in Appendix A.

When comparing our stack (black stars) to the strong-line metallicity
relations of Bian et al. (2018), we find that our stack favors the local
reference calibrations (purple curves), such that when considering the
$\rm{O_{32}}$ and $\rm{O_{2}}$ indices, the local reference relations
reproduce our composite metallicity to within $\lesssim 0.15$ dex at fixed
strong-line ratio. The high-$z$ analog relations of Bian et al. (red curves)
and the calibrations of Maiolino et al. (2008, orange curves), Jones et al.
(2015, green curves), and Curti et al. (2017, cyan curves) are all $\gtrsim
1\sigma$ inconsistent with our stack in the $\rm{O_{32}}$ and $\rm{O_{2}}$
panels, the lone exception being the $\rm{O_{32}}$ relation of Maiolino et al.
This result of inconsistency with these relations has added significance in
that we are analyzing a composite of 16 galaxies rather than an individual
source. Also, our uncertainties here take into account sample variance in
addition to statistical error (see Section 3.3). We note that, in all panels
other than that of the $\rm{O_{2}}$ index, our composite is also highly
inconsistent with the calibrations of Sanders et al. (2021, brown curves). At
metallicities below $\rm{12+log(O/H)}<8.4$, the Sanders et al. (2021)
relations are calibrated with the H II region spectra of dwarf galaxies from
Berg et al. (2012) and the Spitzer Local Volume Legacy survey (Dale et al.,
2009). These dwarf galaxies have been shown in previous studies to not follow
the other strong-line metallicity relations considered in this work, possibly
as a result of biases from the selection methods of the various calibration
samples or an incompleteness in low-metallicity, high-excitation H II regions
(Gburek et al., 2019; Sanders et al., 2020, 2021).

It is interesting that our $z\sim 2.3$ stack generally best agrees with the
local reference calibrations of Bian et al. (2018), particularly when compared
to the findings of Sanders et al. (2020), who compiled O III]
$\lambda\lambda$1661, 1666 and [O III] $\lambda$4363-detected sources at $z>1$
from the literature and the MOSDEF survey and conducted a similar study of
strong-line metallicity diagnostics at the median redshift of the compiled
sample, $z_{\rm{med}}\sim 2.2$. These authors found that at the median
metallicities and median line ratios of the galaxies in their $z>1$ sample
(gray squares in Figure 4; each median limited to galaxies with detections of
the respective line ratio), the high-redshift analog calibrations of Bian et
al. (2018) best reproduced their metallicities at fixed line-ratio. Moreover,
Sanders et al. (2020) found general agreement between their median points and
the Curti et al. (2017) calibrations as well as general agreement with the
calibration sample ($z\sim 0$ SDSS galaxies from Izotov et al. 2006) of the
Jones et al. (2015) relations. Consequently, when considering the results of
our work and those of Sanders et al. (2020), both studies at similar redshift
and metallicity, a tension exists in regard to the evolution with redshift of
strong-line metallicity diagnostics and therefore which calibrations are
reliable at high-$z$.

A likely source of the discrepancy and tension seen between the results of
this work and that of Sanders et al. (2020) lies in how each galaxy sample was
selected. For our dwarf galaxy stacking sample, we did not select galaxies
based on the strength of any particular rest-optical emission-line (see
Section 2.3). This is important particularly when considering [O III]
$\lambda$4363, [O III] $\lambda$5007 ([O III] $\lambda$4959 in our case), and
H$\alpha$. By avoiding selection based on line-strength, we mitigate biases in
our sample such as high sSFRs (SFR /$\,M_{\ast}$) and high excitation,
resulting in a sample which very nearly falls on an extrapolation of the mean
$z\sim 2.3$ $M_{\ast}-\rm{SFR}$ relation (see Section 4.1.2 and Figure 3). In
contrast, the $z\sim 2.2$ sample of Sanders et al. (2020) comprises galaxies
selected on auroral-line detection (either O III] $\lambda\lambda$1661, 1666
or [O III] $\lambda$4363, detected at $\rm{(S/N)_{med}}=6.0_{-2.5}^{+1.5}$)
for $T_{e}$-based metallicity estimation. Due to the faintness of these
auroral lines$\,-\,$[O III] $\lambda$5007 is $\sim 30-100\times$ brighter than
[O III] $\lambda$4363 (Jones et al., 2015)$\,-\,$this detection requirement
preferentially selects younger galaxies with high-excitation, highly-ionized
star-forming regions. These galaxies lie well above the $M_{\ast}-\rm{SFR}$
relation. In fact, the sample of Sanders et al. (2020) lies an average $\sim
0.6$ dex above the best-fit $z\sim 2.3$ $M_{\ast}-\rm{SFR}$ relation of
Sanders et al. (2018) (which is very similar to the fit of the same relation
in Sanders et al. 2021) and has $\rm{O_{32}}$ values $\sim 0.5$ dex higher on
average than typical $z\sim 2.3$ star-forming MOSDEF galaxies at fixed
$M_{\ast}$ (Sanders et al., 2020, Figure 16). Sanders et al. shows that this
auroral-line sample is not representative of typical $z\sim 2.3$ star-forming
galaxies, but rather coincident with $z\sim 2$ extreme emission-line galaxies
(EELGs). In effect, at roughly fixed O/H in Figure 4, we are comparing two
galaxy samples that differ significantly in the ionization state and extremity
of their star-forming regions. As such, both samples are not well-represented
by one single strong-line metallicity calibration.

When considering the sample selection methodologies for our stacked sample and
the $z>1$ auroral-line sample of Sanders et al. (2020) against those for the
calibration samples used to parameterize the relations shown in Figure 4, it
is perhaps not surprising that several relations appear to agree with the
Sanders et al. auroral-line sample while being $>1\sigma$ inconsistent with
our stacked sample. In the low-metallicity regime ($12+\log\rm{(O/H)}\lesssim
8.4$), the locally-calibrated strong-line metallicity relations are generally
defined by individual galaxies with [O III] $\lambda$4363 detections. This is
the case for the relations of Maiolino et al. (2008), Jones et al. (2015), and
Curti et al. (2017). This requirement of an [O III] $\lambda$4363 detection in
individual galaxies, as well as the BPT-related requirement of the Bian et al.
(2018) high-$z$ analog calibration (detailed above), selectively probe
galaxies with ISM conditions more extreme than in typical $z\sim 0$ star-
forming galaxies. Rather, ISM conditions more akin to those in the $z>1$
auroral-line sample are probed, leading to low-metallicity strong-line
calibrations that are likely biased high in $\rm{O_{3}}$ and $\rm{O_{32}}$
(and low in $\rm{O_{2}}$) and that closely predict the Sanders et al. (2020)
$z>1$ sample metallicity at fixed strong-line ratio. Meanwhile, the selection
criterion for the Bian et al. (2018) local reference calibrations (detailed
above) selectively probes galaxies with less extreme star-forming conditions
and better predicts the metallicity at fixed strong-line ratio of our dwarf
galaxy stack, which is not reliant on emission-line detections.151515A more
detailed study of the biases in certain locally-calibrated strong-line
metallicity relations (particularly of Curti et al. 2017), and how those
biases factor into the observed evolution of strong-line metallicity relations
with redshift, can be found in Sanders et al. (2020).

####  4.2.1 The KBSS-LM1 Composite of Steidel et al. (2016)

In this subsection, we briefly discuss the $z\sim 2.4$ composite spectrum of
Steidel et al. (2016) in the context of the strong-line ratio$\,-\,$direct
metallicity relations displayed in Figure 4. This composite, referred to in
Steidel et al. (2016) as “KBSS-LM1,” is derived from the rest-frame far-UV and
optical spectra of 30 star-forming galaxies from the KBSS-MOSFIRE
spectroscopic survey. The galaxies comprising KBSS-LM1 were notably selected
to have emission-line measurements or limits (not detections) of, among other
lines, [O II] $\lambda\lambda$3726, 3729, H$\beta$, [O III]
$\lambda\lambda$4959, 5007, and H$\alpha$. The sample of galaxies was also
selected such that it broadly represents the full KBSS-MOSFIRE sample in SFR,
$M_{\ast}$, and O/H, the latter of which was calculated via the
$T_{e}$-sensitive UV emission-line doublet, O III] $\lambda\lambda$1661, 1666.
The median $M_{\ast}$ and derived $T_{e}$-based metallicity of KBSS-LM1 are,
as reported by Sanders et al. (2020), log($M_{\ast}/\rm{M_{\odot}}$) = 9.8
$\pm$ 0.3 and 12+log(O/H) = 8.14 $\pm$ 0.03, respectively. With these
selection criteria, KBSS-LM1 is therefore similar to our $z\sim 2.3$ dwarf
galaxy composite in the sense that both samples should be representative of
typical $z\sim 2.3$ galaxies at their respective stellar masses. Indeed,
Sanders et al. (2020, Figure 16) find KBSS-LM1 to lie just above the $z\sim
2.3$ star-forming main sequence as well as amongst the typical $z\sim 2.3$
MOSDEF galaxies in $\rm{O_{32}}$ at fixed $M_{\ast}$.

In Figure 4, along with our $z\sim 2.3$ dwarf galaxy stack (black stars) and
$z\sim 2.2$ auroral-line sample of Sanders et al. (2020, gray squares), we
plot KBSS-LM1 as a gray diamond in each panel. As with our dwarf galaxy stack,
albeit at higher O/H owing to its higher $M_{\ast}$, we see that KBSS-LM1 is
best represented by the local reference calibrations of Bian et al. (2018)
instead of the high-redshift analog calibrations, though note the small
statistical error bars (not bootstrapped) of KBSS-LM1. This result suggests
that $z\sim 2.3$ star-forming galaxies with
$8.4\lesssim\log(M_{\ast}/\rm{M_{\odot}})\lesssim 9.8$ that lie on the
$M_{\ast}$-SFR relation will, on average, have their metallicities most
accurately predicted at fixed strong-line ratio via the local reference
strong-line metallicity calibrations of Bian et al. (2018). In contrast,
$z\sim 2.3$ EELGs and galaxies with more extreme ISM conditions, like the
$z\sim 2.2$ auroral-line sample, may require independently-calibrated strong-
line metallicity relations more akin to the high-$z$ analog calibrations of
Bian et al. (2018).

Figure 5: The stellar mass$\,-\,$gas-phase metallicity relation (MZR). Left:
We compare our $z\sim 2.3$ dwarf galaxy stack (green star) at our sample’s
median mass, log($M_{\ast}$/$\rm{M_{\odot}}$)med = 8.29, against the $z\sim
2.2$ direct-method MZR of Sanders et al. (2020, black line) and the
extrapolations of the $z\sim 2.3$ strong-line MZR of Sanders et al. (2021,
blue solid/dashed line) and the $z\sim 2.3$ photoionization-model-based MZR of
Strom et al. (2022, purple solid/dashed line). We find excellent agreement
between our stack and the direct-method MZR, the most comparable of the
relations. For reference, we also include the mean values of the $z>1$
auroral-line sample split into two $M_{\ast}$ bins (black squares; Sanders et
al., 2020), the MOSDEF galaxy stacks from Sanders et al. (2021, blue squares),
and the $z\sim 2.4$ KBSS-LM1 composite of Steidel et al. (2016, orange
diamond). Right: We compare our $z\sim 2.3$ stack (green star; same as in left
panel) to the predicted $z=2$ and $z=3$ MZRs (dark-red and light-red lines,
respectively) of the IllustrisTNG100 simulations (Torrey et al., 2019) and the
$z=2.3$ MZR of the FIRE simulations (red dot-dashed line; Ma et al., 2016),
finding consistency, within errors, with both sets of MZRs. When comparing our
stack after recalculating the sample stellar masses assuming more realistic
non-parametric SFHs (red unfilled square; log($M_{\ast}$/$\rm{M_{\odot}}$)med
= 8.92), we find excellent agreement with the MZR from FIRE. The MZRs from the
left-hand panel are recreated here for reference. We note that for each
display of our stack, the $x$-axis error bar corresponds to the $M_{\ast}$
range of our sample.

###  4.3 The Stellar Mass$\,-\,$Gas-Phase Metallicity Relation

This section explores our $z\sim 2.3$ dwarf galaxy sample in relation to the
low-mass end of the stellar mass$\,-\,$gas-phase metallicity relation (MZR).
The MZR, shown to exist both locally and at high redshift, displays a positive
correlation, parameterized as a power-law, between galaxy stellar mass
($M_{\ast}$) and gas-phase oxygen abundance (O/H) at lower stellar masses
before flattening asymptotically at higher masses ($M_{\ast}\gtrsim 10^{10.5}\
\rm{M_{\odot}}$ locally; Tremonti et al., 2004; Erb et al., 2006; Maiolino et
al., 2008; Andrews & Martini, 2013; Steidel et al., 2014; Curti et al., 2020;
Sanders et al., 2020, 2021; Strom et al., 2022). The MZR also evolves with
time, such that at higher redshifts, the average metallicity of star-forming
galaxies is lower at fixed $M_{\ast}$ than it is locally (Erb et al., 2006;
Maiolino et al., 2008; Zahid et al., 2013, 2014a, 2014b; Steidel et al., 2014;
Sanders et al., 2020, 2021). Here we compare our dwarf galaxy composite
against recent empirical and theoretical parameterizations of the MZR in an
effort to better constrain the low-mass slope and normalization of the
relation at $z\sim 2.3$.

####  4.3.1 Comparison to MZRs Derived from Observations

In the left-hand panel of Figure 5, we plot our $z\sim 2.3$ dwarf galaxy
composite against the $z\sim 2.2$ direct-method MZR of Sanders et al. (2020),
the $z\sim 2.3$ strong-line MZR of Sanders et al. (2021), and the $z\sim 2.3$
MZR of Strom et al. (2022), who calculated their metallicities via
photoionization modelling. The direct-method MZR shown here (black line) is
derived from the Sanders et al. (2020) $z>1$ ($z_{\rm{med}}\sim 2.2$) auroral-
line sample considered in Section 4.2 above. However, as previously mentioned,
this sample of galaxies lies $\langle\Delta\log(\rm{SFR})\rangle\sim 0.6$ dex
above the $z\sim 2.3$ $M_{\ast}\,-\,$SFR relation, which has the effect, via
the FMR, of yielding a MZR biased low in O/H relative to the MZR expected for
typical galaxies that fall on the star-forming main sequence at this redshift.
As such, the direct-method MZR displayed in Figure 5 was parameterized by
Sanders et al. (2020) after applying SFR-corrections to the $z>1$ sample
metallicities via Equation 4. Additionally, Sanders et al. (2020) adjusted
this MZR to account for a low-redshift bias in the low-mass $z>1$ galaxies.
These galaxies (low-mass black square; 5 galaxies) were found to have a median
redshift $\Delta z_{\rm{med}}\sim 0.5$ lower than the high-mass sample
galaxies (high-mass black square; 9 galaxies), and were estimated, on average,
to have a metallicity biased 0.1 dex high due to the redshift evolution of O/H
at fixed $M_{\ast}$ ($d$log(O/H)/$dz\approx-0.2$; Sanders et al., 2020). With
these adjustments, the direct-method MZR reproduced in Figure 5 is an
estimation of the MZR$\,-\,$on the $T_{e}$-based abundance scale$\,-\,$of
typical galaxies lying on the star-forming main sequence at a redshift of
$z\sim 2.2$.

The blue line (and its extrapolation) in the left-hand panel of Figure 5 is a
recreation of the $z\sim 2.3$ strong-line MZR of Sanders et al. (2021), fit as
a power-law to the four lowest-mass bins of MOSDEF galaxies displayed here as
blue squares (The highest-mass bin suffers from incompleteness and was
excluded from the fit (Sanders et al., 2021, Section 2.4)). These MOSDEF
stacks were also used by Sanders et al. to parameterize their
$M_{\ast}-\rm{SFR}$ relation (see Section 4.1.2 and Figure 3). The
metallicities of the MOSDEF stacks were calculated via the high-$z$ analog
calibrations of Bian et al. (2018) and the $\alpha$-element-based strong-line
ratios of $\rm{O_{32}}$, $\rm{O_{3}}$, and log([Ne III] $\lambda$3869/[O II]
$\lambda\lambda$3726, 3729). These relations were chosen so as to use
calibrations that most closely reproduce the excitation sequences, and thus
likely ISM conditions, of the MOSDEF $z\sim 2.3$ star-forming sample.

The $z\sim 2.3$ MZR of Strom et al. (2022), and its extrapolation, are shown
as the purple line in the left-hand panel of Figure 5. This relation was fit
to 195 individual star-forming galaxies from the KBSS survey. The
metallicities in this study were estimated with photoionization models,
described in Strom et al. (2018, 2022), that are able to reproduce the rest-UV
and rest-optical spectroscopic properties of high-$z$, star-forming galaxies.

When comparing our $z\sim 2.3$ dwarf galaxy stack (green star; $x$-axis error
bar represents the $M_{\ast}$ range of the stacking sample) against the
direct-method and strong-line MZRs of Sanders et al. (2020) and Sanders et al.
(2021), respectively, we find that at the median stellar mass of our stacking
sample, log($M_{\ast}$/$\rm{M_{\odot}}$)med = 8.29, and the direct-method
metallicity calculated from the composite spectrum,
$12+\log(\rm{O/H})=7.87^{+0.24}_{-0.22}$, the stack lies virtually on top of
the $z\sim 2.2$ direct-method MZR. Within uncertainties, our stack (and the
mean values of the low- and high-mass bins of the $z>1$ auroral-line sample;
Sanders et al., 2020, black squares) is also consistent with the extrapolation
of the strong-line MZR, lying below this extrapolation by $\sim 0.13$ dex.
Correcting the metallicity of our composite for the slightly-high SFR bias of
the stacking sample (increasing log(O/H) by $\sim 0.06$ dex; see Section
4.1.2) reduces the offset of the stack from the extrapolation of the strong-
line MZR while moving the stack $\sim 0.06$ dex above the direct-method MZR.
This consistency with both MZRs supports their fit power-law parameters, in
particular the slope, which was fit by Sanders et al. (2020) to be
$\beta=0.37$ and by Sanders et al. (2021) to be $\beta=0.30$. However, while
we are comparing to the strong-line MZR of Sanders et al. (2021) due to the
authors’ careful consideration in selecting applicable strong-line
calibrations for high-$z$ star-forming regions, we note that systematic
uncertainties still exist between metallicities calculated directly versus
with strong-line proxies. Therefore, in constraining the slope of the $z\sim
2.3$ MZR, our results and direct-method of metallicity estimation most
favorably suggest the slope fit to the direct-method MZR of Sanders et al.
(2020), $\beta=0.37$.

When considering the $z\sim 2.3$ photoionization-model-based MZR of Strom et
al. (2022), we find that our stack lies $\approx 1\sigma$ below the
extrapolation of this relation, in disagreement with the shallow slope
($\beta=0.14$) proposed by Strom et al. However, we must again take into
account the difference in methods of abundance estimation between the two
studies. As discussed in the literature (e.g., Kewley & Ellison, 2008;
Maiolino & Mannucci, 2019), photoionization models typically overestimate
metallicities by $\sim 0.2$ dex or more compared to direct metallicities. In
general, this is due to a poorly-constrained combination of factors, such as
photoionization models accounting for dust depletion or $T_{e}$-based
metallicities potentially being biased low due to temperature fluctuations or
gradients in star-forming regions leading to nebular spectra dominated by
brighter auroral lines from high-$T_{e}$ zones. While the discrepancy between
direct and theoretical metallicity estimates is stronger at higher
metallicities (e.g., Kewley & Ellison, 2008; Maiolino & Mannucci, 2019; Curti
et al., 2020), Steidel et al. (2016) found that photoionization models predict
an oxygen abundance 0.25 dex higher than their $T_{e}$-based estimate of
12+log(O/H) = 8.14 (0.29 $\rm{Z_{\odot}}$) for KBSS-LM1, a composite of 30
star-forming KBSS galaxies that we briefly discussed in Section 4.2.1. If we
apply this same offset to the $T_{e}$-based metallicity of our dwarf galaxy
composite, our stack will lie very near the extrapolation of the Strom et al.
(2022) MZR. This said, there are too many systematic uncertainties involved to
accurately compare our composite against this MZR at present, and we conclude
that the best estimate for the slope of the $z\sim 2.3$ MZR is $\beta=0.37$,
given by the direct-method MZR of Sanders et al. (2020).

We note that the stellar masses of the samples and studies considered here are
calculated via SED-fitting with consistent assumptions in star-formation
history, IMF, and extinction law. We also note that the KBSS-LM1 composite of
Steidel et al. (2016) lies significantly below the displayed MZRs, including
the direct-method MZR, when plotted at its reported direct metallicity of
12+log(O/H) = $8.14\pm 0.03$. However, we do acknowledge the very small
statistical-only uncertainty of this metallicity.

####  4.3.2 Comparison to MZRs from Cosmological Simulations

In the right-hand panel of Figure 5, we now compare our $z\sim 2.3$ dwarf
galaxy stack against predicted MZRs from well-known cosmological simulations,
the Feedback in Realistic Environments simulations161616FIRE:
https://fire.northwestern.edu/ (FIRE; Hopkins et al., 2014) and The Next
Generation Illustris simulations171717IllustrisTNG: https://www.tng-
project.org/ (IllustrisTNG; Weinberger et al., 2017; Pillepich et al., 2018),
which are the successor to the Illustris simulation (Vogelsberger et al.,
2014a, b; Genel et al., 2014). Our stack, with the median mass
(log($M_{\ast}$/$\rm{M_{\odot}}$) = 8.29) and direct metallicity
($12+\log(\rm{O/H})=7.87^{+0.24}_{-0.22}$) as reported throughout this paper,
is shown by the green star, with the $x$-axis error bar representing the
$M_{\ast}$ range of our stacking sample. The predicted MZR from FIRE, given by
the redshift-dependent, fixed-slope ($\beta=0.35$), gas-phase MZR fitting
function in Ma et al. (2016), is evaluated at $z=2.30$ (the mean redshift of
our stacking sample) and is represented by the red dot-dashed line. The $z=2$
and $z=3$ MZRs from IllustrisTNG (Torrey et al., 2019), specifically from the
TNG100 simulation volume, are displayed by the solid dark-red and light-red
lines, respectively. For reference to MZRs derived from observations, in the
right-hand panel we reproduce, from the left-hand panel, the $z\sim 2.2$
direct-method MZR of Sanders et al. (2020, black line), the $z\sim 2.3$
strong-line MZR of Sanders et al. (2021, blue line), and the $z\sim 2.3$
photoionization-model-based MZR of Strom et al. (2022, purple line).

Properly assessing the predicted MZRs considered in this section requires a
reliable and accurate empirical metallicity estimation method such as the
direct-method, which estimates oxygen abundances through directly probing
physical properties ($T_{e}$ and $n_{e}$) of star-forming regions. We find
that at fixed $M_{\ast}$, the direct metallicity of our stack (green star) is
consistent within uncertainties with both the $z=2.3$ FIRE MZR and $2\leqslant
z\leqslant 3$ MZR of IllustrisTNG. We also observe that the $z\sim 2.2$
direct-method MZR agrees favorably in slope and normalization with the
IllustrisTNG MZR above $M_{\ast}\approx 10^{9}\ \rm{M_{\odot}}$. Below this
mass, the IllustrisTNG MZR deviates away from this slope and the direct-method
MZR toward higher metallicity values, displaying a “bump” in the low-mass
regime. This bump is the result of the minimum wind velocity
($v_{\rm{min}}=350\ \rm{km}\ \rm{s^{-1}}$) enforced in the stellar feedback
models of IllustrisTNG, with $v_{\rm{min}}$ put in place so that the
simulations match the low-end of the galaxy stellar mass function (Pillepich
et al., 2018). As Torrey et al. (2019) explain, while $v_{\rm{min}}$ is not
directly a function of $M_{\ast}$, it is generally set as the wind velocity in
galaxies with $M_{\ast}\lesssim 10^{9}\ \rm{M_{\odot}}$ due to the low dark
matter velocity dispersions in their halos. The higher metallicities at these
masses that we see as the bump in the MZR then arise because these fixed-
velocity winds eject less gas, and therefore fewer metals, than would be the
case if the wind velocity were allowed to be $v<350\ \rm{km}\ \rm{s^{-1}}$.
While further observations are needed to either confirm or deny this bump in
the MZR, as well as evaluate the applicability of the minimum wind velocity
assumption, the direct-method MZR of Sanders et al. (2020) suggests that such
a bump likely does not exist and that wind velocities in low-mass galaxies can
extend lower than $v<350\ \rm{km}\ \rm{s^{-1}}$. This suggestion is reinforced
below when revisiting our dwarf galaxy stacking sample with stellar masses
recalculated under more realistic assumptions.

In this work, in order to facilitate fair comparisons of our empirical results
to those in the literature, when estimating the stellar masses of our stacking
sample via SED-fitting, we made assumptions consistent with those generally
found in the literature; in particular, we assumed constant star formation
histories (SFH). However, the SFHs of galaxies, particularly of galaxies at
high-$z$ owing to their higher-EW emission lines, are likely not well-
described by such simple parameterizations. Instead, stellar masses are likely
more accurate if calculated assuming non-parametric SFHs which can better
reveal the presence of older stellar populations that are hidden in the rest-
UV and rest-optical by brighter, younger stars (Gburek et al., 2019; Tang et
al., 2022; Whitler et al., 2022; Topping et al., 2022). Having these more
realistic, typically larger stellar masses is important when comparing to
simulation results. Therefore, we recalculated the stellar masses of our
stacking sample, via the
PROSPECTOR181818https://prospect.readthedocs.io/en/latest/,191919https://github.com/bd-j/prospector
SED-fitting code (Johnson et al., 2021), assuming non-parametric SFHs. The
resultant median stellar mass of our sample under these assumptions becomes
log($M_{\ast}$/$\rm{M_{\odot}}$) = 8.92, an increase to our fiducial median
stellar mass of 0.63 dex. Our dwarf galaxy stack, shifted to this recalculated
median $M_{\ast}$, is shown in Figure 5 as the red unfilled square. As with
the $x$-axis error bar of our fiducial point (the green star), the $x$-axis
error bar of the red square represents the recalculated range of $M_{\ast}$ in
our stacking sample.

When comparing our $z\sim 2.3$ dwarf galaxy stack to the simulated MZRs
following recalculating the sample stellar masses under the assumption of more
realistic non-parametric SFHs, we find excellent agreement between our
composite and the $z=2.3$ FIRE MZR of Ma et al. (2016), especially if we apply
the metallicity correction ($\Delta\log(\rm{O/H})\sim 0.06$ dex) to our stack
to account for our bias high in SFR. In comparison to the $z=2$ and $z=3$
IllustrisTNG MZRs of Torrey et al. (2019), we find that our composite lies
$\sim 1.6\sigma$ below the $z=2$ MZR and $\sim 1.3\sigma$ below the $z=3$ MZR,
suggesting that a stronger O/H evolution at fixed $M_{\ast}$ with redshift
(like that seen in Figure 7 of Ma et al. 2016 relative to other simulations of
the time$\,-\,$including the original Illustris simulation Torrey et al. 2014)
and/or tuning of the $z=0$ MZR normalization is needed. Like the direct-method
MZR of Sanders et al. (2020), our composite also disagrees with the existence
of a metallicity bump in the MZR at $M_{\ast}\lesssim 10^{9}\ \rm{M_{\odot}}$,
suggesting that the minimum wind velocity assumption in IllustrisTNG should be
revisited.

###  4.4 $M_{\ast}-\rm{SFR}-\rm{O/H}$ Relation at $z\sim 2.3$

The $M_{\ast}-\rm{SFR}-\rm{O/H}$ relation, or fundamental metallicity relation
(FMR), posits that the MZR has a secondary dependence on SFR such that, when
considering all three properties, the scatter in metallicity at fixed
$M_{\ast}$ and SFR is reduced compared to the scatter in metallicity at fixed
$M_{\ast}$ alone (Mannucci et al., 2010; Lara-López et al., 2010; Andrews &
Martini, 2013; Sanders et al., 2018; Curti et al., 2020; Sanders et al.,
2021). Also as a result of this secondary dependence, at fixed $M_{\ast}$, a
higher (lower) than average SFR yields a lower (higher) than average O/H. Of
further interest is that numerous studies have shown that the FMR is redshift-
invariant to within $\sim 0.1$ dex in metallicity out to $z\sim 3.3$ (e.g.,
Mannucci et al., 2010; Henry et al., 2013a, b; Sanders et al., 2018, 2021). In
effect, the observation of the evolution of the MZR over this redshift range,
whereby O/H decreases with increasing redshift at fixed $M_{\ast}$, is
actually the observation at different redshifts of different portions of the
locally-defined FMR since SFR increases with redshift at fixed $M_{\ast}$
(Whitaker et al., 2014; Speagle et al., 2014; Sanders et al., 2021).
Unfortunately, however, this redshift-invariance of the FMR is still a matter
of debate owing to the uncertainties regarding the applicability of locally-
calibrated, strong-line metallicity calibrations at high-redshift (see Section
4.2), widely-used in lieu of hard-to-measure direct metallicities.

Fortunately, in this work we have a direct, $T_{e}$-based oxygen abundance
from our dwarf galaxy composite spectrum with which we can evaluate the
redshift evolution of the FMR. Additionally, given the properties of our
sample, we are able to probe $M_{\ast}-\rm{SFR}-\rm{O/H}$ parameter space that
to-date has been poorly-sampled at $z\sim 2.3$. We do so via the commonly-used
2D planar projection of the 3D FMR, established by Mannucci et al. (2010). In
this projection, the metallicity is a function of the linear combination of
$M_{\ast}$ and SFR, denoted by $\mu_{\alpha}$ and described by the equation

$\mu_{\alpha}=\log(M_{\ast}/\mathrm{M_{\odot}})-\alpha\log(\mathrm{SFR}/\mathrm{M_{\odot}\,yr^{-1}})$
(5)

where $\alpha$ is the parameter which denotes the strength of the SFR-
dependence of the FMR as well as the value which minimizes the scatter in O/H
at fixed $\mu_{\alpha}$. This parameter, $\alpha$, is generally found to be
lower (a weaker SFR dependence) when determined with strong-line metallicities
(e.g., $\alpha=0.32$; Mannucci et al., 2010) and higher (a stronger SFR
dependence) when determined with direct-method metallicities (e.g.,
$\alpha=0.66$; Andrews & Martini, 2013), though recent estimations via strong-
line metallicities by Curti et al. (2020, $\alpha=0.55$) and Sanders et al.
(2021, $\alpha=0.60$) have brought these $\alpha$-estimates into better
agreement. For this work, we use the value of $\alpha=0.63$, which derives
from the SDSS $M_{\ast}-\rm{SFR}$ stacks of Andrews & Martini (2013), with
their direct-method metallicities corrected for diffuse ionized gas (DIG)
contamination by Sanders et al. (2017).

In Figure 6, we show our $z\sim 2.3$ dwarf galaxy stack (green star; $x$-axis
error bar represents the $\mu_{0.63}$ range of the stacking sample) in the
direct-method $\rm{O/H}-\mu_{0.63}$ parameter space, plotted at the
$\mu_{0.63}$ value ($\mu_{0.63}=8.07$) given by the stacking sample’s median
stellar mass (log($M_{\ast}$/$\rm{M_{\odot}}$) = 8.29) and median SFR
(log(SFR) = 0.353). We compare our composite to the Sanders et al. (2017)
best-fit linear representation of the FMR (black line) tracing the $z\sim 0$,
DIG-corrected, Andrews & Martini (2013) $M_{\ast}-\rm{SFR}$ stacks. We also
include in Figure 6 the two mass bins of the Sanders et al. (2020) $z>1$
auroral-line sample (black squares), which these authors found to be
consistent, within $\sim 0.1$ dex at fixed $\mu_{0.63}$, with lower redshift
samples in the $\mu_{0.63}$ direct-method FMR projection. Adding to these
results, we find our dwarf galaxy composite to have a $T_{e}$-based
metallicity that lies $\sim 0.5\sigma$ or $\sim 0.11$ dex below the linear
relation representing the direct-method FMR. We therefore join numerous other
authors in suggesting, via direct-method metallicities, that the FMR evolves
at most by $\sim 0.1$ dex in O/H at fixed $M_{\ast}$ and SFR from $z=0$ to at
least $z\sim 2.3$, though note that within the uncertainties on the
metallicity, our stack is consistent with a redshift-invariant FMR.

Figure 6: The direct-method $M_{\ast}-\rm{SFR}-\rm{O/H}$ fundamental
metallicity relation (FMR), represented here by the O/H$\,-\,\mu_{\alpha}$
planar projection (Mannucci et al., 2010). Referencing Equation 5 for
$\mu_{\alpha}$, we take $\alpha=0.63$, the value determined by Sanders et al.
(2017) using the DIG-corrected, $T_{e}$-based metallicities of the Andrews &
Martini (2013) SDSS $M_{\ast}-\rm{SFR}$ stacks. The best-fit linear relation
by Sanders et al. (2017) to these $z\sim 0$ stacks is shown here by the black
line. Our $z\sim 2.3$ composite is displayed as the green star with
$\mu_{0.63}$ calculated from our sample’s median $M_{\ast}$
(log($M_{\ast}$/$\rm{M_{\odot}}$) = 8.29) and SFR (log(SFR) = 0.353). The
$x$-axis error bar represents our sample’s range in $\mu_{0.63}$. Our stack
lies $\sim 0.11$ dex or $\sim 0.5\sigma$ below the direct-method FMR, though
is consistent with a redshift-invariant FMR within uncertainties. For
reference, we also show the $M_{\ast}$-binned, $z>1$ auroral-line sample of
Sanders et al. (2020, black squares) and the $z\sim 2.4$ KBSS-LM1 composite of
Steidel et al. (2016, orange diamond).

##  5 Summary

In this study, we analyze the median composite spectrum of 16 typical, star-
forming, dwarf galaxies ($7.06\leqslant\log(M_{\ast}/\rm{M_{\odot}})\leqslant
8.93$; $\log(M_{\ast}/\rm{M_{\odot}})_{\rm{median}}=8.29$) at redshifts
$1.7<z<2.6$ ($z_{\rm{mean}}=2.30$) selected independent of the strength of any
particular emission line. These galaxies are gravitationally-lensed by the
foreground clusters Abell 1689, MACS J0717.5+3745, and MACS J1149.5+2223. In
our composite spectrum, we find a $2.5\sigma$ ($4.1\sigma$) detection of the
faint, $T_{e}$-sensitive, [O III] $\lambda$4363 auroral line when considering
our bootstrapped (statistical-only) error spectrum, allowing us to directly
calculate an oxygen abundance from the composite of
$12+\log(\rm{O/H})_{\rm{direct}}=7.87^{+0.24}_{-0.22}$ ($0.15^{+0.11}_{-0.06}\
\rm{Z_{\odot}}$). We summarize the results using this $T_{e}$-based
metallicity, and other conclusions, in this final section.

  1. 1.

To determine how representative our dwarf galaxy sample is of typical, star-
forming, $z\sim 2.3$ dwarf galaxies, we first considered our composite in the
context of the N2-BPT diagram, where we found that our stack lies offset from
the $z\sim 0$ SDSS star-forming sequence in the same parameter space as the
$z\sim 2.3$ star-forming galaxies of larger statistical surveys (e.g., MOSDEF;
Shapley et al., 2015). We also show that our composite lies at higher [O III]
$\lambda$5007/H$\beta$ and lower [N II] $\lambda$6583/H$\alpha$ than any of
the $M_{\ast}$-binned MOSDEF stacks of Sanders et al. (2021); our composite
extends the trend seen with these MOSDEF stacks of lower stellar mass and
metallicity at higher [O III]/H$\beta$ and lower [N II]/H$\alpha$.

  2. 2.

We also considered our stacking sample against an extrapolation of the $z\sim
2.3$ $M_{\ast}-\rm{SFR}$ SFMS of Sanders et al. (2021), finding the sample to
scatter on either side of this relation. Our stacking sample has a median SFR
= 2.25 $\rm{M_{\odot}\,yr^{-1}}$, which lies $\Delta\log(\rm{SFR})\approx
0.19$ dex above this SFMS at fixed $M_{\ast}$
($\log(M_{\ast}/\rm{M_{\odot}})_{\rm{med}}=8.29$), corresponding to a bias in
O/H of $\Delta\log(\rm{O/H})\approx-0.06$ dex via the FMR, well within even
our statistical metallicity uncertainty ($\sigma\sim 0.12$ dex). We conclude
that our stacking sample is not largely biased in SFR or O/H and is thus, on
average, representative of typical, star-forming, $z\sim 2.3$ dwarf galaxies
with stellar masses between $10^{8}\lesssim M_{\ast}/\rm{M_{\odot}}\lesssim
10^{9}$. Our sample serves as an initial extension of representative,
statistical, spectroscopic surveys at $z\sim 2.3$ into the dwarf galaxy mass
regime.

  3. 3.

We analyzed the applicability at $z\sim 2.3$ of several locally-calibrated,
oxygen-based, strong-line metallicity relations from the literature. We find
that at $12+\log(\rm{O/H})_{\rm{direct}}=7.87^{+0.24}_{-0.22}$, our stack lies
in the metallicity-insensitive “turnover” region of the $\rm{O_{3}}$ and
$\rm{R_{23}}$ calibrations, signalling their ineffectiveness for metallicity
estimation of typical, $z\sim 2.3$, dwarf galaxies. When considering the
$\rm{O_{32}}$ and $\rm{O_{2}}$ indices together, our stack’s metallicity is
most accurately reproduced (within $\lesssim 0.15$ dex) at fixed strong-line
ratio by the local reference calibrations of Bian et al. (2018), in agreement
with that seen for the $z\sim 2.4$ KBSS-LM1 composite of Steidel et al.
(2016), who also selected their sample independent of line-strength. We
generally disagree with the conclusions of Sanders et al. (2020), who argue
that their $z\sim 2.2$ auroral-line sample favors the high-$z$ analog
calibrations of Bian et al. (2018). While both samples are at similar redshift
and metallicity, we argue that our discrepancy in conclusion is due to sample
selection effects as well as biases in the low-metallicity strong-line
calibration samples. Indeed, by being selected for having a detection of a
$T_{e}$-sensitive auroral-line, the sample of Sanders et al. (2020) is
strongly biased in SFR, $\rm{O_{32}}$, and $\rm{EW_{0}}$([O III]
$\lambda$5007) relative to typical, $z\sim 2.3$, star-forming galaxies.

  4. 4.

At the median stellar mass of our stacking sample,
log($M_{\ast}$/$\rm{M_{\odot}}$) = 8.29, we compared our composite direct
metallicity, $12+\log(\rm{O/H})=7.87^{+0.24}_{-0.22}$, against the $z\sim 2.2$
direct-method MZR of Sanders et al. (2020), the $z\sim 2.3$ strong-line MZR of
Sanders et al. (2021), and the $z\sim 2.3$ photoionization-model-based MZR of
Strom et al. (2022). After correcting for the slight SFR bias of the stacking
sample, we find that our $z\sim 2.3$ stack lies $\sim 0.06$ dex above the
direct-method MZR and $\sim 0.07$ dex below the strong-line MZR at fixed
$M_{\ast}$, well within our uncertainties. Our stack lies $\approx 1\sigma$
below the Strom et al. MZR. In constraining the slope of the MZR, we defer to
the direct-method MZR, with which we show excellent agreement, as
metallicities for this relation and our composite were calculated
consistently. Therefore, we suggest that the slope of the $z\sim 2.3$ MZR is
that given by this direct-method relation of Sanders et al. (2020),
$\beta=0.37$.

  5. 5.

We also compared our composite, and the MZRs derived from observations, to the
$z=2.30$ MZR from the FIRE simulations (Ma et al., 2016) as well as to the
$z=2$ and $z=3$ MZRs from the IllustrisTNG100 simulations (Torrey et al.,
2019). At the stack’s fiducial median stellar mass,
log($M_{\ast}$/$\rm{M_{\odot}}$) = 8.29, our composite is consistent within
uncertainties with both sets of simulations. However, when recalculating our
sample stellar masses assuming more realistic non-parametric SFHs, the median
stellar mass is increased to log($M_{\ast}$/$\rm{M_{\odot}}$) = 8.92, moving
the stack into excellent agreement with the FIRE MZR and $\sim 1.5\sigma$
below the $2\leqslant z\leqslant 3$ IllustrisTNG MZR. This tension with
IllustrisTNG is in part caused by a “bump” in its MZR deriving from a constant
minimum wind velocity ($v_{\rm{min}}=350\ \rm{km}\ \rm{s^{-1}}$) applied to
galaxies with $M_{\ast}\lesssim 10^{9}\ \rm{M_{\odot}}$. Between our stack
with recalculated masses and the direct-method MZR of Sanders et al. (2020),
we suggest that the low-mass end of the MZR does not contain this bump.

  6. 6.

Our $z\sim 2.3$ dwarf galaxy composite was compared to the locally-defined,
direct-method FMR in order to test the relation’s redshift-invariance. We made
this comparison via the FMR projection proposed by Mannucci et al. (2010) and
given in Equation 5, with $\alpha=0.63$ (Sanders et al., 2017). At
$\mu_{0.63}=8.07$, calculated with the stacking sample’s median stellar mass,
log($M_{\ast}$/$\rm{M_{\odot}}$) = 8.29, and median SFR, log(SFR) = 0.353, we
find our composite to lie $\sim 0.5\sigma$ or $\sim 0.11$ dex below the
Sanders et al. (2017) best-fit linear relation in direct-method
O/H$\,-\,\mu_{0.63}$ space; this relation is fit to the $z\sim 0$ DIG-
corrected stacks of Andrews & Martini (2013). We therefore agree with many in
the literature who suggest that the FMR is redshift-invariant within $\sim
0.1$ dex at fixed $M_{\ast}$ and SFR from $z\sim 0-2.3$.

This study compliments other larger spectroscopic surveys of representative,
star-forming galaxies at $z\sim 2$ and $M_{\ast}\gtrsim 10^{9}\
\rm{M_{\odot}}$ by serving as an initial extension into the dwarf galaxy mass
regime. In analyzing our sample of representative dwarf galaxies, we are able
to use a direct oxygen abundance to provide initial constraints on the low-
mass slope of the high-$z$ MZR and probe scarcely-studied parameter space of
the FMR. While our sample size is small, our work provides a reference point
for future statistical studies of high-$z$ dwarf galaxies with the newly-
operational James Webb Space Telescope, which will greatly increase our
understanding of the processes responsible for galaxy formation and evolution.

This material is based upon work supported by the National Science Foundation
under Grant No. 1617013. The authors wish to recognize and acknowledge the
very significant cultural role and reverence that the summit of Maunakea has
always had within the indigenous Hawaiian community. We are most fortunate to
have the opportunity to conduct observations from this mountain.

##  Appendix A Consistent Fits of the Strong-Line Metallicity Calibrations of
Bian+18

We note that while Bian et al. (2018) provide fits to their high-$z$ analog
stacks for the strong-line indices of $\rm{R_{23}}$, $\rm{O_{32}}$, and
$\rm{O_{3}}$ (with [O III] $\lambda$4959 added to [O III] $\lambda$5007 in the
latter two indices unlike in this work; see Bian et al. (2018) and Sanders et
al. (2021), Footnote 17), they do not provide fits for the $\rm{O_{2}}$ index
or to the local reference stacks. Therefore, we used orthogonal distance
regression and the stated flux and $T_{e}$-based metallicity values given in
Bian et al. (2018) to consistently fit both samples for each line index
considered in Figure 4 and Section 4.2. As in Bian et al. (2018), third-order
polynomials were assumed for the $\rm{R_{23}}$ and $\rm{O_{3}}$ high-$z$
analog fits; we adopted the same functional form for the $\rm{O_{2}}$ high-$z$
analog fit as well. For the fits of the local reference stacks and these same
strong-line indices, we instead assumed second-order polynomial functions due
to there only being 4 bins with estimated direct metallicities compared to 6
bins for the high-$z$ analog sample. For the $\rm{O_{32}}$ index, as in Bian
et al. (2018), we fit a linear functional form to both the local reference and
high-$z$ analog stacks. The coefficients of our fits are given in Table 4
below.

Table 4: Coefficients of the Refit Bian+18 Relations

Ratio | $c_{0}$ | $c_{1}$ | $c_{2}$ | $c_{3}$
---|---|---|---|---
Local Reference Relations
$\rm{O_{32}}$ | 11.1767 | -1.3414 |  |
$\rm{O_{3}}$ | -57.7632 | 14.9708 | -0.9578 |
$\rm{O_{2}}$ | -25.2699 | 5.3889 | -0.2750 |
$\rm{R_{23}}$ | -26.4726 | 6.9012 | -0.4345 |
High-$z$ Analog Relations
$\rm{O_{32}}$ | 14.5895 | -1.7287 |  |
$\rm{O_{3}}$ | 117.8668 | -48.8336 | 6.7607 | -0.3107
$\rm{O_{2}}$ | -80.9802 | 14.1526 | 0.0755 | -0.0723
$\rm{R_{23}}$ | 263.2135 | -101.2448 | 13.0313 | -0.5591

Note. — Ratio definitions given in Section 4.2.

## References