# CECILIA: The Faint Emission Line Spectrum of $z\sim 2-3$ Star-forming
Galaxies

Allison L. Strom Department of Physics and Astronomy, Northwestern University,
2145 Sheridan Road, Evanston, IL 60208, USA  Center for Interdisciplinary
Exploration and Research in Astrophysics (CIERA), Northwestern University,
1800 Sherman Avenue, Evanston, IL 60201, USA  Gwen C. Rudie Carnegie
Observatories, 813 Santa Barbara Street, Pasadena, CA 91101, USA  Ryan F.
Trainor Department of Physics and Astronomy, Franklin & Marshall College, 637
College Avenue, Lancaster, PA 17603, USA  Gabriel B. Brammer Cosmic Dawn
Center (DAWN), Niels Bohr Institute, University of Copenhagen, Jagtvej 128,
København N, DK-2200, Denmark  Michael V. Maseda Department of Astronomy,
University of Wisconsin, 475 N. Charter Street, Madison, WI 53706, USA
Menelaos Raptis Department of Physics and Astronomy, Franklin & Marshall
College, 637 College Avenue, Lancaster, PA 17603, USA  Noah S. J. Rogers
Department of Physics and Astronomy, Northwestern University, 2145 Sheridan
Road, Evanston, IL 60208, USA  Center for Interdisciplinary Exploration and
Research in Astrophysics (CIERA), Northwestern University, 1800 Sherman
Avenue, Evanston, IL 60201, USA  Minnesota Institute for Astrophysics,
University of Minnesota, 116 Church St. SE, Minneapolis, MN, 55455  Charles C.
Steidel Cahill Center for Astronomy and Astrophysics, California Institute of
Technology, MS 249-17, Pasadena, CA 91125, USA  Yuguang Chen (陈昱光) Cahill
Center for Astronomy and Astrophysics, California Institute of Technology,
1200 E California Blvd, MC249-17, Pasadena, CA 91125, USA  David R. Law Space
Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
Allison Strom<EMAIL_ADDRESS>
###### Abstract

We present the first results from CECILIA, a Cycle 1 JWST NIRSpec/MSA program
that uses ultra-deep $\sim 30$ hour G235M/F170LP observations to target
multiple electron temperature-sensitive auroral lines in the spectra of 33
galaxies at $z\sim 1-3$. Using a subset of 23 galaxies, we construct two $\sim
600$ object-hour composite spectra, both with and without the stellar
continuum, and use these to investigate the characteristic rest-optical
($\lambda_{\rm rest}\approx 5700-8500$ Å) spectrum of star-forming galaxies at
the peak epoch of cosmic star formation. Emission lines of eight different
elements (H, He, N, O, Si, S, Ar, and Ni) are detected, with most of these
features observed to be $\lesssim 3\%$ the strength of H$\alpha$. We report
the characteristic strength of three auroral features ([N II]$\lambda 5756$,
[S III]$\lambda 6313$, and [O II]$\lambda\lambda 7322,7332$), as well as other
semi-strong and faint emission lines, including forbidden [Ni
II]$\lambda\lambda 7380,7414$ and permitted O I $\lambda 8449$, some of which
have never before been observed outside of the local universe. Using these
measurements, we find $T_{e}$[N II] $=13630\pm 2540$ K, representing the first
measurement of electron temperature using [N II] in the high-redshift
universe. We also see evidence for broad line emission with a FWHM of
$536^{+45}_{-167}$ km s-1; the broad component of H$\alpha$ is $6.01-28.31$%
the strength of the narrow component and likely arises from star-formation
driven outflows. Finally, we briefly comment on the feasibility of obtaining
large samples of faint emission lines using JWST in the future.

††facilities: JWST (NIRSpec)††software: BPASSv2 (Stanway et al., 2016;
Eldridge et al., 2017), Cloudy (Ferland et al., 2013), GalDNA (Strom et al.,
2018), JWST Calibration Pipeline (Bushouse et al., 2023), grizli (Brammer,
2023), msaexp (Brammer, 2022), PyNeb (Luridiana et al., 2015)

##  1 Introduction

The nebular emission lines originating in galaxies’ star-forming regions are
among the most powerful tools available for investigating the physical
conditions in galaxies at all redshifts (see Kewley et al. 2019 for a review).
Recombination lines of hydrogen provide constraints on star formation rates
(SFRs) and reddening due to dust; weaker recombination lines of heavier
elements such as oxygen are direct tracers of gas-phase enrichment; and
collisionally-excited forbidden transitions of elements like oxygen, nitrogen,
carbon, and sulfur are variously sensitive to electron temperature ($T_{e}$)
and electron density ($n_{e}$), as well as the ionization state and enrichment
of the photoionized gas.

Some of these lines, including H$\alpha$ at 6564 Å and the [O
III]$\lambda\lambda 4960,5008$ doublet, are bright enough to be detected out
to large cosmological distances, even with relatively short exposure times.
Indeed, over the last decade, ground-based near-infrared (NIR) spectrographs
such as Magellan/FIRE (Simcoe et al., 2008, 2010), Keck/MOSFIRE (McLean et
al., 2010, 2012), and VLT/KMOS (Sharples et al., 2013) have led to samples of
1000s of $z\sim 2-3$ galaxies with measurements of such emission lines (e.g.,
Masters et al., 2014; Steidel et al., 2014; Shapley et al., 2015; Sanders et
al., 2015; Wisnioski et al., 2015; Strom et al., 2017), allowing their
metallicity, ionization and excitation properties, and gas density to be
studied in comparable detail to large samples of $z\sim 0$ galaxies (e.g.,
Kauffmann et al., 2003; Brinchmann et al., 2004; Tremonti et al., 2004;
Belfiore et al., 2015; Mingozzi et al., 2020). In the last year, JWST/NIRSpec
and JWST/NIRCam grism observations have extended these efforts to even higher
redshifts ($z\gtrsim 3-6$) by enabling IR spectroscopy out to longer
wavelengths (e.g., Kashino et al., 2023; Kocevski et al., 2023; Oesch et al.,
2023; Shapley et al., 2023; Sun et al., 2023).

Other lines—including metal recombination lines and the $T_{e}$-sensitive
auroral lines of heavy elements, which are both key probes of chemical
enrichment—are faint enough that they are not routinely detected even in
spectra of nearby galaxies. Despite significant investment of observing time
on some of the largest ground-based telescopes in the world, measurements of
auroral [O III]$\lambda 4364$ were only possible for a handful of individual
galaxies at $z\gtrsim 2$ prior to the launch of JWST (Christensen et al.,
2012; James et al., 2014; Sanders et al., 2020). Spectroscopic observations
with JWST promise to yield unprecedented numbers of auroral emission line
measurements in high-$z$ galaxies. The first analyses of the early release
observations (ERO; Pontoppidan et al., 2022) in the SMACS J0723.3–7327 field
reinforced this expectation, with significant detections of [O III]$\lambda
4364$ in several $z\sim 8$ galaxies (Arellano-Córdova et al., 2022; Schaerer
et al., 2022; Taylor et al., 2022; Brinchmann, 2023; Curti et al., 2023; Katz
et al., 2023; Rhoads et al., 2023; Trump et al., 2023; Trussler et al., 2023).
At the same time, many of these studies reported conflicting gas-phase oxygen
abundance (O/H) measurements in the same objects, and it was unclear how
representative this early, very high-$z$ sample might be. Subsequent work has
revisited the issue of auroral line detections in JWST observations of tens of
high-$z$ galaxies (albeit primarily at low to moderate O/H), confirming
suspicions that locally-calibrated metallicity diagnostics are likely
unsuitable for the majority of high-$z$ galaxies (Laseter et al., 2023;
Sanders et al., 2023a). To date, however, consensus regarding how best to
leverage these measurements to, e.g., understand the overall distribution of
chemical enrichment in the early universe has not yet been achieved.

In spite of these challenges, the community has collectively recognized the
goal of using auroral line measurements and the resulting direct-method
metallicities to construct more accurate methods of measuring high-$z$ galaxy
enrichment _in situ_. This is evidenced by the selection of three separate
Cycle 1 JWST programs (PIDs 1879, 1914, and 2593) by the time allocation
committee, with a total investment of over 150 hrs, or $\sim 2.5$% of all the
GO time available in Cycle 1. Here, we report the first results from PID 2593,
also known as CECILIA (Chemical Evolution Constrained Using Ionized Lines in
Interstellar Aurorae; Strom et al. 2021, data accessible via doi:
https://doi.org/10.17909/x66z-p144 (catalog 10.17909/x66z-p144)).

CECILIA was designed to measure auroral [S III]$\lambda 6313$ and [O
II]$\lambda\lambda 7322,7332$ in the spectra of a carefully selected sample of
$z\sim 2-3$ star-forming galaxies, using $\sim 30$ hr G235M/F170LP
observations. Owing to the unique depth of these data, CECILIA is also able to
detect myriad other lines in the galaxies’ rest-optical spectra, some of which
are stronger than any auroral emission line and, thus, more likely to be
observed in more typical integration times with JWST. Consequently, it is
important to understand the expected strength of these faint and semi-strong
emission lines, in order to guide future studies using JWST, as well as with
other current and future facilities.

The remainder of this letter focuses on two $\sim 600$ object-hour rest-
optical composite spectra of $z\sim 2-3$ galaxies observed as part of the
CECILIA survey, with the aim of providing an “atlas” of the characteristic
faint emission line spectrum of high-$z$ galaxies. We describe the CECILIA
survey—including the galaxy sample, the JWST program, and the data
reduction—in Section 2. Section 3 outlines the construction of the composite
spectra and their key features, with a more in-depth discussion of individual
emission lines in Section 4. In Section 5, we close with a summary of our
findings and a brief discussion of implications for future observations of
faint emission lines in $z\gtrsim 2$ galaxies. Throughout the text, we refer
to specific spectral features using their vacuum wavelengths.

##  2 The CECILIA Survey

Figure 1: The HST/WFC3 F140W image of Q2343+125 field targeted by CECILIA is
shown in greyscale, with the approximate 3.6′$\times$3.4′ NIRSpec field of
view overlaid in black. The observed slits are magnified by $\sim$3$\times$
for visibility, with red slits and names indicating the 23 sources included in
this work.

The principal goal of CECILIA is to measure multiple faint rest-optical
auroral lines in the spectra of $z\sim 2-3$ galaxies, which can then be used
to calibrate new high-$z$ metallicity diagnostics. Some of the galaxies
observed as part of CECILIA have preexisting rest-optical spectra obtained
using Keck/MOSFIRE (Steidel et al., 2014; Trainor et al., 2015; Strom et al.,
2017), but even the strongest auroral lines are not routinely detected for
individual galaxies in deep ($\sim 8-10$ hr) observations. Although
JWST/NIRSpec provides greater sensitivity and spectral coverage than ground-
based NIR spectrographs, achieving this goal still pushes the limits of the
observatory. To make the best use of JWST observing time, we first used
detailed photoionization models and existing ground-based rest-ultraviolet
(UV) and rest-optical spectra of the same galaxies to robustly predict the
auroral line strengths. We then used these predictions together with the (pre-
flight) JWST Exposure Time Calculator (ETC)111https://jwst.etc.stsci.edu/ to
identify the depth needed to detect the auroral lines in individual galaxies.
Below, we describe the parent galaxy sample, the emission line predictions,
the design of the NIRSpec program, including exposure time requirements and
microshutter assembly (MSA) design, and the reduction of the JWST data.

###  2.1 Parent Galaxy Sample and Field Selection

CECILIA targets galaxies drawn from the Keck Baryonic Structure Survey (KBSS;
Steidel et al., 2010; Rudie et al., 2012; Trainor et al., 2015; Strom et al.,
2017). KBSS is a large spectroscopic survey of UV-color and narrow-band
Ly$\alpha$ selected galaxies in 15 fields. The survey includes deep $J$, $H$,
and $K$ NIR (rest-optical) spectroscopy from Keck/MOSFIRE, deep optical (rest-
UV) spectroscopy from Keck/LRIS, and imaging in $U_{n}$ through $K_{s}$,
HST/WFC3 F140W, F160W, Spitzer IRAC Ch$1-4$ and MIPS 24 $\mu$m, and narrow-
band (rest-frame) Ly$\alpha$ filters. In total, KBSS comprises $\sim 3500$
$z\sim 1-3$ galaxies with confirmed (rest-UV or rest-optical) spectroscopic
redshifts and rich multiwavelength data. Thanks to the dense sampling of
galaxies with known redshifts and high-quality rest-optical spectroscopy in
KBSS, targeting one of the survey fields allows us to efficiently prioritize
those galaxies predicted to yield auroral line detections, where we can also
ensure that all of the emission lines required to determine $T_{e}$ and
direct-method metallicity are accessible.

The primary targets for CECILIA are UV-color-selected galaxies for which
detecting the auroral [S III]$\lambda 6313$ line is feasible, as this is the
faintest line required to achieve the primary goals of the program (see
Section 2.2). We focused on galaxies with $2.10\leq z\leq 2.68$, where all of
the required emission lines fall within the G235M, G395M, and ground-based NIR
spectral bandpasses. Galaxies were prioritized if they have: somewhat higher
redshifts ($z>2.3$), where NIRSpec is more sensitive at the observed
wavelength of [S III]$\lambda 6313$; smaller than average sizes, increasing
the auroral surface brightness; high SED-based SFRs ($>24~{}\rm
M_{\odot}~{}\rm yr^{-1}$); and/or large observed nebular [O III]$\lambda 5008$
or H$\alpha$ line fluxes ($>7.0\times 10^{-17}$ erg s-1 cm-2). All of these
properties increase the ease of detection with the NIRSpec/MSA. The highest
priority targets were galaxies with detailed emission line models (Section
2.2) whose predicted auroral line surface brightnesses exceeded the detection
threshold of the planned observations; galaxies with models predicting non-
detections were down weighted. Narrow-band selected Ly$\alpha$ emitters (LAEs)
from Trainor et al. (2016) with spectroscopic detections of Ly$\alpha$ and [O
III]$\lambda 5008$ or H$\alpha$ were also prioritized as a way of extending
the galaxy sample to lower stellar masses (M∗) and SFRs.

Of the 15 KBSS fields, we selected the Q2343+125 field due to its high density
of high-priority sources and large catalog of LAEs at $z\approx 2.55$ with
spectroscopic redshifts. Further, this field also has an existing HST/WFC3
F140W mosaic (Figure 1) that provided both the precision astrometry required
for mask design and the galaxy size measurements needed for target
prioritization—without requiring additional (pre-)imaging from space using
JWST or HST.

The CECILIA JWST/NIRSpec observations contain a total sample of 34
galaxies.222One target, Q2343-D27, appears to be a $z=0.0890$ interloper,
based on the JWST/NIRSpec observations. We include 23 of these objects here
(Figure 1), omitting the Ly$\alpha$-selected galaxies that do not have secure
SED models (4 galaxies), galaxies at $z<2$ (4 galaxies), and sources that were
severely impacted by shutter failures in the NIRSpec/MSA (3 galaxies). The
final sample is largely typical of KBSS galaxies, with $\langle z\rangle=2.4$,
masses spanning log(M∗/M${}_{\odot})=8.5-10.7$ and a median value of
log(M∗/M${}_{\odot})=9.7$ (assuming a Chabrier 2003 stellar initial mass
function). Based on H$\alpha$ and H$\beta$ measurements from ancillary MOSFIRE
spectra, the included galaxies have SFRs ranging from $16-42$ M⊙ yr-1, with a
median SFR${}_{\rm H\alpha}=21$ M⊙ yr-1. These are slightly lower than median
values reported in Strom et al. (2017), which were determined in the same
manner, but similar in terms of M∗ to the subsample of KBSS galaxies used to
construct the deep “LM1” composite in Steidel et al. (2016).

Figure 2: The hatched regions show the model predictions for auroral lines in
the rest-optical spectra of typical $z\sim 2-3$ galaxies, separated on the
basis of whether they fall in the G140M bandpass (top panel) or the G235M
bandpass (bottom panel). The width of the hatched regions reflects the typical
range of ionization parameter $U$ in high-$z$ galaxies. The predicted line
fluxes for [O III]$\lambda 4364$ (blue hatched region) are $\sim 1$ dex
fainter than the depth of typical ground-based spectra of individual galaxies,
represented by the distribution of $3\sigma$ upper limits on [O III]$\lambda
4364$ from KBSS (dashed black histogram). Sanders et al. (2020) reported four
ground-based detections of unlensed [O III]$\lambda 4364$ (blue points,
shifted up by 0.24 dex to match the photoionization model abundance scale),
but these galaxies appear atypical. Estimates using the pre-flight ETC
indicated that detecting [O III]$\lambda 4364$ in a representative sample of
high-$z$ galaxies would be prohibitively expensive; the $3\sigma$ limiting
line flux of $6\times 10^{-19}$ erg s-1 cm-2 achievable in a combined 29 hr
G140M exposure (red line) probes $\lesssim 30\%$ of the $z\sim 2-3$ sample
from Strom et al. (2018). Fortunately, the typical predicted line fluxes for
the sum of the [O II]$\lambda\lambda 7322,7332$ lines (purple hatched region)
and [S III]$\lambda 6313$ (orange hatched region) could be detected for
galaxies with a wider range of $U$ and O/H in the same exposure time, due to
the higher sensitivity of NIRSpec in G235M. A $3\sigma$ limiting line flux of
$4.1\times 10^{-19}$ erg s-1 cm-2 reaches $\sim 90\%$ of typical galaxies.

###  2.2 Emission Line Predictions

The expected strengths of the auroral emission lines targeted by CECILIA were
determined using photoionization models designed to reconcile the rest-UV and
rest-optical spectra of $z\sim 2-3$ galaxies. We used a combination of the
Binary Population and Spectral Synthesis models (BPASSv2; Stanway et al.,
2016; Eldridge et al., 2017) and Cloudy photoionization models (Cloudy13;
Ferland et al., 2013) to predict line strengths as a function of gas-phase
metallicity (O/H). We matched the model parameters (gas density $n_{\rm H}$,
stellar Fe/H, and ionization parameter $U$) to the properties of $z\sim 2-3$
KBSS galaxies reported by Strom et al. (2018), which are consistent with the
values reported for other $z\sim 2-3$ samples (e.g., Topping et al., 2020).
The model outputs were then converted to line fluxes using a representative
range of SFRs and dust extinction.

Figure 2 presents the predictions for three of the brightest rest-optical
auroral emission lines as a function of O/H, with the width of the hatched
regions corresponding to the typical range of $U$ in high-$z$ galaxies; lower
ionization galaxies have fainter lines at fixed metallicity. Re-calibrating
strong line metallicity diagnostics and photoionization models at $z\gtrsim 2$
requires measuring auroral lines in galaxies spanning both O/H and $U$, as
both directly influence the strength of nebular emission lines. The top panel
in Figure 2 shows the steep decline in [O III]$\lambda 4364$ with increasing
O/H and implies a limited ability to detect typical $z\sim 2$ galaxies at high
O/H and/or low $U$ using JWST/NIRSpec, even with long exposures in G140M.

In contrast, the bottom panel of Figure 2 shows that a $3\sigma$ line flux
sensitivity of 4.1$\times 10^{-19}$ erg s-1 cm-2 in G235M (corresponding to a
total $\approx 30$ hr exposure time using the pre-flight ETC; see Section
2.3.1) enables the detection of [S III]$\lambda 6313$ _and_ [O
II]$\lambda\lambda 7322,7332$ at virtually all $U$, even in galaxies with
relatively high gas-phase O/H. It is comparatively easier to detect [S
III]$\lambda 6313$ and [O II]$\lambda\lambda 7322,7332$ not only because they
are predicted to be intrinsically brighter than [O III]$\lambda 4364$ in the
same galaxies, but also because of the increasing sensitivity of JWST/NIRSpec
at longer wavelengths. On the basis of these predictions, we elected to obtain
deep spectra of galaxies in a _single_ configuration, in order to maximize the
overall number of auroral lines detected for individual galaxies with a range
of O/H and $U$.

###  2.3 JWST/NIRSpec Program Design

To optimize the efficiency of the JWST program, we generated a large grid of
ETC simulations spanning a range of galaxy sizes, limiting line fluxes, MSA
centering constraints, and redshifts, as well as a comparable grid of MSA
Planning Tool (MPT) simulations that considered the full range of available
centering constraints. In this section, we describe the most salient elements
of the program design.

####  2.3.1 Exposure time requirements

NIRSpec G235M observations of [S III]$\lambda 6313$ and [O II]$\lambda\lambda
7322,7332$ in CECILIA galaxies were modeled using an exponential surface
brightness profile (Sersic index $n=1$) with a projected semi-major axis of
0$\farcs$26 and an axis ratio of $b/a=0.6$, consistent with the measured
morphologies and median sizes of galaxies in our parent sample (Law et al.,
2012).

Pre-flight ETC simulations showed that reaching the required $3\sigma$
limiting line flux of $4.1\times 10^{-19}$ erg s-1 cm-2 for a median-sized
galaxy at $z=2.3$ at the edge of the midpoint tolerance (see Section 2.3.2)
required 29.5 hours of exposure time (20 groups $\times$ 6 integrations
$\times$ 12 exposures) using NRS IRS2 readouts. Our exposure time calculations
assumed “MSA Full Shutter Extraction” and assumed we would need pixel-level
subtraction from A-B pairs. As we discuss below in Section 2.4.1, we have
instead implemented a global background model drawn from slits across the full
MSA, which reduces the overall noise in the final combined data compared to
the conservative assumptions in our original calculations. For the majority of
the sources in our catalog, ETC calculations demonstrated that some of the
background region in each spectrum would be contaminated with light from the
source, and the derived exposure time requirements took this effect into
account.

####  2.3.2 MSA design

The MSA configuration is central to the success of CECILIA, and considerable
experience with ground-based multi-object mask design led us to conduct
extensive trials using different mask parameters in the MSA Planning Tool
(MPT). We experimented with all possible centering constraints, dithering and
nodding options, and three- and five-slitlet length slits. We ran trial masks
spanning the full range of allowable PAs, using small steps in both position
and PA to understand the sensitivity of the optimal configuration to changes
in PA. Based on more than 100 runs of the MPT considering more than 70 million
unique configurations, we determined that many PAs have $<$ 60% as many high-
priority targets as the best masks.

We optimized the MSA centering constraint, which trades exposure time against
sample size, by considering a grid of ETC and MPT runs. Our ETC calculations
spanned the full redshift range of source galaxies and sizes ranging between
the 1st and 3rd quartile of the KBSS size distribution. We considered the S/N
penalty for galaxies at maximal offset in the dispersion direction for each of
the three possible centering restrictions.333The assumed offsets in the ETC
are 0$\farcs$063 for “constrained,” 0$\farcs$076 for “midpoint,” and
0$\farcs$099 for “entire open shutter.” These values use the specified margin
and the pitch of the MSA shutters, thus accounting for the current uncertainty
in the width of the MSA bars. For galaxies with the median size in our sample,
the S/N penalties compared to a perfectly centered target are $7-13$% for
“constrained,” $11-19$% for “midpoint,” and $14-26$% for “entire open
shutter,” where the reported ranges represent different relative angles
between the short axis of the slit and the major axis of the galaxies. MPT
runs showed that “constrained” configurations allowed for only 60% of the
high-priority targets to be placed on a mask compared to the “midpoint”
criteria. Relaxing the centering further via “entire open shutter” constraint
only increased the number of high-priority targets by 7%. Therefore, we
selected the “midpoint” centering constraint for CECILIA observations.

We designed custom software that processed the MPT MSA configurations to check
the wavelength coverage444The post-flight version of MPT now has the ability
to output the wavelength coverage of individual slitlets. (using
MSAViz555https://github.com/spacetelescope/msaviz) and confirmed that primary
targets assigned to a slit on the MSA would have spectral coverage of the
required auroral and nebular lines. This software also considers the known
emission line properties, M∗, and SFRs of target galaxies, which we used to
select a final mask configuration that appropriately sampled the parent sample
to enable an effective metallicity calibration. We selected a default three-
shutter slitlet shape with a three-point nod pattern within the slitlet.

Upon scheduling, we were assigned a Aperture Position Angle, APA = 20.0, with
values from $18.5<$ APA $<20.0$ able to be accommodated within the scheduling
window. At this point, we completed a second set of MPT simulations, including
PA steps of 0.1 degree and $0\farcs 025-0\farcs 01$ position steps to optimize
the PA and final mask. We did not reach convergence,666We define a converged
mask as one where the same optimal mask is returned even when the step size is
decreased. even with angle and position steps much finer than suggested by
JDox, suggesting that significant computational resources would be required to
fully optimize NIRSpec/MSA observations. Based on our simulations, we
ultimately selected an APA = 19.3. Over the 1.5 degree range allowable within
the plan window, the MPT resulted in more than a 30% variation in the number
of high-priority targets, and we advocate for conducting similar PA
optimization to maximize the efficiency of other NIRSpec/MSA programs with low
to moderate density of high priority targets.

Following the selection of pointing and PA, we ran MPT with an expanded
catalog to (1) check for contamination in any of the shutters known to be
stuck open as of June 2022 and (2) open shutters on dark regions of the sky to
sample the background light across the field. These sky slitlets are described
in our modeling of the global sky background in Section 2.4.1. Once the
automated MSA configuration was determined by MPT, the solution was hand-
edited using the MSA configuration editor to (1) elongate slits for high
priority targets where possible, (2) add more background shutters close to
high priority targets to better sample relevant wavelength or field-position
changes in the background, and (3) add high priority targets that did not meet
our centering constraints but could be placed on a mask without conflicting
with other high priority targets. The final MSA design included 34 sources, 23
of which are included in the stacked spectra presented in this letter.

###  2.4 JWST/NIRSpec Data Reduction

The uncalibrated raw G235M data (uncal) frames were processed using the
jwst_level1 pipeline in the grizli777https://github.com/gbrammer/grizli
package version 1.8.9 from Brammer (2023). The level 1 pipeline in grizli uses
the calwebb Science Calibration Pipeline (Bushouse et al. 2023, version
1.10.0, CRDS_CONTEXT $=$ jwst_1100.pmap) for the group_scale correction,
initial flagging of bad pixels and saturated pixels, bias subtraction
(including corrections to the bias using reference pixels), as well as
corrections for detector linearity and persistence and subtraction of the dark
current. Following these calwebb steps, clusters of pixels affected by
snowball cosmic ray events were flagged and the ramp fit was calculated,
including additional processing to detect and remove the effects of cosmic
rays and detector defects. Finally, the gain correction was applied, resulting
in the level 1 processed rate files.

Next, the level 1 processed files were corrected for correlated read noise,
which manifests as vertical banding in the rate files. This $1/f$ noise,
driven by small temperature variations in the ASIC readout electronics, was
modeled and removed using the NSClean algorithm from Rauscher (2023). NSClean
requires the user to create a mask that identifies areas on each of the two
NIRSpec detectors that are unilluminated by source light; these areas are thus
relatively clean tracers of the correlated readout noise. We tested many
different mask design strategies in order to remove as much of the large-scale
vertical banding as possible while also limiting the introduction of
additional high-frequency noise, which we found to be a side effect of the
NSClean algorithm in many cases. We determined the most effective masks for
our program omitted entire rows of pixels in the rectified full-detector image
if any portion of that row was illuminated. Mask designs that omitted only the
limited range of pixels that were illuminated in a given row resulted in
higher levels of high-frequency noise being introduced in the regions of the
detectors that were illuminated by source light.

Following the $1/f$ noise correction by NSClean, we applied the preprocessing
routine steps from msaexp888https://github.com/gbrammer/msaexp version 0.6.11
from Brammer (2022), aside from the $1/f$ noise correction. This routine
repeated the search for snowballs and additional detector defects, which were
also masked. We applied a bias offset correction calculated from the median of
unilluminated pixels in each frame and rescaled the read noise array
associated with each exposure so that it reflected the distribution of the
same unilluminated pixels.

Next, we used msaexp to call the calwebb_spec2 JWST Science Calibration
Pipeline (Bushouse et al. 2023, version 1.10.0), which computed the world
coordinate system reference frame for the data (including the wavelength
calibration), extracted the individual 2D spectra for each slit, and flat-
fielded each 2D spectral cutout. Each spectral cutout was corrected for path
loss assuming the sources uniformly illuminate the slit (i.e., using the
PATHLOSS_UN correction). Note that the current calwebb_spec2 pipeline does not
apply path loss corrections for slits more than three shutters in length, of
which there are many in CECILIA, for reasons described in Section 2.3.2).
Thus, we modified the pipeline to apply the uniform source path loss
correction to all slits. The pipeline correction for the bar shadows produced
by the discretized MSA slitlets was then applied to the data, although, as
described in Section 2.4.1, the pipeline correction left residual bar shadows
on the data and background illumination. The calwebb_spec2 photom step then
provided a final correction to the photometric calibration of the data,
resulting in flux-calibrated 2D spectra for each slit and exposure. Finally,
we used the msaexp drizzle routine to resample the individual 2D spectra onto
a common rectified pixel grid and combine the exposures for each slit with
outlier rejection, using a threshold of 100.

####  2.4.1 Background subtraction and extraction

To correct the data for background light, we opted to use a full-MSA
background solution, rather than a paired exposure differencing algorithm, for
several reasons. First, subtracting a global background model maximizes the
S/N in the final spectra by excluding the shot noise that would be added by
using a low-S/N measure of the background from single adjoining shutters.
Second, the CECILIA targets are extended objects with light from each galaxy
contaminating the shutters above and below the primary shutter. As such, the
typical background algorithms that directly subtract the detected signal above
and below the primary shutter inevitably subtract some source light as well.
This over-subtraction poses a particular issue at the wavelengths of bright
emission lines, which frequently extend well beyond the typical 0$\farcs$6
dither spacing of our observations in the background-subtracted 2D spectra.
Finally, as described below we found that a single global background model
provided a good description of the background across the field, while also
enabling useful checks on the systematics of our observations.

We constructed the global background model by combining data from all the
illuminated shutters in the MSA. Each rectified and drizzled 2D science
spectrum was masked to omit rows corresponding to continuum emission from
target galaxies or from other sources identified in the slit. Pixels
illuminated by extended emission lines were also masked. The full set of
masked science spectra (including those from dedicated sky slitlets, which
were not masked unless they included coincidental sources) were then median
combined into a single 2D background model. The 2D background was averaged in
the spectral direction in order to model the residual bar shadows that were
not fully corrected by the calwebb_spec2 pipeline, and these residual bar
shadows were then removed from the 2D background model. We then averaged the
resulting 2D background model in the spatial direction, weighting each pixel
by the number of spectra contributing to the 2D median at the corresponding
point in order to construct a 1D average background model as a function of
wavelength.

As a cross-check on the consistency of our global background model, we created
similar models from subsets of the observed slits grouped by their position on
the sky (quartiles in right ascension and declination), on the MSA (quadrants
1, 2, 3, and 4), as well as by separating the portions of spectra falling on
each detector (NRS1 and NRS2). The estimated 1D background was consistent
across the field, but we found a small additive offset999The source of this
offset is not clear, but we speculate that it stems from inconsistencies in
the bias estimated from each detector. between the NRS1 and NRS2 detectors. We
therefore applied a compensatory offset to the portions of each slit’s
recorded background spectrum falling on NRS2 before averaging the data from
all slits to create the global background model.

Before subtracting the background from each 2D science spectrum, a slit-
specific version of the global background mode was created that accounts for
the NRS2 offset on the portions of that spectrum that fall on NRS2. This 1D
model is then subtracted from each 2D science spectrum. Notably, the
background model we derived is generally consistent with the predictions of
the JWST Background Tool101010https://jwst-docs.stsci.edu/jwst-other-
tools/jwst-backgrounds-tool (JBT) for our observations. However, a small
additive offset is required to make JBT prediction match the normalization of
our empirical background model, and there are a number of small-scale spectra
features in the empirical background that are unresolved or not included in
the JBT spectrum.

Optimal extraction of the 1D spectra was performed using routines from msaexp.
A spatial profile of the continuum emission for each background-subtracted 2D
spectrum was created by averaging along the wavelength dimension after
weighting by the pipeline-produced 2D weight mask and applying a sigma-
clipping algorithm to mask contaminated pixels and bright emission lines. An
analogous spatial profile of the nebular emission was also created for each
source by averaging the 2D spectrum over small wavelength ranges centered at
the locations of bright emission lines. Each resulting 1D spatial profile was
then fit independently with a Gaussian model. The resulting fits were
typically similar for the continuum and emission line profiles, with the
median profile being 20% wider for the emission lines than the continuum. We
used the Gaussian emission-line spatial model to provide the weights for the
optimal extraction, except in one case where the continuum profile was used
owing to a visibly-poor fit to the emission lines.

Despite the efforts described above, there are still unresolved issues in the
data reduction resulting from known issues with JWST data
products,111111https://jwst-docs.stsci.edu/jwst-calibration-pipeline-
caveats/known-issues-with-jwst-data-products including uncertain variations in
the spectral response as a function of slit position. Likewise, the
unexplained additive offset between the NRS1 and NRS2 detectors, the residual
bar shadows in the pipeline-processed 2D spectra, and the disagreement between
our estimated background and the JBT predictions suggest that there are
systematic effects (perhaps related to detector bias) that are incorrectly
handled by the current pipeline tools and have uncertain downstream effects.
While these uncertainties are not tolerable for the primary goal of
CECILIA—precise abundance determinations of individual galaxies—we expect the
stacking and normalization procedures described in Section 3 likely mitigate
any systematic effects on our composite spectra.

##  3 The Characteristic Rest-optical Spectrum of $\langle z\rangle\sim 2.4$
Galaxies

Figure 3: The median-combined composite spectrum for the CECILIA sample is
shown by the medium blue line, where the individual galaxy spectra are scaled
by the observed continuum at $6800-7000$ Å before being combined. The stack is
then re-scaled so that the continuum in the same wavelength interval matches
the median rest-frame continuum for the constituent galaxies. The 68% CI for
the composite is indicated by the light blue shading, and detected emission
lines of eight different elements (H, He, N, O, Si, S, Ar, and Ni) are
identified by dotted grey lines. The inset panel shows a zoomed-out version of
the composite, centered on H$\alpha$, [N II]$\lambda\lambda 6550,6585$, and [S
II]$\lambda\lambda 6718,6733$, which are the only lines routinely observed in
ground-based observations and shallower JWST spectra of individual high-$z$
galaxies; the grey shaded region indicates the flux range shown in the full
figure, where many fainter lines are visible.

CECILIA contains some of the deepest spectra obtained during Cycle 1, with
$\sim 30$ hr observations of individual galaxies using the NIRSpec/MSA and the
G235M/F170LP configuration. These data offer a unique opportunity to
investigate the spectra of high-$z$ star-forming galaxies, revealing features
that have long remained out of reach of ground-based observations. Given the
uncertainties in the data reduction at the present time, we use composite
spectra as a tool to investigate the nebular emission lines observed in our
data. We have two principal aims: (1) to illustrate the archetypal red rest-
optical ($\lambda_{\rm rest}\approx 5700-8500$ Å) spectrum of a $z\sim 2$
galaxy and (2) determine the typical range of emission line strengths. To
achieve these goals, we construct two composite spectra, one including the
stellar continuum and one only including the nebular emission. In this
section, we describe how the two composite spectra are created, as well as
their key features.

###  3.1 The Total Composite Spectrum

The flux scale of each reduced 1D spectrum is adjusted by comparing the
observed continuum with the best-fit spectral energy distribution (SED) model
of the same galaxy. This strategy has become common practice in analyzing JWST
spectra of high-$z$ galaxies as a way of accounting for uncertainties in the
flux calibration. Specifically, we mask regions of the spectra with large
deviations from the median flux level ($\geq 2\times$ the median absolute
deviation), which excludes not only strong emission lines but also any serious
artifacts remaining in the data due to bad pixels and cosmic rays. We then use
a low-order polynomial to define a multiplicative “slit loss” function for
each object that forces the observed continuum to match the best-fit SED.

After this additional flux correction step, the spectra are shifted into the
rest frame and normalized by the median observed continuum flux in the region
between $\lambda_{\rm rest}=6800-7000$ Å, where there are no emission lines;
this portion of the spectrum is also approximately centered with respect to
the auroral [S III]$\lambda 6313$ and [O II$]\lambda\lambda 7322,7332$ lines.
The spectra are then interpolated onto a common rest-frame wavelength array
and median-combined. The final stack is subsequently rescaled to match the
median rest-frame continuum of all the constituent galaxies between
$\lambda_{\rm rest}=6800-7000$ Å. Uncertainties are estimated by generating
1000 bootstrap-resampled composite spectra and calculating the $68\%$
confidence interval (CI; analogous to asymmetric error bars) at each
wavelength.

Figure 3 shows this composite spectrum (in medium blue) and the corresponding
uncertainties (in light blue) over the range of rest-wavelengths with
continuum S/N $\gtrsim 15$, where we define the S/N as the ratio of the
composite spectrum to half the $68\%$ CI. This requirement results in $\gtrsim
75\%$ ($\geq 17/23$ galaxies) contributing to the final composite at each
wavelength. At the center of the wavelength range where the targeted auroral
lines are found, the stack represents $\sim 690$ object-hours of exposure
time.

Aside from the “strong” H$\alpha$, [N II]$\lambda\lambda 6550,6585$, and [S
II]$\lambda\lambda 6718,6733$ lines (highlighted in the inset panel in Figure
3), no other emission lines are routinely detected in ground-based spectra of
individual $z\sim 2$ galaxies. Lines longward of $\sim 7000$ Å are virtually
inaccessible from the ground at $z\gtrsim 2$, due to a combination of the
rising thermal background in $K$-band and decreasing atmospheric transparency.
In more recent studies of high-$z$ galaxies using JWST/NIRSpec, emission lines
in this wavelength range that are fainter than nebular [N II] and [S II] are
only infrequently observed in individual galaxy spectra (e.g., Cameron et al.,
2023; Shapley et al., 2023; Sanders et al., 2023b)—and even these relatively
strong lines are not always visible in the spectra of some distant galaxies.
In the composite spectrum shown in Figure 3, we identify emission lines from
eight different elements (H, He, N, O, Si, S, Ar, and Ni, denoted by the
vertical dotted grey lines). Many of these have only rarely, if ever, been
observed outside of the nearby universe.

###  3.2 The Nebular Composite Spectrum

Figure 4: The nebular composite spectrum of the CECILIA sample is shown in
medium blue, with the light blue shading representing the 68% CI determined
via bootstrapping. The dark blue curve shows the best-fit model, which
includes emission lines of eight different elements, identified by the dashed
grey lines. The strengths of these lines relative to the narrow component of
H$\alpha$ are reported in Table 1. The residuals from the model are shown in
the bottom panels (medium blue), compared to the uncertainties on the median
stack (light blue).

To quantify the strength of these emission lines, we construct a second,
continuum-subtracted composite spectrum. In this case, after each spectrum is
flux-corrected to match the best-fit SED, the model continuum is subtracted
before the spectrum is shifted into the rest frame. To remove any remaining
irregular wavelength-dependent errors in the continuum subtraction, we
subtract a running median, using a large window ($\Delta\lambda_{\rm rest}\sim
200$ Å) to avoid over-correcting near the emission lines. The spectra are then
normalized by the measured flux in H$\alpha$ and median combined. Because
H$\alpha$ falls in the detector gap for three galaxies, the nebular composite
only includes 18 of the 23 galaxies used to construct the total composite
spectrum. Finally, the resulting composite spectrum is converted to flux units
($\lambda F_{\lambda}$) and rescaled so that the peak flux of (narrow)
H$\alpha$ is 100. Figure 4 shows this composite spectrum (in medium blue),
with the flux limit chosen to facilitate inspection of the semi-strong and
faint features. As before, uncertainties are estimated using bootstrapping
(shown by the light blue shading). The same lines identified in Figure 3 are
marked by dashed grey lines here.

We determine the typical strength of these emission lines by first fitting the
median composite with a model containing 73 emission lines, drawn from the
catalog reported by Esteban et al. (2004), who conducted a detailed analysis
of Very Large Telescope (VLT) UVES (Dekker et al., 2000) echelle
spectrophotometry of the Orion nebula. We select those lines in the wavelength
range sampled by the CECILIA nebular composite that are measured to have a
flux $>0.01\%$ of H$\alpha$ in the Esteban et al. (2004) spectrum. All of the
lines are modeled as single Gaussians, have fixed relative wavelengths (i.e.,
the line centers are not allowed to move relative to one another), and are
required to have the same width. For the strong [N II]$\lambda\lambda
6550,6585$ and semi-strong [O I]$\lambda\lambda 6302,6365$ doublets, which
have relative strengths set by atomic physics, the ratios are fixed at 1:2.96
and 3.15:1, respectively (Froese Fischer & Saha, 1983; Baluja & Zeippen, 1988;
Tachiev & Froese Fischer, 2001). A second Gaussian is included to account for
broad components under the strongest lines (H$\alpha$, [N II]$\lambda\lambda
6550,6585$, and [S II]$\lambda\lambda 6718,6733$) and allowed to be offset in
velocity relative to the narrow components of the same lines; all of the broad
components are required to have the same line width and velocity offset. The
addition of these components significantly improves the residuals from the
model by accounting for excess flux detected near the H$\alpha+[$N II]
complex.

The 1000 bootstrap samples are fit using the same model, and the $68\%$
highest density interval (HDI) for the distributions of measured fluxes are
used to determine uncertainties on the reported line fluxes. Lines are
considered well detected when they have a nonzero flux in $>99\%$ of fits
_and_ the maximum a posteriori (MAP) value for the line flux is $>3\sigma$.
Nineteen emission lines satisfy these criteria and are listed in Table 1. We
also include Si II $\lambda 6373$ (2.9$\sigma$), the weaker [Ni II] line at
7414 Å (2.7$\sigma$), the Paschen line at 8470 Å (3.1$\sigma$, but only
nonzero in 97.5% of the bootstrap stacks), and all of the broad components.

The dark blue curve in the top panels of Figure 4 represents the best-fit
model containing the emission lines in Table 1, with the fit residuals shown
by the medium blue line in the bottom panels. Recall that the peak of (narrow)
H$\alpha$ is set to 100 in the nebular composite, so that the peak of the
semi-strong and faint emission lines corresponds to their strengths relative
to the narrow component of H$\alpha$. The MAP values are also reported in
Table 1, along with the $68\%$ HDI for each line. Because the uncertainties
are calculated via bootstrap, note that these ranges reflect contributions
from both observational uncertainties on the individual line measurements and
physical variation among the objects in our sample.

##  4 Faint Emission Lines in High-redshift Star-forming Galaxies

Table 1: Observed Line Fluxes Relative to H$\alpha$.

Ion | $\lambda_{\rm vac}$ | F$(\lambda)$ | Range | Notes
---|---|---|---|---
| (Å) | (%) | (%) |
Narrow components
${\rm[N~{}II]}$ | 5756.24 | $0.20$ | $0.15-0.25$ | Auroral line
${\rm He~{}I}$ | 5877.27 | $4.12$ | $3.74-4.48$ |
${\rm[O~{}I]}$ | 6302.04 | $2.80$ | $2.36-3.08$ |
${\rm[S~{}III]}$ | 6313.85 | $0.34$ | $0.29-0.41$ | Auroral line
${\rm[O~{}I]}$ | 6365.54 | $0.91$ | $0.77-1.01$ |
${\rm Si~{}II}$ | 6373.12 | $0.12$ | $0.09-0.21$ |
${\rm[N~{}II]}$ | 6549.86 | $2.55$ | $1.86-3.02$ |
${\rm H\alpha}$ | 6564.62 | $100.00$ |  $\cdots$ $\cdots$ |
${\rm[N~{}II]}$ | 6585.27 | $7.59$ | $5.54-8.99$ |
${\rm He~{}I}$ | 6679.99 | $1.21$ | $1.09-1.34$ |
${\rm[S~{}II]}$ | 6718.29 | $7.33$ | $5.55-7.76$ |
${\rm[S~{}II]}$ | 6732.68 | $6.49$ | $4.94-6.93$ |
${\rm He~{}I}$ | 7067.23 | $1.25$ | $1.07-1.51$ |
${\rm[Ar~{}III]}$ | 7137.75 | $2.63$ | $2.22-2.86$ |
${\rm[O~{}II]}$ | 7321.94 | $1.21$ | $1.10-1.38$ | Auroral line
${\rm[O~{}II]}$ | 7332.21 | $0.95$ | $0.88-1.12$ | Auroral line
${\rm[Ni~{}II]}$ | 7379.86 | $0.22$ | $0.15-0.30$ |
${\rm[Ni~{}II]}$ | 7413.65 | $0.17$ | $0.11-0.24$ |
${\rm[Ar~{}III]}$ | 7753.23 | $0.64$ | $0.57-0.72$ |
${\rm P18}$ | 8440.28 | $0.20$ | $0.13-0.27$ |
${\rm O~{}I}$ | 8448.57 | $0.73$ | $0.64-0.87$ |
${\rm P17}$ | 8469.58 | $0.13$ | $0.09-0.18$ |
Broad components
${\rm[N~{}II]}$ | 6549.86 | $0.86$ | $0.15-1.56$ |
${\rm H\alpha}$ | 6564.62 | $12.63$ | $\phantom{0}6.01-28.31$ |
${\rm[N~{}II]}$ | 6585.27 | $2.55$ | $0.45-4.65$ |
${\rm[S~{}II]}$ | 6718.29 | $0.95$ | $0.00-2.48$ |
${\rm[S~{}II]}$ | 6732.68 | $0.81$ | $0.00-1.50$ |

Note. — [F(H$\alpha_{\rm narrow})=100$]

In this section, we highlight individual semi-strong ($\approx 2-3$% of
H$\alpha$) and faint ($\lesssim 1$% of H$\alpha$) emission lines detected in
the CECILIA composite spectra and briefly comment on how they may be used to
study high-$z$ galaxies.

###  4.1 Auroral Lines

Of all the faint lines present in the rest-optical spectra of star-forming
regions, the auroral121212“Auroral” lines are forbidden transitions from the
second excited state to the first excited state of ions of heavy elements and
can be paired with observations of the corresponding “nebular” (first excited
state to ground state) lines to determine $T_{e}$ in low-density gas where
collisional de-excitation does not play a significant role. lines that can be
used to implement the direct method of measuring metallicities have received
the most attention in studies of high-$z$ galaxies. Foremost among these is [O
III]$\lambda 4364$, which falls at $\lambda_{\rm obs}\approx 1.3-1.7~{}\mu$m
for the $z\sim 2-3$ CECILIA galaxies. As described in Section 2.2, CECILIA
instead targets two auroral lines at longer wavelengths that are not only
predicted to be stronger than [O III]$\lambda 4364$, but also fall at
$\lambda_{\rm obs}\gtrsim 2.0~{}\mu$m, where JWST is more sensitive. In total,
three auroral lines fall in the wavelength range sampled by the composite
spectra shown in Figures 3 and 4 and form the basis of our discussion here: [N
II]$\lambda 5756$, [S III]$\lambda 6313$ [O II]$\lambda\lambda 7322,7332$.

Figure 5: The three auroral lines observed in the CECILIA stack are shown in
separate panels with the same flux scale, in the same manner as Figure 4.
Auroral [N II]$\lambda 5756$ (left panel) is the weakest of the three and is
estimated to be $0.15-0.25\%$ the strength of H$\alpha$. Both [S III]$\lambda
6313$ (center panel) and [O II]$\lambda 7322,7332$ (right panel) are
noticeably stronger, and the oxygen lines are clearly the brightest auroral
features in the $\lambda=5700-8500$ Å range, with each being $\sim 1\%$ the
strength of H$\alpha$.

.

The strongest of these is [O II]$\lambda\lambda 7322,7332$, with both lines
observed to be $\sim 1\%$ the strength of the narrow component of H$\alpha$
(right panel of Figure 5). Sanders et al. (2023c) recently reported the
detection of this feature (actually a quadruplet) in two $z=2.18$ galaxies,
which had each been observed for $\sim 15$ hr using Keck/MOSFIRE. Using their
measurements to calculate direct-method oxygen abundances, they find moderate
$12+\log(\rm{O/H})=7.89\pm 0.20$ and $12+\log(\rm{O/H})=8.24\pm 0.27$.
Comparing these abundances to the predictions in the bottom panel of Figure 2,
we see that they lie near the broad peak of the predicted line strengths and,
thus, likely represent only the “tip of the iceberg”: other deep spectroscopic
studies should uncover auroral [O II] lines in galaxies with a wider range of
O/H. This is one of the main goals of the CECILIA program (Section 2.2), which
includes many high-confidence detections of these lines in individual galaxy
spectra that will be investigated in a subsequent paper.

The other auroral line specifically targeted by CECILIA is [S III]$\lambda
6313$ (middle panel of Figure 5), which samples a higher ionization zone than
auroral [O II]. Whereas this line is routinely used to measure abundances in
nearby extragalactic H II regions, it has never been reported in observations
of galaxies outside the local universe. It is significantly detected in the
CECILIA composites, and we have preliminary evidence of its presence in
spectra of individual CECILIA galaxies. However, because of its faintness
($0.29-0.41\%$ the strength of H$\alpha$) and proximity to the comparatively
stronger [O I]$\lambda 6302$ line, [S III]$\lambda 6313$ is unlikely to be
accessible in shallower or lower resolution JWST observations of objects
similar to the CECILIA sample.

Figure 6: Semi-strong and faint emission lines of four elements—O, Si, Ar, and
Ni—are shown in the same manner as Figures 4 and 5, but with a different flux
range used each panel. The top row shows forbidden [O I]$\lambda 6302,6565$
(upper left), with [S III]$\lambda 6313$ and Si II $\lambda 6373$ in the red
wings of each line, respectively, and permitted O I at $8449$ Å, blended with
Pa18 (upper right). The bottom row shows forbidden lines of [Ar
III]$\lambda\lambda 7138,7753$ (lower left) and [Ni II]$\lambda 7380,7414$
(lower right). The stronger [O I] and [Ar III] lines have now been observed in
deep spectra of individual high-$z$ galaxies (e.g., Cameron et al., 2023;
Sanders et al., 2023c), but the permitted O I $\lambda 8449$ line and
forbidden [Ni II] are only rarely observed, even in observations of nearby
galaxies and H II regions.

The third auroral line observed in the composite spectra is [N II]$\lambda
5756$, at $0.15-0.25\%$ the strength of H$\alpha$ (left panel of Figure 5). It
is formally detected at 4.6$\sigma$ in the stack but is likely too faint to be
detected in the individual spectra of high-$z$ galaxies, even using long
exposure times with JWST. Still, for the sample of galaxies where it _is_
possible to detect, it may serve as an important tool for calibrating the
$T_{e}$ relation between different ionization zones (e.g., Garnett, 1992;
Esteban et al., 2009; Croxall et al., 2016; Yates et al., 2020; Rogers et al.,
2021).

We use our measurement of [N II]$\lambda 5756$ to calculate $T_{e}$ and the
corresponding direct-method ionic abundance. First, we use the line strengths
for [S II]$\lambda\lambda 6718,6733$ to determine the electron density and
find $n_{e}\approx 285$ cm-3, which is consistent with values previously
reported for KBSS galaxies (Strom et al., 2017) and other $z\sim 2-3$ galaxy
samples (e.g., Sanders et al., 2016). This density is then combined with the
measurements of nebular and auroral [N II] to calculate $T_{e}$ using the
PyNeb package (Luridiana et al., 2015). However, because the nebular stack
does not contain both H$\alpha$ and H$\beta$, which are required to determine
the Balmer decrement and robustly constrain the reddening, we adopt the
interquartile range in E(B$-$V) for the KBSS parent sample,
E(B$-$V)$=0.06-0.47$ (Strom et al., 2017). Using these values, we find
$T_{e}$[N II]$=13630\pm 2540$ K, where the reported uncertainties also capture
the likely range in reddening for $z\sim 2-3$ galaxies. This temperature can
be used to infer the abundances of low-ionization species, and we ultimately
calculate $12+\log($N+/H${}^{+})=6.33^{+0.18}_{-0.30}$ and
$12+\log($S+/H${}^{+})=5.70^{+0.16}_{-0.26}$ using the nebular lines for both
ions; because the auroral [O II]$\lambda\lambda 7322,7332$ lines are, by
definition, strong functions of $T_{e}$, any O+/H+ abundance determined using
$T_{e}$[N II] would have much larger uncertainties. We cannot confidently
determine the contribution from other common (but unseen) ionization states of
N and S using the nebular composite alone and so do not report total
abundances here, but we plan to revisit the issue of appropriate ionization
correction factors for high-$z$ galaxies in future work.

###  4.2 Other Semi-strong and Faint Lines

In addition to the three auroral emission lines, we also detect eight semi-
strong and faint emission lines from four different heavy elements. Cutouts of
the nebular composite spectrum near these features are shown in Figure 6, in
the same manner as Figures 4 and 5.

The strongest of the lines is [O I]$\lambda 6302$ (upper left panel of Figure
6), which is observed to be $2.36-3.08\%$ the strength of the narrow component
of H$\alpha$ and is significantly stronger than the auroral [S III]$\lambda
6313$ line in its red wing. Its partner line at 6365 Å is $\approx 3.15\times$
weaker as set by atomic physics and is blended with Si II $\lambda 6373$,
which also probes mostly neutral and low-ionization gas. Rather than being an
abundance diagnostic, [O I]$\lambda 6302$ is most commonly used as a way to
identify the principal ionization mechanism in emission line galaxies using a
form of the Baldwin-Philips-Terlevich (BPT) diagram (Baldwin et al., 1981;
Veilleux & Osterbrock, 1987). It can also be useful for identifying
contributions from diffuse ionized gas and shocks (e.g., Tüllmann & Dettmar,
2000; Moy & Rocca-Volmerange, 2002; Zhang et al., 2017). Although widely
studied in the local universe, including in large samples such as the Sloan
Digital Sky Survey (SDSS; Kewley et al., 2006; Law et al., 2021), [O
I]$\lambda 6302$ is not commonly reported for individual $z\gtrsim 2$
galaxies. More recently, however, it has been observed in a handful of
high-$z$ galaxies with moderately deep ground-based or JWST spectroscopy
(e.g., Cameron et al., 2023; Clarke et al., 2023; Sanders et al., 2023c).
Given its typical strength, this semi-strong line should provide an accessible
and promising method for discriminating between AGN and star formation and
probing low-ionization gas in high-$z$ galaxies.

The lower left panel of Figure 6 shows the widely-spaced [Ar
III]$\lambda\lambda 7138,7753$ lines, which, like lines of O++, trace the gas
in the high ionization zone of star-forming regions. In low-$z$ galaxies, the
stronger [Ar III]$\lambda 7138$ line has been used to determine absolute argon
abundances and relative abundance ratios, such as Ar/O (e.g., Berg et al.,
2015; Croxall et al., 2016; Rogers et al., 2021), after accounting for unseen
ionization states of Ar. At $2.22-2.86$% the strength of H$\alpha$, comparable
to [O I]$\lambda 6302$, this line is one of the strongest heavy metal lines
present in the CECILIA nebular composite spectrum, aside from the familiar
strong lines; Sanders et al. (2023c) find similar ratios for their two
$z=2.18$ galaxies. Thus, [Ar III]$\lambda 7138$ is also an attractive target
for spectroscopic studies of galaxy enrichment. Although both Ar and O are
nominally produced by the same mechanism in massive stars, differences in Ar/O
as a function of overall enrichment could reflect a dependence of stellar
nucleosynthesis on metallicity (e.g., Kennicutt et al., 2003; Izotov et al.,
2006).

Permitted O I at 8449 Å (upper right panel of Figure 6) is one of the most
commonly used recombination lines for measuring metallicity in astrophysical
environments where such transitions can be observed (Maiolino & Mannucci,
2019). Typically, however, metal recombination lines are too weak to be useful
diagnostics, even in $z\sim 0$ galaxies and H II regions, so its presence in
the CECILIA composite is unexpected. We consulted the database maintained by
the Atomic Spectroscopy Data Center at the National Institute of Standards and
Technology (NIST), but no other likely candidates for emission lines at the
same rest wavelength were identified. O I $\lambda 8449$ is blended with the
Paschen series line at 8440 Å (P18) in its left wing, and the neighboring
Paschen line at 8470 Å (P17) is also detected in the nebular composite,
relieving concerns that poor wavelength calibration may have led to
misidentifying the line. The anomalous strength of this line in the Orion
Nebula is in fact a longstanding puzzle (e.g., Morgan, 1971; Danziger &
Aaronson, 1974), but Grandi (1975) showed that direct excitation by starlight
is the most likely origin, rather than recombination from O+ or Ly$\beta$
fluorescence. The strength of O I $\lambda 8449$ in the CECILIA composite
relative to H$\alpha$ ($0.64-0.87$%) is comparable to that reported for the
Orion Nebular ($0.56$%; Esteban et al., 2004), suggesting that the same
mechanism may also be dominant in high-$z$ galaxies. If true, the contribution
from direct excitation would severely limit the utility of O I $\lambda 8449$
as a metallicity tracer; García-Rojas et al. (2006) and García-Rojas et al.
(2007) used both O I $\lambda 7771$ and O I $\lambda 8449$ to derive O+
abundances in Galactic H II regions but found an order of magnitude larger
O+/H+ using O I $\lambda 8449$, likely due to this effect.

Also puzzling is the detection of [Ni II]$\lambda\lambda 7380,7414$ (shown in
the lower right panel of Figure 6), which are $0.15-0.30$% and $0.11-0.24$%
the strength of H$\alpha$, respectively. There are comparatively few
references to the observation of this line in astrophysical objects, but a
handful of studies have reported measurements of [Ni II]$\lambda 7380$ and the
corresponding Ni+ abundances in gaseous nebulae in the Milky Way (Dennefeld,
1982; Fesen & Kirshner, 1982; Henry & Fesen, 1988; Esteban et al., 1999). In
many of these cases, [Ni II]$\lambda 7380$ was seen to be much stronger than
expected relative to the associated [Ni II]$\lambda 7414$ line, which is only
marginally detected in the CECILIA composite spectra. Other authors have
explained this by invoking fluorescence by the UV continuum, similar to [Fe
II] lines (Lucy, 1995), and/or collisional excitation in very high density
($n_{e}\approx 10^{6}$ cm-3) gas (Bautista et al., 1996). Similar to O I
$\lambda 8449$, no reasonable alternatives were identified in the NIST
database, but without significantly detecting [Fe II] lines that may also be
impacted by the same physical mechanism, it is difficult to speculate about
its appearance here.

###  4.3 Broad H$\alpha$

Figure 7: The same nebular composite spectrum from Figure 4 is shown near the
strong H$\alpha$, [N II]$\lambda\lambda 6550,6585$, and [S II]$\lambda\lambda
6718,6733$ lines, with a more extended flux range to show the peaks of the
lines. The inset panel is a zoom-in on the grey shaded region, highlighting
the broad components of H$\alpha$ and [N II] (in red). The narrow components
of all three lines are illustrated by the dot-dashed navy curves.

Broad line emission has been observed in both spectra of individual galaxies
and composite spectra of galaxies at $z\sim 2$ and is usually attributed to
galaxy-scale ionized gas outflows (see Section 4.6 of the review by Förster
Schreiber & Wuyts, 2020, and references therein); in contrast, the frequently
brighter, narrow components of emission lines trace galaxies’ star-forming
regions. Both active galactic nuclei (AGN; e.g., Nesvadba et al., 2008; Genzel
et al., 2014; Förster Schreiber et al., 2014; Cresci et al., 2015) and star-
formation (e.g., Genzel et al., 2011; Davies et al., 2019; Freeman et al.,
2019) can generate these outflows, resulting in differences in inferred
outflow velocity (i.e., emission line width), with AGN typically driving
higher velocity outflows than feedback from massive stars.

In order to achieve a good fit to the nebular composite spectrum, we include
two Gaussian components for the strongest lines to account for excess flux
that results in large residuals from a model with only a single (narrow)
component. Based on the results from fitting the 1000 bootstrap samples, these
broad components have a FWHM of $536^{+45}_{-167}$ km s-1 and are consistent
with no velocity offset relative to the narrow components, which have a FWHM
of $288^{+15}_{-20}$ km s-1. The broad H$\alpha$ line is $6.01-28.31$% the
strength of the narrow component of H$\alpha$, with significantly weaker broad
components observed for nebular [N II] and [S II]; this is consistent with the
low end of the range reported for a similar sample of $z\sim 2$ star-forming
galaxies (Freeman et al., 2019). If these components in the CECILIA composite
do indeed reflect the presence of ionized gas outflows, the evidence for broad
(albeit weak) line emission in forbidden transitions and the moderate velocity
width suggest that they are likely driven by star formation. Comparable FWHM
velocities of $\sim 400-500$ km s-1 are observed in deep VLT/SINFONI
(Eisenhauer et al., 2003; Bonnet et al., 2004) spectra of star-forming clumps
at $z\sim 2$ (e.g., Newman et al., 2012; Förster Schreiber et al., 2019).
However, because of the additional median filtering required to remove
remaining fluctuations in the continuum of individual galaxy spectra (Section
3.2), the detailed properties of any broad line emission in the CECILIA stack
could be systematically biased.

##  5 Conclusions

We have reported the first results from the CECILIA program (JWST PID 2593),
which obtained ultra-deep $\sim 30$ hr NIRSpec/G235M observations of 33 star-
forming galaxies at $z\sim 1-3$. Using data for 23 of these galaxies, we
constructed rest-optical composite spectra, both with and without the stellar
continuum, corresponding to exposure times of 690 object-hours and 540 object-
hours, respectively. These composites, shown in Figures 3 and 4, provide one
of the most detailed views to date of star-forming galaxies in the early
universe and function as an atlas of their characteristic rest-optical
emission line spectra.

The principal findings based on our analysis of the stacked spectra are as
follows:

  * •

We significantly detect emission lines of eight different elements (H, He, N,
O, Si, S, Ar, and Ni), including evidence for broad line emission under
H$\alpha$, [N II]$\lambda\lambda 6550,6585$, and [S II]$\lambda\lambda
6718,6733$. The strengths of these lines relative to the narrow component of
H$\alpha$ are reported in Table 1.

  * •

Aside from strong [N II], H$\alpha$, and [S II], which have previously been
studied in large ground-based spectroscopic samples, the majority of emission
lines are $\lesssim 3\%$ the strength of H$\alpha$. Some of these features,
such as [O I]$\lambda 6302$ (shown in the upper left panel of Figure 6), are
now being detected in JWST spectra of individual high-$z$ galaxies, and we
expect other lines with strengths $\gtrsim 2-3$% that of H$\alpha$ to be good
candidates for spectroscopic follow-up of large samples. In addition to the
stronger forbidden [O I] line, these semi-strong lines include the He I line
at $\lambda 5877$ and [Ar III]$\lambda 7138$ (shown in the bottom left panel
of Figure 6).

  * •

The three auroral emission lines present at $\lambda_{\rm rest}\approx
5700-8500$ ([N II]$\lambda 5756$, [S III]$\lambda 6313$, [O II]$\lambda\lambda
7322,7332$, shown in Figure 5) are $\lesssim 1\%$ the strength of H$\alpha$.
Using our measurements of auroral and nebular [N II], we find $T_{e}$[N II]
$=13630\pm 2540$ K, which is the first time a $T_{e}$ has been reported for
high-redshift galaxies using this tracer. Although we have not reported the
significance of detections in individual galaxy spectra in this work, it seems
likely that these auroral lines will remain out of reach of typical
observations of high-$z$ galaxies, particularly those with low SFRs, low
ionization, and/or high metallicity. This only underscores the need for more
accurate line ratio diagnostics for metallicity that make use of the strong
and semi-strong emission lines present in galaxies’ rest-optical spectra.

  * •

We measure broad ($536^{+45}_{-167}$ km s-1 FWHM) line emission under the
strongest lines and a broad component of H$\alpha$ that is $6.01-28.31$% the
strength of the narrow component (Figure 7). These results appear indicative
of star-formation driven outflows. However, we caution that, owing to
remaining uncertainties in the flux calibration (see the discussion in Section
2.4), the appearance of this component should not be over-interpreted. We
defer a more detailed discussion of broad line emission and its connection to
galaxy outflows to a future paper.

JWST is delivering on its promise to provide access to faint emission lines in
the spectra of $\gtrsim 2$ galaxies, evidenced not only by what we have
presented in this let, but also by the many exciting results based on
NIRSpec/MSA and NIRCam grism spectroscopy that have been published over the
last year. Deep observations, such as those obtained as part of CECILIA and
outlined here, will be critical for developing and testing the new tools
necessary to accurately interpret this wealth of data. As known issues with
JWST data products continue to be resolved, it will benefit the extragalactic
community to revisit some of the earliest observations—with the benefit of
hindsight and these new tools—in order to maximize the scientific impact of
these data. To aid in this effort, forthcoming work with CECILIA will focus on
(1) $T_{e}$ measurements and direct-method metallicities for the sample of
galaxies introduced here, as well as (2) new line ratio diagnostics for gas-
phase oxygen abundance.

The authors thank Jane Rigby, Taylor Hutchison, and Marcia Rieke for their
advice regarding the reduction of the JWST data, as well as Jenny Greene for
her input on the scope of the discussion. We are also grateful to the
JWST/NIRSpec team for their ongoing work to support this complex and powerful
instrument. ALS, GCR and RFT acknowledge partial support from the JWST-
GO-02593.008-A, JWST-GO-02593.004-A, and JWST-GO-02593.006-A grants,
respectively. RFT also acknowledges support from the Pittsburgh Foundation
(grant ID UN2021-121482) and the Research Corporation for Scientific
Advancement (Cottrell Scholar Award, grant ID 28289). This work is primarily
based on observations made with NASA/ESA/CSA JWST, associated with PID 2593,
which can be accessed via doi:https://doi.org/10.17909/x66z-p144 (catalog
10.17909/x66z-p144). The data were obtained from the Mikulski Archive for
Space Telescopes (MAST) at the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
under NASA contract NAS 5-03127 for JWST. Some of the data used to generate
the original line flux predictions were obtained at the W.M. Keck Observatory,
which is operated as a scientific partnership between the California Institute
of Technology, the University of California, and NASA. The Observatory was
made possible by the generous financial support of the W.M. Keck Foundation,
and the authors wish to recognize and acknowledge the significant cultural
role and reverence that the summit of Maunakea has within the indigenous
Hawaiian community.

## References

  * Arellano-Córdova et al. (2022) Arellano-Córdova, K. Z., Berg, D. A., Chisholm, J., et al. 2022, ApJ, 940, L23, doi: 10.3847/2041-8213/ac9ab2
  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5, doi: 10.1086/130766
  * Baluja & Zeippen (1988) Baluja, K. L., & Zeippen, C. J. 1988, Journal of Physics B Atomic Molecular Physics, 21, 1455, doi: 10.1088/0953-4075/21/9/007
  * Bautista et al. (1996) Bautista, M. A., Peng, J., & Pradhan, A. K. 1996, ApJ, 460, 372, doi: 10.1086/176976
  * Belfiore et al. (2015) Belfiore, F., Maiolino, R., Bundy, K., et al. 2015, MNRAS, 449, 867, doi: 10.1093/mnras/stv296
  * Berg et al. (2015) Berg, D. A., Skillman, E. D., Croxall, K. V., et al. 2015, ApJ, 806, 16, doi: 10.1088/0004-637X/806/1/16
  * Bonnet et al. (2004) Bonnet, H., Conzelmann, R., Delabre, B., et al. 2004, in SPIE Conference Series, Vol. 5490, Advancements in Adaptive Optics, ed. D. Bonaccini Calia, B. L. Ellerbroek, & R. Ragazzoni, 130–138, doi: 10.1117/12.551187
  * Brammer (2022) Brammer, G. 2022, msaexp: NIRSpec analyis tools, 0.3.4, Zenodo, Zenodo, doi: 10.5281/zenodo.7299500
  * Brammer (2023) —. 2023, grizli, 1.8.2, Zenodo, Zenodo, doi: 10.5281/zenodo.7712834
  * Brinchmann (2023) Brinchmann, J. 2023, MNRAS, doi: 10.1093/mnras/stad1704
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151, doi: 10.1111/j.1365-2966.2004.07881.x
  * Bushouse et al. (2023) Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023, JWST Calibration Pipeline, 1.11.0, Zenodo, doi: 10.5281/zenodo.8067394
  * Cameron et al. (2023) Cameron, A. J., Saxena, A., Bunker, A. J., et al. 2023, A&A, 677, A115, doi: 10.1051/0004-6361/202346107
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Christensen et al. (2012) Christensen, L., Laursen, P., Richard, J., et al. 2012, MNRAS, 427, 1973, doi: 10.1111/j.1365-2966.2012.22007.x
  * Clarke et al. (2023) Clarke, L., Shapley, A., Sanders, R. L., et al. 2023, arXiv e-prints, arXiv:2305.07781, doi: 10.48550/arXiv.2305.07781
  * Cresci et al. (2015) Cresci, G., Mainieri, V., Brusa, M., et al. 2015, ApJ, 799, 82, doi: 10.1088/0004-637X/799/1/82
  * Croxall et al. (2016) Croxall, K. V., Pogge, R. W., Berg, D. A., Skillman, E. D., & Moustakas, J. 2016, ApJ, 830, 4, doi: 10.3847/0004-637X/830/1/4
  * Curti et al. (2023) Curti, M., D’Eugenio, F., Carniani, S., et al. 2023, MNRAS, 518, 425, doi: 10.1093/mnras/stac2737
  * Danziger & Aaronson (1974) Danziger, I. J., & Aaronson, M. 1974, PASP, 86, 208, doi: 10.1086/129584
  * Davies et al. (2019) Davies, R. L., Förster Schreiber, N. M., Übler, H., et al. 2019, ApJ, 873, 122, doi: 10.3847/1538-4357/ab06f1
  * Dekker et al. (2000) Dekker, H., D’Odorico, S., Kaufer, A., Delabre, B., & Kotzlowski, H. 2000, in SPIE Conference Series, Vol. 4008, Optical and IR Telescope Instrumentation and Detectors, ed. M. Iye & A. F. Moorwood, 534–545, doi: 10.1117/12.395512
  * Dennefeld (1982) Dennefeld, M. 1982, A&A, 112, 215
  * Eisenhauer et al. (2003) Eisenhauer, F., Abuter, R., Bickert, K., et al. 2003, in SPIE Conference Series, Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, ed. M. Iye & A. F. M. Moorwood, 1548–1561, doi: 10.1117/12.459468
  * Eldridge et al. (2017) Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058, doi: 10.1017/pasa.2017.51
  * Esteban et al. (2009) Esteban, C., Bresolin, F., Peimbert, M., et al. 2009, ApJ, 700, 654, doi: 10.1088/0004-637X/700/1/654
  * Esteban et al. (2004) Esteban, C., Peimbert, M., García-Rojas, J., et al. 2004, MNRAS, 355, 229, doi: 10.1111/j.1365-2966.2004.08313.x
  * Esteban et al. (1999) Esteban, C., Peimbert, M., Torres-Peimbert, S., García-Rojas, J., & Rodríguez, M. 1999, ApJS, 120, 113, doi: 10.1086/313172
  * Ferland et al. (2013) Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mexicana Astron. Astrofis., 49, 137
  * Fesen & Kirshner (1982) Fesen, R. A., & Kirshner, R. P. 1982, ApJ, 258, 1, doi: 10.1086/160043
  * Förster Schreiber & Wuyts (2020) Förster Schreiber, N. M., & Wuyts, S. 2020, ARA&A, 58, 661, doi: 10.1146/annurev-astro-032620-021910
  * Förster Schreiber et al. (2014) Förster Schreiber, N. M., Genzel, R., Newman, S. F., et al. 2014, ApJ, 787, 38, doi: 10.1088/0004-637X/787/1/38
  * Förster Schreiber et al. (2019) Förster Schreiber, N. M., Übler, H., Davies, R. L., et al. 2019, ApJ, 875, 21, doi: 10.3847/1538-4357/ab0ca2
  * Freeman et al. (2019) Freeman, W. R., Siana, B., Kriek, M., et al. 2019, ApJ, 873, 102, doi: 10.3847/1538-4357/ab0655
  * Froese Fischer & Saha (1983) Froese Fischer, C., & Saha, H. P. 1983, Phys. Rev. A, 28, 3169, doi: 10.1103/PhysRevA.28.3169
  * García-Rojas et al. (2007) García-Rojas, J., Esteban, C., Peimbert, A., et al. 2007, Rev. Mexicana Astron. Astrofis., 43, 3, doi: 10.48550/arXiv.astro-ph/0610065
  * García-Rojas et al. (2006) García-Rojas, J., Esteban, C., Peimbert, M., et al. 2006, MNRAS, 368, 253, doi: 10.1111/j.1365-2966.2006.10105.x
  * Garnett (1992) Garnett, D. R. 1992, AJ, 103, 1330, doi: 10.1086/116146
  * Genzel et al. (2011) Genzel, R., Newman, S., Jones, T., et al. 2011, ApJ, 733, 101, doi: 10.1088/0004-637X/733/2/101
  * Genzel et al. (2014) Genzel, R., Förster Schreiber, N. M., Rosario, D., et al. 2014, ApJ, 796, 7, doi: 10.1088/0004-637X/796/1/7
  * Grandi (1975) Grandi, S. A. 1975, ApJ, 196, 465, doi: 10.1086/153426
  * Henry & Fesen (1988) Henry, R. B. C., & Fesen, R. A. 1988, ApJ, 329, 693, doi: 10.1086/166413
  * Izotov et al. (2006) Izotov, Y. I., Stasińska, G., Meynet, G., Guseva, N. G., & Thuan, T. X. 2006, A&A, 448, 955, doi: 10.1051/0004-6361:20053763
  * James et al. (2014) James, B. L., Pettini, M., Christensen, L., et al. 2014, MNRAS, 440, 1794, doi: 10.1093/mnras/stu287
  * Kashino et al. (2023) Kashino, D., Lilly, S. J., Matthee, J., et al. 2023, ApJ, 950, 66, doi: 10.3847/1538-4357/acc588
  * Katz et al. (2023) Katz, H., Saxena, A., Cameron, A. J., et al. 2023, MNRAS, 518, 592, doi: 10.1093/mnras/stac2657
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055, doi: 10.1111/j.1365-2966.2003.07154.x
  * Kennicutt et al. (2003) Kennicutt, Jr., R. C., Bresolin, F., & Garnett, D. R. 2003, ApJ, 591, 801, doi: 10.1086/375398
  * Kewley et al. (2006) Kewley, L. J., Groves, B., Kauffmann, G., & Heckman, T. 2006, MNRAS, 372, 961, doi: 10.1111/j.1365-2966.2006.10859.x
  * Kewley et al. (2019) Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511, doi: 10.1146/annurev-astro-081817-051832
  * Kocevski et al. (2023) Kocevski, D. D., Onoue, M., Inayoshi, K., et al. 2023, ApJ, 954, L4, doi: 10.3847/2041-8213/ace5a0
  * Laseter et al. (2023) Laseter, I. H., Maseda, M. V., Curti, M., et al. 2023, arXiv e-prints, arXiv:2306.03120, doi: 10.48550/arXiv.2306.03120
  * Law et al. (2012) Law, D. R., Steidel, C. C., Shapley, A. E., et al. 2012, ApJ, 745, 85, doi: 10.1088/0004-637X/745/1/85
  * Law et al. (2021) Law, D. R., Ji, X., Belfiore, F., et al. 2021, ApJ, 915, 35, doi: 10.3847/1538-4357/abfe0a
  * Lucy (1995) Lucy, L. B. 1995, A&A, 294, 555
  * Luridiana et al. (2015) Luridiana, V., Morisset, C., & Shaw, R. A. 2015, A&A, 573, A42, doi: 10.1051/0004-6361/201323152
  * Maiolino & Mannucci (2019) Maiolino, R., & Mannucci, F. 2019, A&A Rev., 27, 3, doi: 10.1007/s00159-018-0112-2
  * Masters et al. (2014) Masters, D., McCarthy, P., Siana, B., et al. 2014, ApJ, 785, 153, doi: 10.1088/0004-637X/785/2/153
  * McLean et al. (2010) McLean, I. S., Steidel, C. C., Epps, H., et al. 2010, in SPIE Conference Series, Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 77351E, doi: 10.1117/12.856715
  * McLean et al. (2012) McLean, I. S., Steidel, C. C., Epps, H. W., et al. 2012, in SPIE Conference Series, Vol. 8446, doi: 10.1117/12.924794
  * Mingozzi et al. (2020) Mingozzi, M., Belfiore, F., Cresci, G., et al. 2020, A&A, 636, A42, doi: 10.1051/0004-6361/201937203
  * Morgan (1971) Morgan, L. A. 1971, MNRAS, 153, 393, doi: 10.1093/mnras/153.4.393
  * Moy & Rocca-Volmerange (2002) Moy, E., & Rocca-Volmerange, B. 2002, A&A, 383, 46, doi: 10.1051/0004-6361:20011727
  * Nesvadba et al. (2008) Nesvadba, N. P. H., Lehnert, M. D., De Breuck, C., Gilbert, A. M., & van Breugel, W. 2008, A&A, 491, 407, doi: 10.1051/0004-6361:200810346
  * Newman et al. (2012) Newman, S. F., Shapiro Griffin, K., Genzel, R., et al. 2012, ApJ, 752, 111, doi: 10.1088/0004-637X/752/2/111
  * Oesch et al. (2023) Oesch, P. A., Brammer, G., Naidu, R. P., et al. 2023, MNRAS, 525, 2864, doi: 10.1093/mnras/stad2411
  * Pontoppidan et al. (2022) Pontoppidan, K. M., Barrientes, J., Blome, C., et al. 2022, ApJ, 936, L14, doi: 10.3847/2041-8213/ac8a4e
  * Rauscher (2023) Rauscher, B. J. 2023, arXiv e-prints, arXiv:2306.03250, doi: 10.48550/arXiv.2306.03250
  * Rhoads et al. (2023) Rhoads, J. E., Wold, I. G. B., Harish, S., et al. 2023, ApJ, 942, L14, doi: 10.3847/2041-8213/acaaaf
  * Rogers et al. (2021) Rogers, N. S. J., Skillman, E. D., Pogge, R. W., et al. 2021, ApJ, 915, 21, doi: 10.3847/1538-4357/abf8b9
  * Rudie et al. (2012) Rudie, G. C., Steidel, C. C., Trainor, R. F., et al. 2012, ApJ, 750, 67, doi: 10.1088/0004-637X/750/1/67
  * Sanders et al. (2023a) Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2023a, arXiv e-prints, arXiv:2303.08149, doi: 10.48550/arXiv.2303.08149
  * Sanders et al. (2023b) —. 2023b, ApJ, 955, 54, doi: 10.3847/1538-4357/acedad
  * Sanders et al. (2015) Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2015, ApJ, 799, 138, doi: 10.1088/0004-637X/799/2/138
  * Sanders et al. (2016) —. 2016, ApJ, 816, 23, doi: 10.3847/0004-637X/816/1/23
  * Sanders et al. (2020) Sanders, R. L., Shapley, A. E., Reddy, N. A., et al. 2020, MNRAS, 491, 1427, doi: 10.1093/mnras/stz3032
  * Sanders et al. (2023c) Sanders, R. L., Shapley, A. E., Clarke, L., et al. 2023c, ApJ, 943, 75, doi: 10.3847/1538-4357/aca9cc
  * Schaerer et al. (2022) Schaerer, D., Marques-Chaves, R., Barrufet, L., et al. 2022, A&A, 665, L4, doi: 10.1051/0004-6361/202244556
  * Shapley et al. (2023) Shapley, A. E., Reddy, N. A., Sanders, R. L., Topping, M. W., & Brammer, G. B. 2023, ApJ, 950, L1, doi: 10.3847/2041-8213/acd939
  * Shapley et al. (2015) Shapley, A. E., Reddy, N. A., Kriek, M., et al. 2015, ApJ, 801, 88, doi: 10.1088/0004-637X/801/2/88
  * Sharples et al. (2013) Sharples, R., Bender, R., Agudo Berbel, A., et al. 2013, The Messenger, 151, 21
  * Simcoe et al. (2008) Simcoe, R. A., Burgasser, A. J., Bernstein, R. A., et al. 2008, in Proc. SPIE, Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II, 70140U, doi: 10.1117/12.790414
  * Simcoe et al. (2010) Simcoe, R. A., Burgasser, A. J., Bochanski, J. J., et al. 2010, in SPIE Conference Series, Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 773514, doi: 10.1117/12.856088
  * Stanway et al. (2016) Stanway, E. R., Eldridge, J. J., & Becker, G. D. 2016, MNRAS, 456, 485, doi: 10.1093/mnras/stv2661
  * Steidel et al. (2010) Steidel, C. C., Erb, D. K., Shapley, A. E., et al. 2010, ApJ, 717, 289, doi: 10.1088/0004-637X/717/1/289
  * Steidel et al. (2016) Steidel, C. C., Strom, A. L., Pettini, M., et al. 2016, ApJ, 826, 159, doi: 10.3847/0004-637X/826/2/159
  * Steidel et al. (2014) Steidel, C. C., Rudie, G. C., Strom, A. L., et al. 2014, ApJ, 795, 165, doi: 10.1088/0004-637X/795/2/165
  * Strom et al. (2021) Strom, A. L., Rudie, G. C., Chen, Y., et al. 2021, CECILIA: A direct-method metallicity calibration for Cosmic Noon through the Epoch of Reionization, JWST Proposal. Cycle 1, ID. #2593
  * Strom et al. (2018) Strom, A. L., Steidel, C. C., Rudie, G. C., Trainor, R. F., & Pettini, M. 2018, ApJ, 868, 117, doi: 10.3847/1538-4357/aae1a5
  * Strom et al. (2017) Strom, A. L., Steidel, C. C., Rudie, G. C., et al. 2017, ApJ, 836, 164, doi: 10.3847/1538-4357/836/2/164
  * Sun et al. (2023) Sun, F., Egami, E., Pirzkal, N., et al. 2023, ApJ, 953, 53, doi: 10.3847/1538-4357/acd53c
  * Tachiev & Froese Fischer (2001) Tachiev, G., & Froese Fischer, C. 2001, Canadian Journal of Physics, 79, 955, doi: 10.1139/p01-059
  * Taylor et al. (2022) Taylor, A. J., Barger, A. J., & Cowie, L. L. 2022, ApJ, 939, L3, doi: 10.3847/2041-8213/ac959d
  * Topping et al. (2020) Topping, M. W., Shapley, A. E., Reddy, N. A., et al. 2020, MNRAS, 499, 1652, doi: 10.1093/mnras/staa2941
  * Trainor et al. (2015) Trainor, R. F., Steidel, C. C., Strom, A. L., & Rudie, G. C. 2015, ApJ, 809, 89, doi: 10.1088/0004-637X/809/1/89
  * Trainor et al. (2016) Trainor, R. F., Strom, A. L., Steidel, C. C., & Rudie, G. C. 2016, ApJ, 832, 171, doi: 10.3847/0004-637X/832/2/171
  * Tremonti et al. (2004) Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898, doi: 10.1086/423264
  * Trump et al. (2023) Trump, J. R., Arrabal Haro, P., Simons, R. C., et al. 2023, ApJ, 945, 35, doi: 10.3847/1538-4357/acba8a
  * Trussler et al. (2023) Trussler, J. A. A., Adams, N. J., Conselice, C. J., et al. 2023, MNRAS, 523, 3423, doi: 10.1093/mnras/stad1629
  * Tüllmann & Dettmar (2000) Tüllmann, R., & Dettmar, R. J. 2000, A&A, 362, 119, doi: 10.48550/arXiv.astro-ph/0008340
  * Veilleux & Osterbrock (1987) Veilleux, S., & Osterbrock, D. E. 1987, ApJS, 63, 295, doi: 10.1086/191166
  * Wisnioski et al. (2015) Wisnioski, E., Förster Schreiber, N. M., Wuyts, S., et al. 2015, ApJ, 799, 209, doi: 10.1088/0004-637X/799/2/209
  * Yates et al. (2020) Yates, R. M., Schady, P., Chen, T. W., Schweyer, T., & Wiseman, P. 2020, A&A, 634, A107, doi: 10.1051/0004-6361/201936506
  * Zhang et al. (2017) Zhang, K., Yan, R., Bundy, K., et al. 2017, MNRAS, 466, 3217, doi: 10.1093/mnras/stw3308