We present results on the emissionline properties of z1.47.5 starforming galaxies in the Assembly of Ultradeep Restoptical Observations Revealing Astrophysics AURORA Cycle 1 JWSTNIRSpec program. Based on its depth continuous wavelength coverage from 15 microns and medium spectral resolution R1000 AURORA includes detections of a large suite of nebular emission lines spanning a broad range in rest wavelength. We investigate the locations of AURORA galaxies in multiple different emissionline diagrams including traditional BPT diagrams of OIIIHbeta vs. NIIHalpha SIIHalpha and OIHalpha and the ionizationmetallicity diagram of OIIIOII O32 vs. OIIIOIIHbeta R23. We also consider a bluer restframe ionizationmetallicity diagram introduced recently to characterize zgt10 galaxies NeIIIOII vs. NeIIIOIIHdelta as well as longerwavelength diagnostic diagrams extending into the restframe nearIR OIIIHbeta vs. SIIISII S32 and HeIPagamma and SIIIPagamma vs. FeIIPabeta. With a significant boost in signaltonoise and large representative samples of individual galaxy detections the AURORA emissionline diagrams presented here definitively confirm a physical picture in which chemicallyyoung alphaenhanced massive stars photoionize the ISM in distant galaxies with a harder ionizing spectrum at fixed nebular metallicity than in their z0 counterparts. We also uncover previously unseen evolution prior to z2 in the OIIIHbeta vs. NIIHalpha diagram which motivates deep NIRSpec observations at even higher redshift. Finally we present the first statistical sample of restframe nearIR emissionline diagnostics in starforming galaxies at high redshift. In order to truly interpret restframe nearIR line ratios including FeII we must obtain better constraints on dust depletion in the highredshift ISM.
['Astrophysics - Astrophysics of Galaxies']
The AURORA Survey A New Era of Emissionline Diagrams with JWSTNIRSpec
{'The AURORA Survey: A New Era of Emission-line Diagrams with JWST /NIRSpec': 'ALICE E. SHAPLEY, 1 RYAN L. SANDERS, 2 MICHAEL W. TOPPING, 3 NAVEEN A. REDDY, 4 DANIELLE A. BERG, 5 RYCHARD J. BOUWENS, 6 GABRIEL BRAMMER, 7 ADAM C. CARNALL, 8 FERGUS CULLEN, 8 ROMEEL DAVÉ, 8 JAMES S. DUNLOP, 8 RICHARD S. ELLIS, 9 N. M. FÖRSTER SCHREIBER, 10 STEVEN R. FURLANETTO, 1 KARL GLAZEBROOK, 11 GARTH D. ILLINGWORTH, 12 TUCKER JONES, 13 MARISKA KRIEK, 6 DEREK J. MCLEOD, 8 ROSS J. MCLURE, 8 DESIKA NARAYANAN, 14 PASCAL OESCH, 15 ANTHONY J. PAHL, 16 MAX PETTINI, 17 DANIEL SCHAERER, 15 DANIEL P. STARK, 3 CHARLES C. STEIDEL, 18 MENGTAO TANG, 3 LEONARDO CLARKE, 1 CALLUM T. DONNAN, 8 AND EMILY KEHOE 1  \n1 Department of Physics & Astronomy, University of California, Los Angeles, 430 Portola Plaza, Los Angeles, CA 90095, USA 2 Department of Physics and Astronomy, University of Kentucky, 505 Rose Street, Lexington, KY 40506, USA 3 Steward Observatory, University of Arizona, 933 N Cherry Avenue, Tucson, AZ 85721, USA 4 Department of Physics & Astronomy, University of California, Riverside, 900 University Avenue, Riverside, CA 92521, USA 5 Department of Astronomy, The University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX 78712, USA 6 Leiden Observatory, Leiden University, NL-2300 RA Leiden, Netherlands 7 Niels Bohr Institute, University of Copenhagen, Lyngbyvej 2, DK2100 Copenhagen Ø, Denmark 8 Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh, EH9 3HJ, UK 9 Department of Physics & Astronomy, University College London. Gower St., London WC1E 6BT, UK 10 Max-Planck-Institut für extraterrestrische Physik (MPE), Giessenbachstr.1, D-85748 Garching, Germany 11 Centre for Astrophysics and Supercomputing, Swinburne University of Technology, P.O. Box 218, Hawthorn, VIC 3122, Australia 12 Department of Astronomy and Astrophysics, UCO/Lick Observatory, University of California, Santa Cruz, CA 95064, USA 13 Department of Physics and Astronomy, University of California Davis, 1 Shields Avenue, Davis, CA 95616, USA 14 Department of Astronomy, University of Florida, 211 Bryant Space Sciences Center, Gainesville, FL, USA 15 Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland 16 The Observatories of the Carnegie Institution for Science, 813 Santa Barbara Street, Pasadena, CA 91101, USA 17 Institute of Astronomy, Madingley Road, Cambridge CB3 OHA, UK 18 Cahill Center for Astronomy and Astrophysics, California Institute of Technology, MS 249-17, Pasadena, CA 91125, USA', 'ABSTRACT': 'Wepresent results on the emission-line properties of z = 1 . 4 -7 . 5 star-forming galaxies in the Assembly of Ultradeep Rest-optical Observations Revealing Astrophysics (AURORA) Cycle 1 JWST /NIRSpec program. Based on its depth, continuous wavelength coverage from 1-5 µ m, and medium spectral resolution ( R ∼ 1000), AURORA includes detections of a large suite of nebular emission lines spanning a broad range in rest wavelength. We investigate the locations of AURORA galaxies in multiple different emission-line diagrams, including traditional \'BPT" diagrams of [OIII] λ 5007/H β vs. [NII] λ 6583/H α , [SII] λλ 6717 , 6731/H α , and [OI] λ 6300/H α , and the \'ionization-metallicity" diagram of [OIII] λ 5007/[OII] λ 3727 ( O 32) vs. ([OIII] λ 5007+[OII] λ 3727)/H β ( R 23). We also consider a bluer rest-frame \'ionization-metallicity" diagram introduced recently to characterize z > 10 galaxies: [NeIII] λ 3869/[OII] λ 3727 vs. ([NeIII] λ 3869+[OII] λ 3727)/H δ ; as well as longer-wavelength diagnostic diagrams extending into the rest-frame near-IR: [OIII] λ 5007/H β vs. [SIII] λλ 9069 , 9532/[SII] λλ 6717 , 6731 ( S 32); and He I λ 1 . 083 µ m/Pa γ and [SIII] λ 9532/Pa γ vs. [FeII] λ 1 . 257 µ m/Pa β . With a significant boost in signalto-noise and large, representative samples of individual galaxy detections, the AURORA emission-line diagrams presented here definitively confirm a physical picture in which chemically-young, α -enhanced, massive stars photoionize the ISM in distant galaxies with a harder ionizing spectrum at fixed nebular metallicity than in their z ∼ 0 counterparts. We also uncover previously unseen evolution prior to z ∼ 2 in the [OIII] λ 5007/H β vs. [NII] λ 6583/H α diagram, which motivates deep NIRSpec observations at even higher redshift. Finally, we present the first statistical sample of rest-frame near-IR emission-line diagnostics in star-forming galaxies at high redshift. In order to truly interpret rest-frame near-IR line ratios including [FeII] λ 1 . 257 µ m, we must obtain better constraints on dust depletion in the high-redshift ISM.', '1. INTRODUCTION': 'Recombination and collisionally-excited emission lines from nearby star-forming regions and galaxies provide the key to understanding the interplay between massive stars, gas, dust, and heavy elements in the ionized interstellar medium. In the Orion Nebula, for example, several hundred optical emission lines have been measured (Esteban et al. 2004), enabling an exquisitely detailed view of the density and temperature structure in the gas, as well as the mixture of chemical elements. Similar analyses can be extended to extragalactic H II regions and low-redshift star-forming galaxies (e.g., Izotov et al. 2006; Berg et al. 2015; MéndezDelgado et al. 2023).  \nWith the advent of near-IR spectrographs on 8-10 meter telescopes on the ground and the HST /WFC/IR grism in space, it was possible to extend rest-optical nebular emissionline measurements to samples of star-forming galaxies at z ∼ 1 -3 (e.g., Steidel et al. 2014; Shapley et al. 2015; Sanders et al. 2016; Momcheva et al. 2016; Kashino et al. 2017; Förster Schreiber et al. 2019; Sanders et al. 2021). In such studies, the set of detected emission lines has typically been limited to the brightest rest-optical lines, including [OII] λ 3727, H β , [OIII] λλ 5007 , 4959, H α , [NII] λ 6583, and [SII] λλ 6717 , 6731 (but see, e.g., Sanders et al. 2023b; Clarke et al. 2023; Jeong et al. 2020; Backhaus et al. 2022). Furthermore, the strongly wavelength-dependent and bright sky background due to the forest of OH vibration-rotation transitions at 1 -2 . 5 µ mleads to significant challenges in obtaining robust measurements from the ground for all but the very brightest of the rest-optical emission lines (e.g., H α , [OIII] λ 5007).  \nDespite these observational hurdles, significant progress was achieved over the past decade based on multi-object near-IR spectroscopic surveys out to z ∼ 3. Relations among metallicity, mass, star-formation rate (SFR), dynamical properties, and dust attenuation have been constructed for starforming galaxies, demonstrating the evolution at fixed stellar mass towards higher SFR, lower metallicity, and lower ratio of rotational-to-random dynamical support, with increasing redshift (Shapley et al. 2015; Shivaei et al. 2015; Kashino et al. 2017; Wisnioski et al. 2019; Förster Schreiber & Wuyts 2020; Shapley et al. 2022). Furthermore, galaxies seem to follow the same manifold of stellar mass, SFR, and metallicity at z = 0 -3, suggesting the same equilibrium is in place among inflows, outflows, and star formation (Sanders et al. 2021). Finally, an essential characteristic of distant starforming galaxies emerging from these spectroscopic studies is their chemical youth, encapsulated by evidence for non-solar, α -enhanced chemical abundance patterns. This evidence comes in many forms, including a comparison of nebular oxygen and stellar iron abundances (e.g., Steidel et al. 2016; Topping et al. 2020; Cullen et al. 2021; Stan-  \net al. 2024); the loci of high-redshift galaxies in multiple rest-optical emission-line diagnostic diagrams (e.g. Steidel et al. 2014; Shapley et al. 2019; Sanders et al. 2020b; Jeong et al. 2020; Clarke et al. 2023); and photoionization modeling of strong nebular emission-line ratios incorporating direct metallicity estimates (e.g. Sanders et al. 2020a).  \nThe infrared spectroscopic capabilities of JWST represent a transformational advance for the study of ionized gas in star-forming galaxies at high redshift. The redshift frontier for rest-optical emission-line measurements has been extended from z ∼ 3 out to z ∼ 12 (e.g., Shapley et al. 2023a; Sanders et al. 2023a, 2024; Cameron et al. 2023; Nakajima et al. 2023; Roberts-Borsani et al. 2024; Castellano et al. 2024). Furthermore, the unprecedented wavelength baseline for high signal-to-noise (S/N) infrared spectroscopy has highlighted the promise of extending diagnostic nebular studies into the rest-frame near-IR for galaxies observed during the epoch of peak star formation in the universe ( z ∼ 2 -3) (Brinchmann 2023; Calabrò et al. 2023). At the other end of the rest-frame optical, blue rest-optical emission-line diagnostics (e.g., Bunker et al. 2023) have been obtained for the most extreme galaxies at z > 10. While not the focus of the current work, JWST has also enabled unprecedented rest-UV spectroscopic diagnostics of galaxies at z ∼ 4 -12 (e.g., Tang et al. 2023; Williams et al. 2023; Castellano et al. 2024; Topping et al. 2024). Thus far, these new emission-line diagnostic diagrams are populated by extremely small high-redshift samples, with typically modest S/N. Alternatively, data from stacked measurements is plotted, which increases the S/N of faint features at the cost of obscuring the variation within a given dataset (e.g, Sanders et al. 2023a).  \nThe Assembly of Ultradeep Observations Revealing Astrophysics (AURORA) Cycle 1 JWST /NIRSpec program was designed to detect faint auroral lines from ionized oxygen, sulfur, and nitrogen for z > 1 . 4 galaxies, which can be used to infer electron temperatures and, in turn, direct chemical abundances. As such, AURORA consists of long exposures with JWST /NIRSpec using each of the three mediumresolution ( R ∼ 1000) gratings. Because of their superior depth, the AURORA data yield significantly higher S/N measurements than in previous work, including emission-line excitation sequences for individual high-redshift galaxies spanning multiple diagnostic emission-line ratio diagrams (e.g., [OIII] λ 5007/H β vs. [NII] λ 6583/H α ). Along with deep integration times, AURORA was also designed to obtain continuous wavelength coverage from 1 -5 µ m. The AURORA dataset thus yields both rest-frame blue-optical emission-line diagnostics highlighted for interpreting unique z > 10 spectra, and also the first statistical sample of rest-frame near-IR diagnostics for star-forming galaxies at z > 1. We demonstrate how these diagnostics provide novel insights into the  \nchemical abundance patterns and interplay between massive stars and gas in distant star-forming galaxies.  \nIn §2, we describe the AURORA observations, sample, data reduction, flux calibration including slit-loss corrections, and measurements. In §3, we present results on familiar emission-line diagnostic diagrams, extended towards higher redshift, and novel emission-line diagrams spanning a broader wavelength baseline and measured for the first time with actual statistical power. In §4, we consider the implications of these new measurements and discuss future directions. Throughout, we adopt cosmological parameters of H 0 = 70 km s -1 Mpc -1 , Ω m = 0 . 30, and Ω Λ = 0 . 7, and a Chabrier (2003) IMF.', '2.1. AURORA NIRSpec Observations': 'We use NIRSpec Micro Shutter Assembly (MSA) data from the AURORA program (Program ID: 1914; Co-PIs Shapley and Sanders). The AURORA NIRSpec observations we analyzed consisted of 2 pointings, one in the COSMOS field (PA=71.82 degrees) and one in GOODSN (PA=70.34 degrees). Both pointings utilized the grating/filter combination of G140M/F100LP, G235M/F170LP, and G395M/F290LP, which provide a spectral resolution of R ∼ 1000 over the wavelength range approximately 1 -5 µ m. In each pointing, the G140M/F100LP combination was observed for 44204 seconds (30 exposures of 20 groups and the NRSIRS2 readout mode), the G235M/F170LP combination for 28886 seconds (30 exposures of 65 groups and the NRSIRS2RAPID readout mode), and the G395M/F290LP combination for 15056 seconds (6 exposures of 85 groups and the NRSIRS2RAPID readout mode). The readout modes for each grating/filter were chosen to conform to the data volume allocation. Based on the NIRSpec Exposure Time Calculator, the exposure time for each grating was designed to yield uniform ∼ 5 σ emission-line sensitivity of 10 -18 erg s -1 cm -2 . A 3-point nod pattern was adopted for each observation, and each MSA \'slit" consisted of 3 microshutters.', '2.2. AURORA Target Sample': 'The COSMOS and GOODS-N masks included, respectively, 46 and 51 galaxy targets. Once the PA for each MSA pointing was assigned to our program, we used the NIRSpec MSA Planning Tool (MPT) to optimize the number of targets predicted to yield a significant detection for at least one auroral oxygen emission line (i.e., [OIII] λ 4363 and/or [OII] λλ 7320 , 7330). In the considerable time that elapsed between the submission of our Cycle 1 proposal in November 2020 and its execution in late 2023/early 2024, updated empirical information became available regarding the sensitivity of NIRSpec. Scaling from NIRSpec R ∼ 1000 measurements from the CEERS Early Release Science (ERS) program  \n(Finkelstein et al. 2023), we therefore lowered the threshold for \'detection" to 5 × 10 -19 erg s -1 cm -2 (3 σ ). In identifying potential auroral-line targets, we utilized publicly-available multi-wavelength photomety and redshifts from the 3D-HST survey (Skelton et al. 2014), rest-optical emission-line catalogs of the MOSFIRE Deep Evolution Field (MOSDEF) survey (Kriek et al. 2015), and the Extreme Emission-line Galaxy (EELG) catalogs of Tang et al. (2019). We used the set of measured strong rest-frame optical nebular emission lines in the MOSDEF and 3D-HST EELG catalogs, the relationships between strong-line ratios and [OIII] electron temperatures ( T 3 and T 2 for the ionization zones traced, respectively, by [OIII] λ 5007 and [OII] λ 3727) derived in Sanders et al. (2017) and Sanders et al. (2021), and the relationships among the strengths of several auroral lines measured in local H II regions (e.g., Berg et al. 2020) to predict T 3 and T 2, and, accordingly the strengths of the following auroral features: [OIII] λ 4363, [SIII] λ 6312, [OII] λλ 7320 , 7330, [NII] λ 5755, and [SII] λ 4068 , 4076. To these samples, we added galaxies with preJWST spectroscopic measurements of [OIII] λ 5007 flux greater than or equal to 10 -16 erg s -1 cm -2 , or with SEDs indicating a photometric excess at the wavelength of [OIII] λ 5007 corresponding to a line flux greater than or equal to 10 -16 erg s -1 cm -2 . For objects with f ([OIII] λ 5007) ≥ 10 -16 erg s -1 cm -2 , we predicted a detection of [OIII] λ 4363 at ≥ 3 σ if T 3 ≥ 10 , 000K, based on our limiting line sensitivity.  \nThese criteria resulted in 249 (210) potential auroral-line targets in COSMOS (GOODS-N), ranging in redshift from 1 . 2 < z ≤ 4 . 65. These potential targets were prioritized using the following scheme (wherein we note that galaxies at z ≥ 1 . 63 receive higher priority given that [OII] λ 3727 is covered by NIRSpec in this redshift range):  \n- 1. z ≥ 1 . 63 and both [OIII] λ 4363 and [OII] λλ 7320 , 7330 have estimated flux ≥ 7 × 10 -19 erg s -1 cm -2\n- 2. z ≥ 1 . 63 and (both [OIII] λ 4363 and [OII] λλ 7320 , 7330 have estimated flux between 5 × 10 -19 erg s -1 cm -2 and 7 × 10 -19 erg s -1 cm -2 ) or (only one of [OIII] λ 4363 or [OII] λλ 7320 , 7330 have estimated flux ≥ 7 × 10 -19 erg s -1 cm -2 ); or z < 1 . 63 and both [OIII] λ 4363 and [OII] λλ 7320 , 7330 have estimated flux ≥ 7 × 10 -19 erg s -1 cm -2\n- 3. z ≥ 1 . 63 and [OIII] λ 5007 flux ≥ 1 × 10 -16 erg s -1 cm -2 from spectroscopy or estimated to have [OIII] λ 5007 flux ≥ 1 × 10 -16 erg s -1 cm -2 from photometric excess; z < 1 . 63 and (both [OIII] λ 4363 and [OII] λλ 7320 , 7330 have estimated flux between 5 × 10 -19 erg s -1 cm -2 and 7 × 10 -19 erg s -1 cm -2 ) or (only one of [OIII] λ 4363 or [OII] λλ 7320 , 7330 have estimated flux ≥ 7 × 10 -19 erg s -1 cm -2 )  \n- 4. z ≥ 1 . 63 and only one of [OIII] λ 4363 or [OII] λλ 7320 , 7330 have estimated flux between 5 × 10 -19 erg s -1 cm -2 and 7 × 10 -19 erg s -1 cm -2 ; z < 1 . 63 and [OIII] λ 5007 flux ≥ 1 × 10 -16 erg s -1 cm -2 from spectroscopy or estimated to have [OIII] λ 5007 flux ≥ 1 × 10 -16 erg s -1 cm -2 from photometric excess  \nUsing the MPT, we were able to assign slits to 16 (20) auroral-line targets on the COSMOS (GOODS-N) mask, and optimized the mask design such that the lines of interest (auroral features and their strong nebular counterparts) fell on the NIRSpec detectors. These 36 targets span a redshift range of 1 . 39 ≤ z ≤ 4 . 41, with a median of z med = 2 . 33; all but 3 are at z ≥ 1 . 6 and therefore have coverage of [OII] λ 3727.  \nBased on the science goals of the AURORA program, we included additional sources in each of the COSMOS and GOODS-N AURORA masks. In COSMOS, we made use of the photometric redshift catalog from the PRIMER survey (J. Dunlop et al. 2024, in prep.), which is defined down to a F356W limit of m AB , 356W = 28. We added 7 sources at z phot > 6, 1 quiescent galaxy candidate at z ∼ 2, and 22 galaxies at 1 . 5 ≤ z phot ≤ 5 and m AB , 356W ≤ 27. In GOODS-N, we filled the mask with photometric candidates and spectroscopicallyconfirmed galaxies from the FRESCO survey (two H α emitters at 5 ≤ z ≤ 6, one [OIII] emitter at 7 . 2; and 3 quiescent galaxy candidates at 2 ≤ z < 3); 4 z > 6 photometric candidates from the JADES survey (Hainline et al. 2023) and 6 the literature (Jung et al. 2020; Finkelstein et al. 2015; Bouwens et al. 2015); and 15 galaxies at 1 . 5 ≤ z phot ≤ 5 from the publicly available GOODS-N photometric redshift catalog in the Dawn JWST Archive 1 (Heintz et al. 2024). In total we assigned slits to 97 AURORA targets.', '2.3. Data Reduction': "We obtained final two-dimensional (2D) spectra using a data reduction pipeline composed of standard STScI tools in addition to custom routines developed for the AURORA survey. Each raw uncalibrated exposure was processed using the calwebb\\_detector1 routine within the jwst pipeline 2 . This step first removed the detector bias and dark current signal, and corrected for detector gain. We then applied a linearity correction to the images, and derived a count rate for each pixel from the non-destructive readout groups. In calculating the count rates, we mitigated significant jumps in subsequent readout groups, such as those resulting from cosmic rays and including large events identified as 'snowballs.' Next, we corrected each of the resulting count rate files for 1 / f noise using the nsclean package (Rauscher 2024). We extracted the 2D spectrum traced out from each slit from the corrected  \nrate files, and applied a flat-field correction, wavelength solution, and absolute photometric calibration to each spectrum using reference files from the updated Calibration Reference Data System (CRDS) context ( jwst\\_1193.pmap ) based on in-flight estimates. As each grating, detector, and nod position was observed twice with the same pointing center and MSA configuration, we compared the repeated exposures to identify any transient artifacts (e.g., remaining cosmic rays or hot pixels) not excluded by earlier stages of the reduction. The identified pixels were flagged as to not be utilized for the final spectrum. Finally, the cutout spectrum from each slit was rectified, and interpolated on to a common wavelength grid before being combined following the three-dither nodding pattern described above, yielding the final 2D spectrum for each object.  \nOne-dimensional (1D) science and error spectra were optimally extracted from the rectified 2D spectra. During extraction, the spatial profile was based on the brightest emission line in each grating (in 88% of grating spectra), or the integrated continuum profile if no emission lines were detected (in 10.5% of grating spectra). Otherwise, in a very small minority of cases ( ∼ 1%) a blind extraction was performed based on the spatial profiles in the other grating(s). Any artifacts or cosmic rays in the extracted trace not removed during the 2D reduction were manually flagged during the 1D extraction and masked during spectral fitting.", '2.4. Slit Loss Correction': "The combination of a small NIRSpec MSA microshutter width (0 . '' 20), potentially off-center location of the target within the microshutter, and a wavelength-dependent PSF necessitates accurate corrections for light falling outside the microshutter and 1D extraction window. These 'slit losses' were computed as follows (see Reddy et al. 2023 for further details). Using the NIRCam F115W imaging, a 12 '' × 12 '' subimage centered on the target was extracted. The corresponding segmentation map was used to mask all pixels belonging to unrelated objects within the subimage. The subimage was then rotated based on the position angle of the NIRSpec MSA observations. This rotated subimage was then convolved with Gaussian kernels to produce the expected light profiles at longer wavelengths. The full widths at half maximum (FWHMs) of these kernels were computed by subtracting in quadrature the FWHM of the JWST PSF at 1 . 15 µ m from the FHWMs of the JWST PSFs at λ = 1 . 2 -5 . 3 µ m, in increments of ∆ λ = 0 . 1 µ m.  \nFor 4 and 2 objects in COSMOS and GOODS-N, respectively, without NIRCam F115W imaging, the HST/WFC3 F160W image was used to generate a Sérsic fit describing the intrinsic shape of the galaxy. The parameters of the Sérsic fit were obtained from van der Wel et al. (2014), and a model was created only if the fit was deemed acceptable  \nbased on the criteria given in that study. This model was then convolved with the JWST PSFs at λ = 0 . 6 -5 . 3 µ m, in increments of ∆ λ = 0 . 1 µ m, as determined from the JWST WebbPSF software 3 . Finally, a point source was assumed for 5 targets in GOODS-N without F115W imaging or an acceptable Sérsic model fit to their F160W light profiles.  \nThe convolved images of the target at each wavelength were then shifted according to the expected position of the target within the microshutter and then masked. The masking accounts for the 0 . '' 20 microshutter width, the 0 . '' 07 gap between adjacent microshutters along the cross-dispersion axis, and the center and size of the window used to extract the 2D spectrograms to 1D. The 1D spectra were then divided by the fraction of light passing through the slit and extraction apertures as a function of wavelength.", '2.5. Flux Calibration': 'The procedures for flux calibration are described in detail in Sanders et al. 2024 (in prep.). Briefly, the final flux calibration for AURORA spectra proceeded in two stages, addressing first the relative grating-to-grating calibration, and then the overall absolute calibration across all gratings.  \nThe depth of the AURORA observations results in continuum and/or emission lines detected in the overlapping regions for G140M+G235M, G235M+G395M, or both, for a majority of AURORA targets. Overlapping measurements were used to improve the grating-to-grating flux calibration, relative to the output of the NIRSpec pipeline (Section 2.3). The G235M grating was adopted as the reference relative to which we assessed the calibration of the G140M and G395M spectra. The median G140M/G235M and G395M/G235M flux ratios in overlapping wavelength regions were assessed separately for each AURORA mask, and found to be near unity for all grating pairs except G395M/G235M in COSMOS, which was found to be 0.80. With a small number of exceptions, G140M or G395M spectra were scaled by the median G140M/G235M or G395M/G235M flux ratio, respectively, determined from the sample of line and continuum measurements in overlapping spectral regions. Exceptions consisted of outliers, whose grating flux ratios lay outside the > 1 σ scatter of the distribution around the median. For such objects, we forced the grating flux levels to agree in the region of overlap. In order to assess the overall robustness of the grating-to-grating flux calibration, we measured fluxes for 156 emission lines detected at S/N ≥ 3 in the overlap region of neighboring gratings. For these, we found a median flux offset of 0.1% and an intrinsic scatter of 8% after accounting for measurement uncertainties.  \nFigure 1. AURORAredshift distribution. This distribution includes 36 auroral-line targets spanning a range 1 . 4 ≤ z ≤ 4 . 5; 18 photometric or emission-line targets at z ≥ 5; 4 quiescent-galaxy targets at 2 ≤ z ≤ 3; and 37 filler targets at 1 . 5 ≤ z ≤ 5 based on photometric redshifts.  \n<!-- image -->  \nThe absolute flux calibration for AURORA is obtained by scaling the spectra to match the available multi-band photometry. The spectrum of each target in each grating was passed through JWST /NIRCam HST /WFC3 filter transmission curves that fully overlapped the grating of interest, to obtain synthetic photometric measurements. The ratio of photometric flux measurements from imaging and NIRSpec spectra was calculated for all filters in which both types of measurements had S/N ≥ 3. We applied a single scale factor to the NIRSpec spectra to force the median of the spectrumto-imaging flux ratios to be unity. The median scale factor was 1.35, with a standard deviation of 0.17 dex.', '2.6. Emission Line Measurements': 'We measured redshifts and emission-line fluxes for 95 out of 97 galaxies targeted in the AURORA sample. The distribution of redshifts measured from NIRSpec spectra is shown in Figure 1. We used single Gaussian fits for widelyseparated lines; adjacent lines such as [NII] λ 6548, H α , and [NII] λ 6583 were fit simultaneously with multiple Gaussians; and closely spaced lines that are blended and unresolved at R ∼ 1000 are fit with a single Gaussian. As a first pass, the underlying continuum was taken to be the best-fit SED model to the raw multi-wavelength photometry. We then used these initial emission-line fluxes to correct the photometry for the contributions from nebular emission, and repeated the SED fitting. We then repeated the process of emission-line mea-  \nFigure 2. Example AURORA spectrum, for GOODSN-30564 ( z = 2 . 4828). This galaxy was targeted based on the prediction that auroral emission lines ([OIII] λ 4363 and [OII] λλ 7320 , 7330) would yield significant detections in its AURORA spectrum, and both features are successfully detected. Flux density is plotted in black, while the error spectrum is plotted in red. Black labels indicate emission lines analyzed in Section 3, while grey labels correspond to additional visible emission lines whose analysis we save for future work. The spectral windows plotted here are zoomed in to highlight features discussed in this work, but represent only a subset of the rest wavelengths covered for GOODSN-30564 by our NIRSpec set up.  \n<!-- image -->  \nFigure 3. SFR(H α ) vs. M ∗ . Red, green, blue, and magenta symbols are used, respectively, for the 1 . 4 ≤ z < 2 . 7, 2 . 7 ≤ z < 4 . 0, 4 . 0 ≤ z < 5 . 0, and z ≥ 5 . 0 samples. SFR(H α ) was estimated from the inferred dust-corrected H α luminosity, using the conversion factor from Hao et al. (2011). Galaxies flagged as AGNs or quiescent are not included here.  \n<!-- image -->  \nsurements as in the first pass, but using the updated SED fit as the estimate of the continuum. This same best-fit SED model is used for stellar mass estimates, as described in the next section. The second pass of spectral fitting represents our final set of emission-line measurements, and (as in the first pass) automatically accounts for any underlying stellar absorption affecting nebular Balmer emission lines. Figure 2 shows the final flux-calibrated spectrum of one of the auroralline targets from the AURORA sample, GOODSN-30564 ( z = 2 . 4828), which boasts the detection of more than 50 individual emission lines (only a subset of which are shown within the wavelength ranges of Figure 2).', '2.7. SED Fitting': 'The spectral energy distribution (SED) of each galaxy was fit using the a large suite of publicly available photometry from JWST/NIRCam, and HST/WFC3 and ACS, with photometric catalogs downloaded from the DAWN JWST Archive. The COSMOS catalog includes the following photometric bands: (from HST /ACS and WFC3) F435W, F606W, F814W, F850LP, F105W, F125W, F140W, F160W, (and from JWST /NIRCam) F090W, F115W, F150W, F200W, F277W, F356W, F410M, and F444W. The GOODS-N catalog includes the following bands: (from HST /ACS and WFC3) F435W, F606W, F775W, F814W, F850LP, F105W, F125W, F140W, F160W, (and from JWST /NIRCam) F090W, F115W,  \nF150W, F182M, F200W, F210M, F277W, F335M, F356W, F410M, and F444W. In COSMOS, the footprint of our mask almost entirely overlaps the PRIMER imaging area, but 4 sources do not have PRIMER coverage. For such sources, we used photometry from the 3DHST catalog, including groundbased imaging, HST/WFC3 and ACS, and Spitzer/IRAC (Skelton et al. 2014). Furthermore, for one source in COSMOS, the photometry in the DJA catalog is affected by bad pixels at wavelengths longer than 2 µ m. Accordingly, we used measurements from the catalog described in Donnan et al. (2024) (which are consistent with the DJA catalog at wavelengths that overlap), including the HST + JWST bands F435W, F606W, F814W, F090W, F115W, F150W, F200W, F277W, F356W, F410M, and F444W.  \nFor SED modeling, we used the program FAST (Kriek et al. 2009), assuming the flexible stellar population synthesis (FSPS) models of Conroy et al. (2009) and a Chabrier (2003) IMF, and fixed the galaxy redshift to the value determined from AURORA NIRSpec spectra. Following the methodology of Reddy et al. (2018), we utilized two combinations of metallicity and extinction/attenuation curves for SED modeling. These include 1.4 solar metallicity ( Z ⊙ = 0 . 014) coupled with the Calzetti et al. (2000) attenuation curve (hereafter \'1.4 Z ⊙ + Calzetti"), and 0.27 solar models accompanied by the SMC extinction curve of Gordon et al. (2003) (hereafter \'0.27 Z ⊙ +SMC"). We also assumed star-formation histories of the form SFR ( t ) ∝ t × exp ( -t /τ ) (i.e., delayedτ models). Here, t is time since the onset of star formation and τ is the characteristic star-formation timescale. The choice of 1.4 Z ⊙ + Calzetti or 0.27 Z ⊙ +SMC was assigned to each galaxy based on the combination that yielded the lower χ 2 . Accordingly, 76 galaxies were modeled assuming 1.4 Z ⊙ + Calzetti, and 18 galaxies with 0.27 Z ⊙ +SMC. One filler target in GOODS-N at z = 1 . 81 does not have a robust SED.  \nFinally, all relevant photometric bands were corrected for the contributions from strong nebular emission lines using the method described in Sanders et al. (2021). Photometry was additionally corrected for contributions from nebular continuum emission based on predictions from grids of Cloudy photoionization models (Ferland et al. 2017), tied to the measured H β line flux (see Sanders et al. 2024, in prep., for full details).  \nWe compared the stellar masses inferred from the SED modelling described above to those obtained by fitting the uncorrected broadband photometry using PROSPECTOR (Johnson et al. 2021), assuming a non-parametric starformation history (SFH), and incorporating nebular emission in the SED models (see Topping et al. 2024, in prep., for a full description). Briefly, for this model setup, we assumed a SFH comprising eight independent time bins constrained by a continuity prior (see, e.g., Tacchella et al. 2022). In addition, we fixed the nebular and stellar metallicities to be  \nidentical, and allowed the ionization parameter to vary from log( U ) = -4 to -1, with a uniform prior. The stellar masses inferred from these non-parametric fits that consider nebular emission are in good agreement with those from the stellaronly FAST models, such that PROSPECTOR yielded 0.09 dex larger masses on average, with a scatter between the masses from two methods of 0.45 dex for the full sample.  \nTo place the AURORA sample in context, we plot the SFRs and stellar masses of our targets in Figure 3. Stellar masses are derived as described above. In this figure, SFR is estimated from dust-corrected Balmer-emission-line fluxes. For the majority of galaxies, the H α /H β ratio was used to calculate E ( B -V )gas, assuming the Cardelli et al. (1989) dust extinction curve and an intrinsic ratio of H α /H β = 2 . 79. In a small number of cases where H α or H β was not measured, we used, respectively, the ratio H β /H γ (assuming an intrinsic ratio of H β /H γ = 2 . 11) or H α /H γ (assuming an intrinsic ratio of H α /H γ = 5 . 90) to estimate E ( B -V )gas. The value of E ( B -V )gas was used to infer the dust-corrected H α luminosity, and the conversion from Hao et al. (2011), adjusted to the Chabrier IMF (SFR( M ⊙ yr -1 ) = 10 -41 . 33 L (H α )(ergs -1 )), was used to translate H α luminosity to instantaneous SFR.', '2.8. AURORA Emission-line Sample': 'We investigate several emission-line-ratio diagrams in this work. Each diagram is based on features spanning a different range of rest-frame wavelength. Given the wide redshift range of our targets and the wavelength gaps that are present in the spectra of some AURORA galaxies due to the manner in which the spectra fell across the two NIRSpec detectors, the subset of galaxies plotted in each diagram with coverage of all required lines differs slightly. Furthermore, as the focus of this work is the properties of star-forming galaxies, we separately identify both AGNs and quiescent galaxies within the AURORA sample. As in previous work, AGNs are identified spectroscopically on the basis of broad Balmerline emission and/or [NII]/H α> 0 . 5 (Shapley et al. 2023a) (5 galaxies). Quiescent galaxies are identified either based on their photometric SEDs (4 galaxies) or else spectroscopically on the basis of strong stellar absorption lines and a lack of strong emission lines arising from star formation (1 galaxy). Two objects are identified as both AGNs and quiescent. Accordingly, 87 AURORA galaxies remain in our star-forming sample, and appear in subsequent plots.', '2.9. z ∼ 0 Comparison Samples': "As in past work, in order to understand the evolving emission-line ratios of distant galaxies, we selected comparison samples of both z ∼ 0 galaxies and H II regions. The z ∼ 0 galaxy sample is drawn from the Sloan Digital Sky Survey (SDSS) Data release 7 (DR7; Abazajian et al. 2009). Stellar-absorption-corrected emission-line measure-  \ns are drawn from the MPA-JHU catalog of measurements for DR7 4 . We limited the SDSS sample to galaxies at 0 . 04 ≤ z ≤ 0 . 10 to reduce aperture effects. SDSS does not have coverage of the [SIII] λλ 9069 , 9532 doublet, which we measure for AURORA galaxies at z ≤ 4 . 5. Therefore, when investigating the relationship between [OIII] λ 5007 and [SIII]/[SII], we also fold in data from the Mapping Nearby Galaxies at APO (MaNGA) integral field spectroscopic survey (Bundy et al. 2015), which extends out to 1 . 04 µ m. We include the properties of 1150 MaNGA targets at z < 0 . 08, such that [SIII] is covered, and sum a subset of pixels in the central 3 '' for each galaxy to simulate a SDSS fiber spectrum (see Sanders et al. 2020b, for details).  \nFor some of the emission-line ratio diagrams presented here - in particular, those that include [OI] and [SII] emission - galaxies and H II regions occupy different regions of parameter space (Zhang et al. 2017). Such differences arise because the SDSS emission-line spectra of z ∼ 0 galaxies include the contributions of both H II regions and diffuse ionized gas (DIG). As demonstrated in Sanders et al. (2017) and Shapley et al. (2019), due to their ∼ 2 orders of magnitude higher surface density of star formation ( Σ SFR), and the anticorrelation between Σ SFR and DIG fractional emission-line contribution to the integrated spectrum (Oey et al. 2007), the integrated spectra of z > 2 galaxies likely contain minimal contributions from DIG. Therefore, it is essential to compare the emission-line ratios of z > 2 galaxies with those of z ∼ 0 H II regions. For this comparison, we assemble measurements for a sample of 298 H II regions from the CHemical Abundances Of Spirals (CHAOS) survey (Berg et al. 2015, 2020; Croxall et al. 2015, 2016; Rogers et al. 2021).", '3. RESULTS': 'JWST /NIRSpec boasts a transformative advance in both sensitivity and wavelength coverage. These developments facilitate new insights into \'classical" emission-line diagrams such as [OIII] λ 5007/H β vs. [NII] λ 6583/H α , [OIII] λ 5007/H β vs. [SII] λλ 6717 , 6731/H α , [OIII] λ 5007/H β vs. [OI] λ 6300/H α , and [OIII] λ 5007/[OII] λ 3727 (O32) vs. ([OII] λ 3727+[OIII] λλ 4959 , 5007)/H β (R23). Notably, for the first time we construct these diagrams including ≥ 3 σ detections for the vast majority (i.e., ≥ 80%) of the individual galaxies in the AURORA sample with coverage of the relevant lines. In past work (e.g., Shapley et al. 2015; Strom et al. 2017; Cameron et al. 2023), the corresponding all-line detection fraction was significantly lower. For example, in the CEERS NIRSpec analysis of Sanders et al. (2023a), it was ≤ 50% in diagrams involving [NII] or [SII], and ≤ 15% when [OI] was being considered. In the face of such low  \nFigure 4. \'Classical" Emission-line Diagnostic Diagrams. In each panel, galaxies from the AURORA sample are plotted with colored symbols. Red, green, blue, and magenta are used, respectively, for the 1 . 4 ≤ z < 2 . 7, 2 . 7 ≤ z < 4 . 0, 4 . 0 ≤ z < 5 . 0, and z ≥ 5 . 0 samples. Star-forming galaxies with ≥ 3 σ detections of all 4 BPT emission lines are plotted with circles. Limits (3 σ ) are plotted as arrows. The 2D greyscale histogram and cyan dashed curve correspond, respectively, to the distribution of z ∼ 0 SDSS galaxies and CHAOS H II regions. Top left: [NII] BPT diagram. The black dotted curve is the \'maximum starburst" line from Kewley et al. (2001), while the black solid curve is an empirical AGN/star-formation threshold from Kauffmann et al. (2003). Top Right: [SII] BPT diagram. The black dotted curve is the \'maximum starburst" line from Kewley et al. (2001). Bottom: [OI] BPT diagram. The black dotted curve is the \'maximum starburst" line from Kewley et al. (2001).  \n<!-- image -->  \ndetection fractions, composite spectra were used to obtain an unbiased estimate of the average sample properties. Because of the high detection fraction in AURORA, there is no need to rely on composite spectra.  \nIn addition, it is possible to explore new parameter spaces, including emission lines at bluer rest-frame optical wavelengths that are useful for characterizing the highest-redshift galaxies (e.g., Bunker et al. 2023). Such features (e.g., [NeIII] λ 3869, H δ ) did not typically yield sufficient S/N for analysis in high-redshift galaxy spectra, preJWST . Finally, the extension of wavelength coverage out to 5 µ m enables a systematic investigation of emission lines at redder rest-frame wavelengths for galaxies at z ≤ 3, such as [SIII] λλ 9069 , 9532 (Sanders et al. 2020b), He I λ 1 . 083 µ m, [FeII] λ 1 . 257 µ m, and the Paschen series (Brinchmann 2023; Calabrò et al. 2023). We begin this section by revisiting \'classical" emission-line diagrams, and then turn to more novel parameter spaces at both shorter and longer restframe wavelengths. For closely-spaced emission-line ratios ([OIII] λ 5007/H β , [NII] λ 6583/H α , [SII] λλ 6717 , 6731/H α , [OI] λ 6300/H α ) and those in the rest-frame near-IR, we do not apply a dust correction. However, for O 32, R 23, [NeIII] λ 3869/[OII] λ 3727, ([NeIII] λ 3869+[OII] λ 3727)/H δ , and S 32, we use dust-corrected line ratios based on E ( B -V )gas, inferred as described in Section 2.7.', '3.1. 3.1.1. BPT Diagrams': 'We begin by revisiting familiar emission-line diagnostic diagrams, introduced by Baldwin et al. (1981) (hereafter BPT) and Veilleux & Osterbrock (1987) to classify extragalactic emission-line spectra as gas photoionized by young or evolved stars, shocks, or an AGN. Figure 4 includes plots of [OIII] λ 5007/H β vs. [NII] λ 6583/H α (hereafter, the \'[NII] BPT diagram", top left), [OIII] λ 5007/H β vs. [SII] λλ 6717 , 6731/H α (hereafter, the \'[SII] BPT diagram", top right), and [OIII] λ 5007/H β vs. [OI] λ 6300/H α (hereafter, the \'[OI] BPT diagram", bottom). In each panel, the z ∼ 0 SDSS sample is plotted as a 2D greyscale histogram in the background, and the median sequence of nearby extragalactic H II regions from the CHAOS survey (e.g., Berg et al. 2015) is indicated as a dashed cyan curve. In the AURORA sample, we plot both detections and \'useful" limits (i.e., in which at least one line ratio has detections in both numerator and denominator).  \nThe AURORA sample recovers the well-known offset in the [NII] BPT diagram for z > 1 galaxies toward higher [NII] λ 6583/H α and [OIII] λ 5007/H β relative to the sequence of z ∼ 0 star-forming galaxies (Steidel et al. 2014; Shapley et al. 2015). Prior to JWST , the long-wavelength cut-off at ∼ 2 . 4 µ m for ground-based near-IR spectroscopy limited studies of this diagram to galaxies at z < 2 . 6. Early work  \nFigure 5. O 32 vs. R 23 diagram, corrected for dust. Symbols for AURORA galaxies are as in Figure 4. The z ∼ 0 SDSS histogram only includes galaxies classified as star-forming using the criteria of Kauffmann et al. (2003), given that AGNs do not neatly segregate in this diagram.  \n<!-- image -->  \nbased on composite JWST /NIRSpec spectra suggested a similar, constant offset at z > 2 (Sanders et al. 2023a), yet this result was based on a single stacked datapoint at each of z = 2 . 0 -2 . 7, z = 2 . 7 -4 . 0, z = 4 . 0 -5 . 0, and z = 5 . 0 -6 . 5. The AURORA sample includes 35 (17) star-forming galaxies at z = 1 . 4 -2 . 7 ( z = 2 . 7 -4 . 0) with individual [NII] BPT measurements (as well as a small additional sample at z > 4), enabling the first robust evolutionary measurement beyond z ∼ 2. In contrast to the initial stacked result from Sanders et al. (2023a), we find that the z = 2 . 7 -4 . 0 sample (green points) is even further offset on average (by ∼ 0 . 1 dex higher in typical [OIII] λ 5007/H β at fixed [NII] λ 6583/H α ) from the z ∼ 0 star-forming sequence than the z = 2 . 0 -2 . 7 sample (red points). As shown in Figure 3, the z = 2 . 0 -2 . 7 and z = 2 . 7 -4 . 0 samples significantly overlap in terms of SFR and M ∗ . Therefore, this relative offset in the [NII] BPT diagram suggests specific differences in the physical properties of H II regions as opposed to global sample differences, but larger sample numbers are still required to overcome potential sample variance.  \nThe sample of AURORA measurements at z > 4 (including 4 galaxies with all lines detected, and 7 with upper limits in [NII] λ 6583/H α ) is too small to draw firm conclusions regarding continued evolution in the [NII] BPT diagram. Early work with JWST /NIRSpec suggests that the typical [NII] λ 6583/H α line ratios for galaxies with stellar masses ranging from 10 8 -10 10 M ⊙ are quite modest at z > 4, ranging from 0.03 to 0.06 (Shapley et al. 2023b; Cameron et al.  \n2023). However, within AURORA, we observe a large dynamic range in [NII] λ 6583/H α , even at z > 6. We now require much larger samples in the [NII] BPT diagram at z > 4, with similar S/N to that obtained in AURORA.  \nNext we turn to the [SII] BPT diagram. As discussed in Sanders et al. (2017) and Shapley et al. (2019), z ∼ 0 SDSS galaxies and H II regions occupy different excitation sequences in this space, because of the added contribution of diffuse ionized gas to SDSS integrated spectra. This separation between the histogram of SDSS galaxies and the median sequence of CHAOS H II regions is clear in the top right panel of Figure 4. As justified in Section 2.9, z ∼ 2 galaxies should be compared with z ∼ 0 H II regions in the [SII] BPT diagram (Shapley et al. 2019). Based on their even higher Σ SFR, z > 2 galaxies should as well. Using the correct z ∼ 0 comparison dataset (i.e., the median H II region sequence), we find, as in Shapley et al. (2019), that z > 2 galaxies are significantly offset towards higher [OIII] λ 5007/H β and [SII] λλ 6717 , 6731/H α .  \nThe bottom panel of Figure 4 shows the [OI] BPT diagram. Detections of [OI] λ 6300 were only achieved recently for the first time for a small sample of star-forming galaxies at z > 1 (Clarke et al. 2023; Sanders et al. 2023b), based on ultra-deep ground-based spectroscopy with Keck/MOSFIRE. Stacked JWST/NIRSpec spectra have been used to detect [OI] λ 6300 at high redshift as well (Cameron et al. 2023; Strom et al. 2023). The AURORA sample contains individual [OI] λ 6300 detections of 46 star-forming galaxies, representing both an order-of-magnitude improvement in sample size and a significant (factor of several) increase in the typical [OI] λ 6300 detection S/N. Zhang et al. (2017) demonstrate that H II regions and DIG are characterized by distinct sequences of [OI] λ 6300/H α emission. In turn, we argue that high-redshift galaxies should be compared with z ∼ 0 H II regions in the [OI] BPT diagram, in analogy with our above approach to the [SII] BPT diagram. As in Clarke et al. (2023), the distribution of z > 1 . 4 AURORA galaxies is significantly offset towards higher in the [OI] λ 6300/H α relative to local H II regions.', '3.1.2. Ionization-Metallicity Diagram': "A common tool for investigating the combination of ionization parameter and metallicity in star-forming galaxies is the so-called O 32 vs. R 23 diagram (Figure 5). As described previously, because of the wavelength separation between [OII] λ 3727 and [OIII] λ 5007, we estimate dust-corrected line ratios assuming the Cardelli et al. (1989) dust curve. In the O 32 vs. R 23 diagram, the z ∼ 0 SDSS comparison sample only includes galaxies classified as star-forming according to the [OIII] λ 5007/H β vs. [NII] λ 6583/H α diagram. We note that z ∼ 0 galaxies and AGNs don't separate as cleanly in this space, with the AGNs identified in the z ∼ 0 [NII] BPT  \nFigure 6. [NeIII] λ 3869/[OII] λ 3727 vs. ([NeIII] λ 3869+[OII] λ 3727)/H δ , corrected for dust. Symbols for AURORA galaxies are as in Figure 4. The z ∼ 0 SDSS histogram only includes galaxies classified as star-forming using the criteria of Kauffmann et al. (2003), given that AGNs do not neatly segregate in this diagram. Also plotted are GN-z11 ( z = 10 . 6, cyan diamond) from Bunker et al. (2023), RXCJ2248-ID ( z = 6 . 1, purple pentagon) from Topping et al. (2024), and the stacked spectra at z = 5 . 5 -9 . 5 (orange squares) from Roberts-Borsani et al. (2024).  \n<!-- image -->  \ndiagram simply adding scatter roughly symmetrically around the star-forming sequence.  \nBecause useful O 32 vs. R 23 measurements can be obtained with coverage up to only ∼ 5000Å, we are able to plot the majority ( > 70%) of the z > 5 population in the AURORA sample here, in contrast to the BPT diagrams of the previous subsection. The AURORA sample shows a clear progression towards higher average O 32, in bins of increasing redshift. Also, there are are 20 galaxies with log( R 23) ≥ 1 (28% of AURORA galaxies with all the relevant lines detected), sharing almost no local counterparts (Sanders et al. 2016; Strom et al. 2017). We discuss these highR 23 galaxies further in Section 4.1.", '3.2. A Bluer Ionization-Metallicity Diagram for z > 10 Galaxies': 'JWST /NIRSpec dramatically extended the redshift range over which diagnostic rest-optical lines can be observed. However, at z > 10, where galaxies are now almost routinely being spectroscopically confirmed, even NIRSpec does not extend to strong, commonly-studied features such as [OIII] λ 5007, H α , and H β . Accordingly, in the context of understanding the physical properties of the ISM in the z = 10 . 6 galaxy, GN-z11, Bunker et al. (2023) introduce an emission-  \ne diagnostic diagram based on bluer rest-frame optical features. In the joint space of ionization and metallicity, the [NeIII] λ 3869/[OII] λ 3727 ratio serves the role of O 32, and, in place of R 23, there is ([NeIII] λ 3869+[OII] λ 3727)/H δ . These emission lines are all close in wavelength, and their ratios therefore relatively insensitive to dust attenuation. Bunker et al. place GN-z11 in the context of photoionization models and local galaxies to infer a low metallicity and high ionization parameter, and indirectly compare with a small sample of z > 5 galaxies with [NeIII]/[OII] and R 23 measurements (assuming a fixed ratio of H δ /H β ). GN-z11 appears to have similar line ratios to those of the stacked z ∼ 5 . 5 -9 . 5 NIRSpec prism (i.e., low-resolution) spectra of Roberts-Borsani et al. (2024). Thus far, however, there have been no direct comparisons between GN-z11 with individual galaxy spectra at high redshift and medium resolution.  \nIn Figure 6, we plot [NeIII] λ 3869/[OII] λ 3727 vs. ([NeIII] λ 3869+[OII] λ 3727)/H δ for the AURORA sample, along with the line ratios for GN-z11, the lowresolution stacks from Roberts-Borsani et al. (2024), and RXCJ2248-ID, a lensed z = 6 . 1 galaxy with an extremely high S/N NIRSpec R ∼ 1000 spectrum (Topping et al. 2024). Spanning from z = 1 . 4 to z = 7 . 0, at one extreme the AURORA star-forming galaxies overlap with the high [NeIII] λ 3869/[OII] λ 3727, low ([NeIII] λ 3869+[OII] λ 3727)/H δ tail of the SDSS z ∼ 0 distribution. At the other extreme, they extend towards GN-z11 and the z ∼ 5 . 5 -9 . 5 stacks. Within the AURORA sample, there exists a trend towards higher average [NeIII] λ 3869/[OII] λ 3727 and lower average([NeIII] λ 3869+[OII] λ 3727)/H δ at increasing redshift. Accordingly, GN-z11 and the z ∼ 5 . 5 -9 . 5 stacks are most consistent with the z > 5 portion of the AURORA sample. There are two exceptions at lower redshift: an auroralline target, GOODSN-30274 ( z = 1 . 800, red point), and a filler target GOODSN-917107 ( z = 4 . 773, blue point). These two galaxies have similar line ratios to those of GN-z11. Notably, the superior data quality of GOODSN-30274, including the detection of auroral [OIII] λ 4363 emission, may provide insights into the detailed nature of GN-z11. RXCJ2248ID has even more extreme [NeIII] λ 3869/[OII] λ 3727 and lower average([NeIII] λ 3869+[OII] λ 3727)/H δ than GN-z11, with [NeIII] λ 3869/[OII] λ 3727= 18 . 0 ± 5 . 7. Understanding the extreme phase of galaxy formation represented by RXCJ2248-ID constitutes an important goal for studies of the galaxy population during the reionization epoch.', '3.3. Longer-Wavelength Diagnostics: [SIII] and into the Near-IR': 'In addition to enabling studies of rest-optical emission lines for galaxies at z > 10, NIRSpec opens up a previously  \nunexplored window into the rest-frame near-IR spectroscopic properties of galaxies up to z ∼ 3.', '3.3.1. S 32 Ratio': 'The ratio of [SIII] λλ 9069 , 9532/[SII] λλ 6717 , 6731 ( S 32) provides a complementary probe of the ionization parameter and has been used to probe the ionization state of the ISM in z ∼ 0 galaxies (Kumari et al. 2021). As discussed in Sanders et al. (2020b), in concert with the [OIII] λ 5007/H β ratio, S 32 can distinguish among different scenarios for the evolution of physical conditions in high-redshift star-forming regions. Specifically, it is possible to differentiate between the effects of an evolving ionizing spectrum or ionization parameter at fixed nebular metallicity, both of which have been suggested in the literature as possible explanations for the observed emission-line ratios in distant galaxies. In groundbased near-IR spectroscopic observations (e.g., Sanders et al. 2020b), [SIII] λλ 9069 , 9532 was only detected in the K band for a handful of individual z ∼ 1 . 5 star-forming galaxies, and stacked composite spectra. In AURORA, we detect the S 32 ratio for 55 individual star-forming galaxies at z = 1 . 4 -4 . 5 (95% of the targets with coverage of the relevant lines). In Figure 7, we plot [OIII] λ 5007/H β vs. the dust-corrected S 32 ratio for AURORA galaxies, along with galaxy-integrated measurements from the MaNGA survey (Bundy et al. 2015) and the median sequence of H II regions from the CHAOS survey. As discussed previously, the CHAOS H II regions comprise the appropriate DIG-free comparison sample for AURORA, in order to infer the redshift evolution of the internal conditions in H II regions. The AURORA sample is unambiguously offset from the CHAOS sequence, towards lower S 32 at fixed [OIII] λ 5007/H β , within the range of overlapping [OIII] λ 5007/H β .', '3.3.2. Rest-frame Near-IR Line Ratios': 'Both Brinchmann (2023) and Calabrò et al. (2023) emphasized the potential of JWST /NIRSpec for detecting even longer rest-frame wavelength features with diagnostic power. Both works feature the [FeII] λ 1 . 257 µ m/Pa β ratio. Brinchmann (2023) pairs this ratio with He I λ 1 . 083 µ m/Pa γ and Calabrò et al. (2023) combines it with [SIII] λ 9532/Pa γ . Photoionization models with an input stellar ionizing spectrum suggest a typical maximum [FeII] λ 1 . 257 µ m/Pa β = 0 . 2 (Izotov & Thuan 2016). Higher [FeII] λ 1 . 257 µ m/Pa β occurs with photoionization by the harder ionizing spectrum of an AGN, or else by fast interstellar shocks (Calabrò et al. 2023; Brinchmann 2023). In the [SIII] λ 9532/Pa γ vs. [FeII] λ 1 . 257 µ m/Pa β diagram, Calabrò et al. (2023) determine the functional form for a curve discriminating between star-forming galaxies and AGNs (analogous to the curves in rest-optical BPT diagrams; e.g., Kauffmann et al. 2003).  \nIn Figure 8, we plot both He I λ 1 . 083 µ m/Pa γ (left) and [SIII] λ 9532/Pa γ (right) vs. [FeII] λ 1 . 257 µ m/Pa β for AU-  \nFigure 7. [OIII] λ 5007/H β vs. S 32 (dust-corrected). Symbols for AURORA galaxies are as in Figure 4. The grey z ∼ 0 2D histogram is based on data from the MaNGA survey, which has coverage of [SIII] λλ 9069 , 9532, and only includes galaxies classified as starforming using the criteria of Kauffmann et al. (2003). The cyan dashed curve represents the median sequence of H II regions from the CHAOS survey.  \n<!-- image -->  \nRORA galaxies at z ≤ 3. Local samples of such measurements are small (e.g., Izotov & Thuan 2016; Riffel et al. 2019), so we pair the AURORA measurements with Cloudy photoionization models (Ferland et al. 2017), generated assuming a BPASS binary stellar population with constant star formation over 10 8 yr for the ionizing spectrum (Eldridge et al. 2017). Model ionization parameter runs from log( U ) = -3 . 6 to -2 . 4; model gas-phase metallicity extends from Z neb = 0 . 05 -1 . 0 Z ⊙ ; and model stellar metallicity ranges from Z ∗ = 0 . 001 -0 . 014. Given that the Cloudy models were run assuming a solar abundance pattern, we scale the [FeII] λ 1 . 257 µ m flux by a factor Z ∗ / Z neb to reflect the well-established α -enhancement in the chemical abundance patterns of high-redshift star-forming galaxies. Furthermore a dust depletion factor of 20 was assumed for Fe (consistent with the model assumptions in Calabrò et al. 2023 and highredshift measurements in Pettini et al. 2002). The AURORA dataset of 23 star-forming galaxies with individual detections of all relevant lines represents an order of magnitude increase in the z > 1 sample size for this diagram. The 3 systems presented in Brinchmann (2023) have significantly larger error bars, and the sample presented in Calabrò et al. (2023) is dominated by individual non-detections. Furthermore, while the 3 z = 1 -3 systems from Brinchmann (2023) were presented as evidence for shock excitation, we find that the AU-  \nRORA sample can easily be explained by photoionization by massive stars. For example, the median [FeII] λ 1 . 257 µ m/Pa β is 0.2, which corresponds to [FeII] λ 1 . 257 µ m/H β =0 . 03, well within the range of stellar photoionization model predictions from Izotov & Thuan (2016) and Calabrò et al. (2023). We note that a significant systematic uncertainty in the comparison between measurements and models is associated with the assumed dust depletion factor for Fe. Independent constraints on the dust depletion of Fe and its dependence on metallicity are essential for robust model comparisons.', '4.1. The JWST/NIRSpec Case for α -Enhancement': 'Aphysical picture was beginning to emerge from a decade of ground-based near-IR spectroscopic studies, explaining the rest-optical nebular emission properties of star-forming galaxies observed during the epoch of peak star-formation in the universe. With JWST /NIRSpec, it is possible to elucidate this physical picture with significantly increased wavelength coverage (and in turn redshift range) and S/N (and in turn detections for representative individual, rather than stacked, measurements). As first proposed in Steidel et al. (2016), and expanded upon in, e.g., Strom et al. (2017), Shapley et al. (2019), Sanders et al. (2020b), Jeong et al. (2020), Runco et al. (2021), Topping et al. (2020), Clarke et al. (2023), and Cullen et al. (2021), star-forming galaxies at z ∼ 2 -3 have rest-optical emission line properties suggestive of photoionization by chemically-young (i.e., α -enhanced) massive stars. Accordingly, these stars have O/Fe ratios roughly 2-5 times higher than the ratio representative of the solar abundance pattern. This abundance pattern in massive stars corresponds to a harder ionizing spectrum (lower stellar iron abundance) at fixed nebular metallicity (gas-phase oxygen abundance).  \nThe AURORA data presented here confirm this physical picture in stunning detail. As described in Section 3.1.1 and Figure 4, the AURORA sample is not only offset with respect to the star-forming sequence in the [NII] BPT diagram, but also relative to the sequence of local H II regions in the [SII] and [OI] BPT diagrams. In Shapley et al. (2019) and Clarke et al. (2023), we investigated the [SII] and [OI] BPT diagrams, with significantly lower S/N, and, in the case of the [OI] diagram in particular, a significantly smaller sample of individual measurements. The depth of the AURORA spectra yields detections for the majority (indeed, almost all) of galaxies targeted, and therefore the samples of individual galaxies plotted are not only larger but also representative . These representative samples indicate shifts in the emissionline diagnostic diagrams that are consistent with photoionization by a harder ionizing spectrum at fixed nebular metallicity, relative to what is observed in low-redshift star-forming galaxies.  \n<!-- image -->  \nFigure 8. Rest-frame Near-IR Emission-line Ratios. Symbols for AURORA galaxies are as in Figure 4. Plotted in each panel are Cloudy (Ferland et al. 2017) model grids, assuming a BPASS (Eldridge et al. 2017) stellar population model with a constant star-formation history over 10 8 yr for the ionizing spectrum. A dust depletion factor of 20 is assumed for Fe (Calabrò et al. 2023; Pettini et al. 2002). Stellar metallicity (i.e., hardness of the ionizing spectrum) is indicated by the color of the model points, with Z ∗ = 0 . 001 , 0 . 004 , 0 . 008 , and 0 . 014 plotted, respectively, in dark blue, purple, turquoise, and orange. For connected model points of fixed ionization parameter ( U ) and stellar metallicity, increasing gasphase metallicity ( Z neb), ranging from 0.05 to 1.0 Z ⊙ , is indicated with increasing point size. Arrows are used to indicate the direction in which U and Z neb increase. The [FeII] λ 1 . 257 µ m line flux is also scaled by the factor Z ∗ / Z neb to reflect the α -enhanced chemical abundance patterns in distant star-forming galaxies. Left: He I λ 1 . 083 µ m/Pa γ vs. [FeII] λ 1 . 257 µ m/Pa β . In addition to the AURORA points and Cloudy model grids, we indicate three z ∼ 1 -3 galaxies from Brinchmann (2023) using black squares. Right: [SIII] λ 9532/Pa γ vs. [FeII] λ 1 . 257 µ m/Pa β . The same points are plotted as in the left panel, as well as the curve from Calabrò et al. (2023) intended to discriminate between star-forming galaxies and AGNs.  \n<!-- image -->  \nFurthermore, the high S/N of individual AURORA measurements enables us to discern differences between the average [NII] BPT properties of z ∼ 2 and z ∼ 3 galaxies for the first time. We find that our z ∼ 3 subsample is even more offset on average from the z ∼ 0 star-forming in the [NII] BPT diagram than our z ∼ 2 sequence is - a difference previously impossible to detect based on shallower, stacked measurements (Sanders et al. 2023a). This ongoing evolution before z ∼ 2 may provide important, complementary clues regarding the early star-formation histories of distant galaxies (Cullen et al. 2021; Kriek et al. 2019). Given the design of the AURORA program, its target redshift distribution is skewed towards z = 1 . 4 -4 (Figure 1). A clear extension of the current work is to assemble larger samples of JWST /NIRSpec spectra of z ∼ 4 -6 . 5 star-forming galaxies with the same wavelength coverage and at least comparable depth. Such data will enable a similarly in-depth investigation of the evolution in the [NII], [SII], and [OI] BPT diagrams back to cosmic dawn.  \nThe S 32 ratio is less familiar in the literature, based on the longer rest-frame wavelength of the [SIII] emission lines. Sanders et al. (2020b) demonstrated how the combination of [OIII] λ 5007/H β and S 32 line ratios could be used to  \nunambiguously determine the physical conditions in distant star-forming galaxies in relation to their local counterparts. Specifically, based on Cloudy photoionization modeling, Sanders et al. (2020b) showed that a shift in the [OIII] λ 5007/H β vs. S 32 perpendicular to the excitation sequence of z ∼ 0 H II regions (the appropriate comparison sample for z > 1 . 4 galaxies) could only be explained in terms of a harder ionizing spectrum at fixed nebular metallicity - not a larger ionization parameter (e.g., Hirschmann et al. 2017; Kashino et al. 2017). This earlier work was based on a sample of only 10 individual detections at z ∼ 1 . 5 with modest S/N, as well as stacked spectra containing ∼ 10 -30 z ∼ 1 . 5 galaxies. In contrast, the AURORA sample contains 55 individual S 32 detections at z = 1 . 4 -4 . 5 with significantly higher S/N - and robustly recovers the offset towards lower S 32 at fixed [OIII] λ 5007/H β indicative of a harder ionizing spectrum at fixed nebular metallicity.  \nAnother result from ground-based near-IR spectroscopic studies of z ∼ 2 -3 star-forming galaxies is the measurement of log( R 23) values in excess of unity (Sanders et al. 2016; Strom et al. 2017; Runco et al. 2021). There are virtually no corresponding measurements of local star-forming galaxies  \nat such high log( R 23) values (but see, e.g., Izotov et al. 2021). However, typical error bars on individual ground-based datapoints still make them consistent with the local star-forming sequence in the O 32 vs. R 23 diagram. With AURORA, the error bars on individual datapoints are small enough to demonstrate that 15 out of 20 of the galaxies in our sample with log( R 23) > 1 also have the lower 68th-percent confidence bound on log( R 23) greater than 1. Such high log( R 23) values with few if any local counterparts must be understood if R 23 is to be used as a robust metallicity indicator at z ≥ 4, as in early JWST/ NIRSpec work (Nakajima et al. 2023). Notably, hardening the ionizing spectrum at fixed nebular metallicity in photoionization models produces higher peak R 23 values at fixed O 32 in photoionization models (e.g. Steidel et al. 2016; Sanders et al. 2020a). Detailed photoionization modeling of the ionization parameter and ionizing spectrum of individual AURORA galaxies with direct metallicity estimates, as in Sanders et al. (2020a), will enable further examination of this question.', '4.2. Dust Depletion of Refractory Elements': 'As presented in Section 3.3.2, JWST /NIRSpec enables measurements of rest-frame near-IR diagnostic line ratios for AURORA galaxies at z ≤ 3. The range of [FeII] λ 1 . 257 µ m/Pa β line ratios in the AURORA sample appears can be explained in terms of photoionization by massive stars - not shocks or an AGN (Brinchmann 2023). However, the main source of uncertainty in comparing these new rest-frame near-IR line ratios with photoionization models is the assumption of the degree of dust depletion for the refractory element, iron. We assumed a factor of 20 depletion for iron, roughly consistent with previous work (Calabrò et al. 2023; Pettini et al. 2002). However, better observational constraints on the level of dust depletion for iron in distant starforming galaxies is an essential component of interpreting these new line ratios, building on the work of Pettini et al. (2002), Jones et al. (2018), and De Cia et al. (2016). In particular, a description of dust depletion incorporating its metallicity dependence (Izotov et al. 2006; Rodríguez & Rubin 2005) would represent an advance over the simple description presented here and in other recent work (e.g., Calabrò et al. 2023).', '4.3. Why Emission-line Diagrams?': "NIRSpec provides a significantly improved data quality for individual galaxy spectra compared to previous instrumentation. NIRSpec spectra from the AURORA program include detections of faint auroral lines for roughly half the sample ( ∼ 50 galaxies) enabling direct metallicities (Sanders et al. 2024, in prep.). In some cases, these spectra feature a large suite of up to tens of nebular emission lines spanning from the rest-frame blue optical into the rest-frame near-IR  \n(e.g., Figure 2), which can be simultaneously modeled to infer the physical properties of the ionized ISM. This transformational advance in S/N raises a question of the very relevance of the emission-line diagnostic diagrams forming the focus of this paper. However, while the AURORA data demonstrate what is possible to achieve for individual galaxy spectra with NIRSpec with 24-hour exposure times, the majority of NIRSpec spectra collected to date do not have the same depth. These shallower NIRSpec spectra (Bunker et al. 2023; D'Eugenio et al. 2024; Shapley et al. 2023a) do contain detections of subsets of the strongest rest-frame optical emission lines, and therefore can be situated in emission-line diagnostic diagrams. Accordingly, the AURORA dataset for which both detailed modeling of individual spectra and plots of emission-line diagnostic diagrams are possible - provides a physical basis for interpreting the emission-line diagnostic diagrams of the much larger samples of lower-S/N JWST /NIRSpec spectra that have been and will be collected in the years to come. In the current work, we present the AURORA emission-line diagrams. In future work, we will present detailed modeling of the ISM physical properties implied by individual AURORA spectra. Combining these complementary approaches is essential for charting the evolution of galaxy chemical enrichment and ISM physical conditions over cosmic time.", 'ACKNOWLEDGEMENTS': 'This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-03127 for JWST. The specific observations analyzed can be accessed via DOI: 10.17909/hvne-7139. We also acknowledge support from NASA grant JWST-GO-01914. FC acknowledges support from a UKRI Frontier Research Guarantee Grant (PI Cullen; grant reference: EP/X021025/1). ACC thanks the Leverhulme Trust for their support via a Leverhulme Early Career Fellowship. CTD, DJM, RJM, and JSD acknowledge the support of the Science and Technology Facilities Council. JSD also acknowledges the support of the Royal Society through a Royal Society Research Professorship. RD acknowledges support from the Wolfson Research Merit Award program of the U.K. Royal Society. KG acknowledges support from the Australian Research Council Laureate Fellowship FL180100060. MK acknowledges funding from the Dutch Research Council (NWO) through the award of the Vici grant VI.C.222.047 (project 2010007169). PO acknowledges the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number MB22.00072, as well as from the Swiss National Science Foundation (SNSF)  \nthrough project grant 200020\\_207349. AJP was generously supported by a Carnegie Fellowship through the Carnegie Observatories. Finally, we thank members of the JADES  \nteam for assistance with target selection in the GOODS-N field.', 'REFERENCES': 'Zhang, K., Yan, R., Bundy, K., et al. 2017, MNRAS, 466, 3217'}