CECILIA: The Faint Emission Line Spectrum of z ∼ 2–3 Star-forming Galaxies

We present the first results from Chemical Evolution Constrained Using Ionized Lines in Interstellar Aurorae (CECILIA), a Cycle 1 JWST NIRSpec/MSA program that uses ultra-deep ∼30 hr G235M/F170LP observations to target multiple electron temperature-sensitive auroral lines in the spectra of 33 galaxies at z ∼ 1–3. Using a subset of 23 galaxies, we construct two ∼600 object-hour composite spectra, both with and without the stellar continuum, and use these to investigate the characteristic rest-optical (λ rest ≈ 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 ≲3% the strength of Hα. We report the characteristic strength of three auroral features ([N ii]λ5756, [S iii]λ6313, and [O ii]λ λ7322, 7332), as well as other semi-strong and faint emission lines, including forbidden [Ni ii]λ λ7380, 7414 and permitted O i λ8449, some of which have never before been observed outside of the local Universe. Using these measurements, we find T e [N ii] = 13,630 ± 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−167+45 km s−1; the broad component of Hα 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.


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 ele-ments 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.
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 groundbased telescopes in the world, measurements of auroral [O III]λ4364 were only possible for a handful of individual galaxies at z ≳ 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.37327field reinforced this expectation, with significant detections of [O III]λ4364 in several z ∼ 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(PIDs , 1914(PIDs , and 2593) ) by the time allocation committee, with a total investment of over 150 hrs, or ∼ 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: 10.17909/x66z-p144).
CECILIA was designed to measure auroral [S III]λ6313 and [O II]λλ7322, 7332 in the spectra of a carefully selected sample of z ∼ 2 − 3 star-forming galaxies, using ∼ 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 ∼ 600 object-hour rest-optical composite spectra of z ∼ 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 ≳ 2 galaxies.Throughout the text, we refer to specific spectral features using their vacuum wavelengths.

THE CECILIA SURVEY
The principal goal of CECILIA is to measure multiple faint rest-optical auroral lines in the spectra of z ∼ 2 − 3 galaxies, which can then be used to calibrate new highz 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 (∼ 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)1 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.

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 narrowband Lyα selected galaxies in 15 fields.The survey includes deep J, H, and K NIR (rest-optical) spectroscopy from Keck/MOSFIRE, deep optical (rest-UV) spec-troscopy from Keck/LRIS, and imaging in U n through K s , HST /WFC3 F140W, F160W, Spitzer IRAC Ch1−4 and MIPS 24 µm, and narrow-band (rest-frame) Lyα filters.In total, KBSS comprises ∼ 3500 z ∼ 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-colorselected galaxies for which detecting the auroral [S III]λ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 ≤ z ≤ 2.68, where all of the required emission lines fall within the G235M, G395M, and groundbased NIR spectral bandpasses.Galaxies were prioritized if they have: somewhat higher redshifts (z > 2.3), where NIRSpec is more sensitive at the observed wave-length of [S III]λ6313; smaller than average sizes, increasing the auroral surface brightness; high SED-based SFRs (> 24 M ⊙ yr −1 ); and/or large observed nebular [O III]λ5008 or Hα line fluxes (> 7.0 × 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α emitters (LAEs) from Trainor et al. (2016) with spectroscopic detections of Lyα and [O III]λ5008 or Hα 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 ≈ 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. 2We include 23 of these objects here (Figure 1), omitting the Lyα-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 NIR-Spec/MSA (3 galaxies).The final sample is largely typical of KBSS galaxies, with ⟨z⟩ = 2.4, masses spanning log(M * /M ⊙ ) = 8.5 − 10.7 and a median value of log(M * /M ⊙ ) = 9.7 (assuming a Chabrier 2003 stellar initial mass function).Based on Hα and Hβ measurements from ancillary MOSFIRE spectra, the included galaxies have SFRs ranging from 16−42 M ⊙ yr −1 , with a median SFR Hα = 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).

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 ∼ 2 − 3 galaxies.We used (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]λ4364 in a representative sample of highz galaxies would be prohibitively expensive; the 3σ limiting line flux of 6 × 10 −19 erg s −1 cm −2 achievable in a combined 29 hr G140M exposure (red line) probes ≲ 30% of the z ∼ 2 − 3 sample from Strom et al. (2018).Fortunately, the typical predicted line fluxes for the sum of the [O II]λλ7322, 7332 lines (purple hatched region) and [S III]λ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σ limiting line flux of 4.1 × 10 −19 erg s −1 cm −2 reaches ∼ 90% of typical galaxies.
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 H , stellar Fe/H, and ionization parameter U ) to the properties of z ∼ 2 − 3 KBSS galaxies reported by Strom et al. (2018), which are consistent with the values reported for other z ∼ 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 ≳ 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]λ4364 with increasing O/H and implies a limited ability to detect typical z ∼ 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σ line flux sensitivity of 4.1×10 −19 erg s −1 cm −2 in G235M (corresponding to a total ≈ 30 hr exposure time using the pre-flight ETC; see Section 2.3.1)enables the detection of [S III]λ6313 and [O II]λλ7322, 7332 at virtually all U , even in galaxies with relatively high gas-phase O/H.It is comparatively easier to detect [S III]λ6313 and [O II]λλ7322, 7332 not only because they are predicted to be intrinsically brighter than [O III]λ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 .

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.

Exposure time requirements
NIRSpec G235M observations of [S III]λ6313 and [O II]λλ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. ′′ 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σ limiting line flux of 4.1×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 × 6 integrations × 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.

MSA design
The MSA configuration is central to the success of CE-CILIA, 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. 3For 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 coverage4 (using MSAViz5 ) 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. ′′ 025 − 0. ′′ 01 position steps to optimize the PA and final mask.We did not reach convergence,6 even with angle and position steps much finer than suggested by JDox, suggesting that significant computational resources would be required to fully optimize NIR-Spec/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.

JWST/NIRSpec Data Reduction
The uncalibrated raw G235M data (uncal) frames were processed using the jwst level1 pipeline in the grizli7 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 msaexp 8 version 0.6.11from 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.

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 oversubtraction poses a particular issue at the wavelengths of bright emission lines, which frequently extend well beyond the typical 0. ′′ 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 addi-tive offset9 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 Tool10 (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,11 including uncertain variations in the spectral response as a function of slit position.Likewise, the unexplained addi-tive 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.

THE CHARACTERISTIC REST-OPTICAL
SPECTRUM OF ⟨Z⟩ ∼ 2.4 GALAXIES CECILIA contains some of the deepest spectra obtained during Cycle 1, with ∼ 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 (λ rest ≈ 5700 − 8500 Å) spectrum of a z ∼ 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.

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 (≥ 2× 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 be- tween λ 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]λ6313 and [O II]λλ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 λ 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 ≳ 15, where we define the S/N as the ratio of the composite spectrum to half the 68% CI.This requirement results in ≳ 75% (≥ 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 ∼ 690 object-hours of exposure time.
Aside from the "strong" Hα, [N II]λλ6550, 6585, and [S II]λλ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 ∼ 2 galaxies.Lines longward of ∼ 7000 Å are virtually inaccessible from the ground at z ≳ 2, due to a combination of the rising thermal background in K-band and decreasing atmospheric transparency.In more recent studies of highz 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.

The Nebular Composite Spectrum
To quantify the strength of these emission lines, we construct a second, continuum-subtracted composite spectrum.In this case, after each spectrum is fluxcorrected 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 (∆λ rest ∼ 200 Å) to avoid over-correcting near the emission lines.The spectra are then normalized by the measured flux in Hα and median combined.Because Hα 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 (λF λ ) and rescaled so that the peak flux of (narrow) Hα 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α 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]λλ6550, 6585 and semistrong [O I]λλ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α, [N II]λλ6550, 6585, and [S II]λλ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α + [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σ.Nineteen emission lines satisfy these criteria and are listed in Table 1.We also include Si II λ6373 (2.9σ), the weaker [Ni II] line at 7414 Å (2.7σ), the Paschen line at 8470 Å (3.1σ, 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α 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α.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.

FAINT EMISSION LINES IN HIGH-REDSHIFT STAR-FORMING GALAXIES
In this section, we highlight individual semi-strong (≈ 2 − 3% of Hα) and faint (≲ 1% of Hα) emission lines detected in the CECILIA composite spectra and briefly comment on how they may be used to study highz galaxies.

Auroral Lines
Of all the faint lines present in the rest-optical spectra of star-forming regions, the auroral12 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]λ4364, which falls at λ obs ≈ 1.3 − 1.7 µm for the z ∼ 2−3 CECILIA galaxies.As described in Section 2.2, CECILIA instead targets two auroral lines at longer  1.The residuals from the model are shown in the bottom panels (medium blue), compared to the uncertainties on the median stack (light blue).wavelengths that are not only predicted to be stronger than [O III]λ4364, but also fall at λ obs ≳ 2.0 µ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: The strongest of these is [O II]λλ7322, 7332, with both lines observed to be ∼ 1% the strength of the narrow component of Hα (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 ∼ 15 hr using Keck/MOSFIRE.Using their measurements to calculate direct-method oxygen abundances, they find moderate 12 + log(O/H) = 7.89 ± 0.20 and 12 + log(O/H) = 8.24 ± 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 CE-CILIA is [S III]λ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 CE-CILIA 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α) and proximity to the comparatively stronger [O I]λ6302 line, [S III]λ6313 is unlikely to be accessible in shallower or lower resolution JWST observations of objects similar to the CECILIA sample.
The third auroral line observed in the composite spectra is [N II]λ5756, at 0.15 − 0.25% the strength of Hα (left panel of Figure 5).It is formally detected at 4.6σ 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]λ5756 to calculate T e and the corresponding direct-method ionic abundance.First, we use the line strengths for [S II]λλ6718, 6733 to determine the electron density and find n e ≈ 285 cm −3 , which is consistent with values previously reported for KBSS galaxies (Strom et al. 2017) and other z ∼ 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α and Hβ, 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 ± 2540 K, where the reported uncertainties also capture the likely range in reddening for z ∼ 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]λλ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.

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]λ6302 (upper left panel of Figure 6), which is observed to be 2.36 − 3.08% the strength of the narrow component of Hα and is significantly stronger than the auroral [S III]λ6313 line in its red wing.Its partner line at 6365 Å is ≈ 3.15× weaker as set by atomic physics and is blended with Si II λ6373, which also probes mostly neutral and lowionization gas.Rather than being an abundance diagnostic, [O I]λ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]λ6302 is not commonly reported for individual z ≳ 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 widelyspaced [Ar III]λλ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]λ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α, comparable to [O I]λ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]λ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 ∼ 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 λ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β fluorescence.The strength of O I λ8449 in the CECILIA composite relative to Hα (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 λ8449 as a metallicity tracer; García-Rojas et al. (2006)   Also puzzling is the detection of [Ni II]λλ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α, 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]λ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]λ7380 was seen to be much stronger than expected relative to the associated [Ni II]λ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 ≈ 10 6 cm −3 ) gas (Bautista et al. 1996).Similar to O I λ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.

Broad Hα
Broad line emission has been observed in both spectra of individual galaxies and composite spectra of galaxies at z ∼ 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α line is 6.01 − 28.31% the strength of the narrow component of Hα, 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 ∼ 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 ∼ 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 ∼ 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.

CONCLUSIONS
We have reported the first results from the CECILIA program (JWST PID 2593), which obtained ultra-deep ∼ 30 hr NIRSpec/G235M observations of 33 starforming galaxies at z ∼ 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 restoptical 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α, [N II]λλ6550, 6585, and [S II]λλ6718, 6733.The strengths of these lines relative to the narrow component of Hα are reported in Table 1.
• Aside from strong [N II], Hα, and [S II], which have previously been studied in large groundbased spectroscopic samples, the majority of emission lines are ≲ 3% the strength of Hα.Some of these features, such as [O I]λ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 ≳ 2 − 3% that of Hα 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 λ5877 and [Ar III]λ7138 (shown in the bottom left panel of Figure 6).
• The three auroral emission lines present at λ rest ≈ 5700 − 8500 ([N II]λ5756, [S III]λ6313, [O II]λλ7322, 7332, shown in Figure 5) are ≲ 1% the strength of Hα.Using our measurements of auroral and nebular [N II], we find T e [N II] = 13630 ± 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α 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 ≳ 2 galaxies, evidenced not only by what we have presented in this let, but also by the many exciting results based on NIR-Spec/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 observationswith 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.This work is primarily based on observations made with NASA/ESA/CSA JWST, associated with PID 2593, which can be accessed via doi: 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.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)

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

Figure 2 .
Figure 2. The hatched regions show the model predictions for auroral lines in the rest-optical spectra of typical z ∼ 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]λ4364 (blue hatched region) are ∼ 1 dex fainter than the depth of typical ground-based spectra of individual galaxies, represented by the distribution of 3σ upper limits on [O III]λ4364 from KBSS (dashed black histogram).Sanders et al. (2020) reported four ground-based detections of unlensed [O III]λ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]λ4364 in a representative sample of highz galaxies would be prohibitively expensive; the 3σ limiting line flux of 6 × 10 −19 erg s −1 cm −2 achievable in a combined 29 hr G140M exposure (red line) probes ≲ 30% of the z ∼ 2 − 3 sample fromStrom et al. (2018).Fortunately, the typical predicted line fluxes for the sum of the [O II]λλ7322, 7332 lines (purple hatched region) and [S III]λ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σ limiting line flux of 4.1 × 10 −19 erg s −1 cm −2 reaches ∼ 90% of typical galaxies.

Figure 3 .
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α, [N II]λλ6550, 6585, and [S II]λλ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.

Figure 4 .
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α are reported in Table1.The residuals from the model are shown in the bottom panels (medium blue), compared to the uncertainties on the median stack (light blue).

Figure 5 .
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]λ5756 (left panel) is the weakest of the three and is estimated to be 0.15 − 0.25% the strength of Hα.Both [S III]λ6313 (center panel) and [O II]λ7322, 7332 (right panel) are noticeably stronger, and the oxygen lines are clearly the brightest auroral features in the λ = 5700 − 8500 Å range, with each being ∼ 1% the strength of Hα..

Figure 6 .
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]λ6302, 6565 (upper left), with [S III]λ6313 and Si II λ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]λλ7138, 7753 (lower left) and [Ni II]λ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 λ8449 line and forbidden [Ni II] are only rarely observed, even in observations of nearby galaxies and H II regions.

Figure 7 .
Figure 7.The same nebular composite spectrum from Figure 4 is shown near the strong Hα, [N II]λλ6550, 6585, and [S II]λλ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α and [N II] (in red).The narrow components of all three lines are illustrated by the dot-dashed navy curves.
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).

Table 1 .
Observed Line Fluxes Relative to Hα.
and García-Rojas et al. (2007) used both O I λ7771 and O I λ8449 to derive O + abundances in Galactic H II regions but found an order of magnitude larger O + /H + using O I λ8449, likely due to this effect.