The exceptional spectra of the most luminous z gt 10 sources observed so far have challenged our understanding of early galaxy evolution requiring a new observational benchmark for meaningful interpretation. As such we construct spectroscopic templates representative of highredshift starforming populations using 482 confirmed sources at z  5.012.9 with JWSTNIRSpec prism observations and report on their average properties. We find z  511 galaxies are dominated by blue UV continuum slopes   2.3 to 2.7 and reduced Balmer indices characteristic of dustpoor and young systems with a shift towards bluer slopes and younger ages with redshift. The evolution is mirrored by ubiquitous C III detections across all redshifts restframe equivalent widths of 514  which increase in strength towards early times. Restframe optical lines reveal elevated ratios O32  731 R23  58 and Ne3O2  12 and subsolar metallicities logOH  7.37.9 typical of ionization conditions and metallicities rarely observed in z  0 populations. Within our sample we identify 57 Ly emitters which we stack and compare to a matched sample of nonemitters. The former are characterized by more extreme ionizing conditions with enhanced C III C IV and He II  O III line emission younger stellar populations from Balmer jumps and a more pristine interstellar medium seen through bluer UV slopes and elevated restframe optical line ratios. The novel comparison illustrates important intrinsic differences between the two populations with implications for Ly visibility. The spectral templates derived here represent a new observational benchmark with which to interpret highredshift sources lifting our constraints on their global properties to unprecedented heights and extending out to the earliest of cosmic times.
['Galaxy evolution', 'High-redshift galaxies', 'Reionization', 'Metallicity', 'Stellar ages', 'Early universe', '594', '734', '1383', '1031', '1581', '435', 'Astrophysics - Astrophysics of Galaxies']
Between the Extremes A JWST Spectroscopic Benchmark for Highredshift Galaxies Using 500 Confirmed Sources at z  5
{'Between the Extremes: A JWST Spectroscopic Benchmark for High-redshift Galaxies Using ∼ 500 Confirmed Sources at z ⩾ 5': "Guido Roberts-Borsani , 1 Tommaso Treu , 2 Alice Shapley , 2 Adriano Fontana , 3 Laura Pentericci , 3 Marco Castellano , 3 Takahiro Morishita , 4 Pietro Bergamini , 5, 6 and Piero Rosati 6, 7  \n1 Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland 2 Department of Physics and Astronomy, University of California, Los Angeles, 430 Portola Plaza, Los Angeles, CA 90095, USA 3 INAF - Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone, Rome, Italy 4 IPAC, California Institute of Technology, MC 314-6, 1200 E. California Boulevard, Pasadena, CA 91125, USA 5 Dipartimento di Fisica, Universit'a degli Studi di Milano, Via Celoria 16, I-20133 Milano, Italy 6 INAF - Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via Gobetti 93/3, I-40129 Bologna, Italy 7 Dipartimento di Fisica e Scienze della Terra, Universit'a degli Studi di Ferrara, Via Saragat 1, I-44122 Ferrara, Italy", 'ABSTRACT': 'The exceptional spectra of the most luminous z > 10 sources observed so far have challenged our understanding of early galaxy evolution, requiring a new observational benchmark for meaningful interpretation. As such, we construct spectroscopic templates representative of high-redshift, star-forming populations, using 482 confirmed sources at z = 5 . 0 -12 . 9 with JWST/NIRSpec prism observations, and report on their average properties. We find z = 5 -11 galaxies are dominated by blue UV continuum slopes ( β = -2 . 3 to -2 . 7) and reduced Balmer indices, characteristic of dust-poor and young systems, with a shift towards bluer slopes and younger ages with redshift. The evolution is mirrored by ubiquitous C III ] detections across all redshifts (rest-frame equivalent widths of 5 -14 ˚ A), which increase in strength towards early times. Rest-frame optical lines reveal elevated ratios ( O 32 = 7 -31, R 23 = 5 -8, and Ne 3 O 2 = 1 -2) and subsolar metallicities (log (O/H)= 7 . 3 -7 . 9), typical of ionization conditions and metallicities rarely observed in z ∼ 0 populations. Within our sample, we identify 57 Ly α emitters, which we stack and compare to a matched sample of nonemitters. The former are characterized by more extreme ionizing conditions with enhanced C III ], C IV , and He II +[O III ] line emission, younger stellar populations from Balmer jumps, and a more pristine interstellar medium seen through bluer UV slopes and elevated rest-frame optical line ratios. The novel comparison illustrates important intrinsic differences between the two populations, with implications for Ly α visibility. The spectral templates derived here represent a new observational benchmark with which to interpret high-redshift sources, lifting our constraints on their global properties to unprecedented heights and extending out to the earliest of cosmic times.  \nKeywords: galaxies: high-redshift, galaxies: ISM, galaxies: star formation, cosmology: dark ages, reionization, first stars', '1. INTRODUCTION': "Determining the properties of sources that emerged in the early Universe and transformed it into a fully ionized one represents one of the major milestones of modern extragalactic astronomy, and a prerequisite toward charting the evolution of galaxies from Cosmic Dawn to the present day. The most fundamental of these are the  \nages, masses, and metallicities of their stellar populations, the ionizing conditions and metal and dust contents of their interstellar medium (ISM), and the presence or otherwise of supermassive black holes. Stellar and emission-line features at rest-frame optical wavelengths ( λ ≃ 3800-6700 ˚ A) underpin most of these properties and are critical in determining if and how they evolve between sources at the end of the Reionization Era ( z ≲ 5 . 3 -5 . 8; Eilers et al. 2018; Bosman et al. 2022; Zhu et al. 2023) and those firmly within it ( z ≳ 7 . 5; Treu et al. 2013; Mason et al. 2019b). Direct probes of  \nthe rest-frame optical spectrum in representative z ⩾ 5 sources, therefore, represents a necessary step forward in deducing whether Reionization Era sources reflect distinct populations compared to their lower-redshift counterparts.  \nUntil recently, inferences of such properties in z ∼ 5 -13 sources - probing the aftermath, the conclusion, and the unfolding of the reionization process - have relied primarily on photometric studies of the rest-frame ultraviolet (UV) with imaging from the Hubble Space Telescope (HST) and (when available) the rest-frame optical from the Spitzer Space Telescope (e.g., Labb'e et al. 2013; Smit et al. 2014; Bouwens et al. 2015; Strait et al. 2020; Roberts-Borsani et al. 2022; Yang et al. 2022; Whitler et al. 2023b), revealing properties rarely seen in the local Universe. Analyses of individual sources at z ∼ 6 -9 have painted a picture of generally compact (a few kiloparsecs or even subkiloparsec), young ( ∼ 1 -50 Myr) and star-forming ( ∼ 1-100 M ⊙ /yr) sources with subsolar ( ∼ 0 . 2 -0 . 8 Z ⊙ ) metallicities, low dust contents ( A v < 1 mag), and extreme rest-frame optical line emission ([O III ] λλ 4960,5008 ˚ A + H β equivalent widths (EWs) in excess of > 1000 ˚ A). Stacks of Spitzer photometry (e.g., Labb'e et al. 2013; Stefanon et al. 2023), yielding higher signal-to-noise ratio (SNR) measurements out to z ∼ 10, reaffirmed such a picture and provided an early benchmark for the interpretation of high-redshift sources.  \nThe advent of > 2 µ m filters from the James Webb Space Telescope's (JWST) Near Infrared Camera (NIRCam) has extended the redshift frontier for galaxy formation studies, with the enhanced sensitivity and narrower bandpasses revealing bluer, compact, and starforming sources with intense line emission out to z ∼ 11 -13 (e.g., Topping et al. 2022, 2023; Cullen et al. 2023; Endsley et al. 2023; Laporte et al. 2023; Leethochawalit et al. 2023; Ono et al. 2023; Santini et al. 2023; Withers et al. 2023). However, such analyses have generally relied on photometric redshifts, spectral energy distribution (SED) modeling using synthetic templates, and scaling relations calibrated to low-redshift samples.  \nSpectroscopy remains the ultimate tool for confirmation of galaxy redshifts and the interpretation of their global properties. While spectroscopic studies over statistical and magnitude-limited z ⩾ 5 Lyman break galaxy (LBG) samples have already been attempted with ground-based facilities (e.g., Treu et al. 2013; Jung et al. 2018; Pentericci et al. 2018; Mason et al. 2019b; Endsley et al. 2021a), these have typically targeted the rest-frame UV Ly α emission line, which is resonantly scattered from the line of sight in a neutral medium (with exceptions over smaller samples; e.g., Stark et al.  \n2017; Laporte et al. 2017a; Mainali et al. 2018, or gas tracers in the far-infrared; e.g., Bouwens et al. 2021). For the few luminous sources where spectroscopic redshifts could be determined, some characterization of their rest-frame optical properties (e.g., metallicities, specific star formation rates (SFRs), [O III ]+H β EWdistributions) has been possible with Spitzer/IRAC photometry (Hashimoto et al. 2018; Laporte et al. 2017b; Endsley et al. 2021a,b; Roberts-Borsani et al. 2023a). Direct spectroscopic characterization of the rest-frame optical has however, until recently, remained beyond the reach of z ⩾ 5 sources and achievable only in groundbased surveys of lower-redshift ( z ∼ 1 -4) populations (e.g., the MOSDEF and VANDELS surveys; Kriek et al. 2015; McLure et al. 2018). As such, constraints over the full rest-frame UV to optical of high redshifts sources remain sparse, and their average spectra and properties unknown.  \nJWST's Near Infrared Spectrograph (NIRSpec) has proven remarkably efficient at confirming galaxy candidates out to z ∼ 14 and affirm their underlying properties across cosmic time (e.g., Arrabal Haro et al. 2023a; Curtis-Lake et al. 2023; Harikane et al. 2023b; Mascia et al. 2023a; Morishita et al. 2023a; Nakajima et al. 2023; Roberts-Borsani et al. 2023b; Sanders et al. 2023a; Tang et al. 2023; Williams et al. 2022; Curti et al. 2023; Carniani et al. 2024; Roberts-Borsani et al. 2024). The revelation of spectacular and perplexing spectra of GN-z11 (Bunker et al. 2023a) and GHZ2 (Castellano et al. 2024), however, has placed our already limited understanding of galaxy evolution at the earliest of times in doubt. Whether these sources comprise atypical or somewhat representative galaxies for their luminosities and epochs remains an open question, one where comparisons to photoionization models or lowredshift analogs have proved challenging. Indeed, comparisons to observations of high-redshift sources offer a more direct framework with which to interpret such observations, however the generally low SNR of individual spectra, relatively small samples sizes, and scatter in their inferred properties has made such comparisons challenging. Determining the 'typical' spectrum of a high-redshift source (e.g., Shapley et al. 2003) with which to establish a representative benchmark and template to place individual sources into context, therefore, represents an ambitious challenge and one which would significantly enhance our capabilities in the interpretation of highz observations.  \nThe primary objectives of this paper are to (i) study the average spectra of a large number of high-redshift galaxies in order to understand their properties at a population level, while (ii) simultaneously providing high  \nSNR composite templates of observed galaxy spectra to the community. The paper is organized as follows. In Section 2 we present the spectroscopic and photometric data sets used in this study, and present our sample selection and stacking methods in Section 3. In Section 4 we present stacks of galaxy spectra as a function of global properties, and report on their main spectral features and line ratios. Section 5 presents the derivation of global galaxy properties from those spectra and places them into context with comparisons to known objects and literature studies. Finally, we summarize our findings and present our conclusions in Section 6. All stacked spectra presented in this study are available on the GLASS-JWST website 1 , while redshifts are listed in Tables 3-4 and line fluxes in Table 5 in the Appendix.", '2.1. NIRSpec Prism Spectroscopy': "Given its large wavelength coverage, sensitivity to both line and continuum emission, and the number of observed sources, this analysis focuses on observations and objects where NIRSpec prism spectra are available. Our data sets derive primarily from publicly available programs or data releases (at the time of writing) using the multishutter array (MSA), and include observations over the Abell 2744, RXJ2129, and MACS 0647 lensing clusters through the GLASS-JWST Director's Discretionary Time (DDT) proposals (PID 2756; see Morishita et al. 2023a and Roberts-Borsani et al. 2023a), GO 3073 (Castellano et al. 2024), UNCOVER (PID 2561; see Bezanson et al. 2022), DDT 2767 (see Williams et al. 2022), and GO 1433 (see Hsiao et al. 2023), respectively, as well as the CANDELS-EGS, CANDELS-UDS, GOODS-North, and GOODS-South blank fields through CEERS (PID 1345; see Bagley et al. 2023; Finkelstein et al. 2023) and its associated DDT (PID 2750; see Arrabal Haro et al. 2023b), GO 2565, and the JADES survey (GTO 1181 and GTO 1211 for GOODS-North, GTO 1180, GTO 1210, and GO 3215 for GOODS-South; Bunker et al. 2023b), respectively. In total, the data sets amount to 13 programs over seven unique fields, each affording low resolution ( R ≃ 30 -300) spectroscopy with a continuous wavelength coverage from ∼ 0 . 6 -5 . 3 µ m. A list of the data sets along with their targeted fields is presented in Table 1.  \nWith the exception of the JADES GTO 1210 program, we download, reduce, and extract all public NIRSpec/MSA spectra following the procedure described  \nin Morishita et al. (2023a). In short, the countrate ( rate.fits ) exposures are downloaded from the MAST archive and reduced with the msaexp 2 code, which adopts the official STScI pipeline modules in addition to custom routines (including flat-fielding, snowball masking, 1/ f removal, and others). For a given object, the code then extracts the 2D spectrum from each exposure, rectifies and drizzles it to a common pixel grid, and background subtracts it using its associated dithers before combining these via an inverse-varianceweighted mean with outlier rejection. Individual exposures affected by catastrophic events (e.g., transient optical shorts or significant contamination) are discarded prior to combining. The final, rectified 2D spectrum is then used together with an optimal extraction procedure (Horne 1986) to extract a 1D spectrum with associated uncertainties. In the case of the JADES GTO 1210 spectra, we use the data products provided by the collaboration. All spectra are visually inspected for quality control and the 1D extraction repeated or optimized if required.", '2.2. Uniformly Processed HST and JWST Photometry': "We also leverage deep JWST/NIRCam and/or HST/ACS+WFC3 imaging for spectrophotometric scaling and verification. The images are downloaded from the Grizli v6 Image Release 3 , a set of publicly available NIRCam and HST images deriving from public programs that were uniformly reduced over a consistent astrometric reference frame with the STScI pipelines and the Grizli code (Brammer 2023). The subset of NIRCam images used here derive from the programs listed in Table 1 and the ERS 1324 program (GLASS-JWST; Treu et al. 2022), while the HST data derive from legacy surveys such as the Frontier Fields (Lotz et al. 2017), RELICS (Coe et al. 2019), and CANDELS (Grogin et al. 2011) programs.  \nTo derive photometric catalogs we use Grizli 's implementation of the Python wrapper for SExtractor , SEP (Bertin & Arnouts 1996; Barbary 2016), in conjunction with an infrared detection image derived from an inverse-variance-weighted mean of the NIRCam F115W, F277W and F444W images, where available. In some cases, the reduced NIRSpec spectra lie outside of an NIRCam footprint, and in such cases we adopt the photometry from the HST images (using a stacked WFC3/F125W+F160W infrared image as reference). For the Extended Groth Strip (EGS) field, we ex-  \nTable 1. The JWST/NIRSpec Prism Data Sets Used in This Study and the Number of Confirmed z ⩾ 5 Sources Identified Here at the Time of Writing.  \nNote -The integration times refer to the median (left) and maximum (right) integration times (in hours) for the parent sample of z ⩾ 5 sources used in this study.  \ntend the covered area with public CANDELS images 4 (ACS/F606W+F814W and WFC3/F125W+F160W), since these extend slightly beyond the Grizli images. Although Spitzer/IRAC data also exist for a large fraction of the total area considered, we opt not to use them given the comparatively shallower depth and poorer spatial resolution, which could lead to significant confusion in crowded fields. For galaxies residing in the JADESDeep survey, we utilize the publicly released catalogs from Rieke & the JADES Collaboration (2023) and Eisenstein et al. (2023). In all cases we derive and/or utilize aperture photometry measured in fixed, 0 '' . 5 diameter circular apertures. The choice of aperture size is motivated to encapsulate a maximum of target flux while minimizing potential contamination - we find the vast majority of our galaxy sample (see Section 3) are adequately encompassed without leaving room for contamination from nearby compact sources or light from extended objects in e.g., lensing clusters.", '3.1. Redshift Measurements of z ⩾ 5 Sources': 'We verify the redshifts of all NIRSpec/prism galaxies identified and extracted in Section 2.1 using the 1D spectrum and the redshift-fitting module of msaexp , which evaluates a χ 2 goodness of fit as a function of redshift using a series of EAzY (Brammer 2021) galaxy and line templates at the nominal NIRSpec/prism resolution. We fit each galaxy using a redshift range between z = 0 -15 and visually inspect the resulting solu-  \nions. In total, our procedure yields 439 spectra at z ⩾ 5 with secure redshifts, many of which have already been reported in the literature (e.g., Roberts-Borsani et al. 2023a; Tang et al. 2023; Fujimoto et al. 2023; Williams et al. 2022; Arrabal Haro et al. 2023a; Wang et al. 2022; Morishita et al. 2023a; Mascia et al. 2023a; Nakajima et al. 2023; Harikane et al. 2023a; Napolitano et al. 2024) while a large number (approximately 69%) represent novel confirmations. We further supplement our sample of galaxies with a further 48 sources confirmed at z ⩾ 5 in the JADES-Deep survey and listed in Bunker et al. (2023b), increasing the total to 487 sources. For those sources where multiple observations or multiple lensed images exist (e.g., Roberts-Borsani et al. 2023b; Hsiao et al. 2023), we combine the (magnification-corrected, where applicable) observations via an inverse-varianceweighted mean, thus counting them as a single observations or source. In total, this results in a final sample of 482 confirmed sources.  \nA requirement of our selection is to be able to scale sources for an overall slit loss in the rest-frame UV (see Section 3.2 below). A number of objects in our sample suffer from a significant gap in the relevant portion of their NIRSpec spectrum, a lack of public photometry, or insufficient continuum flux to be able to do so. Therefore, those sources are discarded from our primary analysis and our fiducial sample is reduced to 454 unique sources. The full list of spectroscopic redshifts for our fiducial sample is tabulated in Table 3, along with several key spectrophotometric properties. Finally, in Table 4 we also report the redshifts for those sources where a photometric scaling could not be performed, as well as for an additional 34 objects for which the redshift is not  \nsecurely constrained - while the msaexp code favors a highz solution, these sources either lack strong spectral features or SNR to unambiguously confirm their redshifts and are thus not included in any subsequent analysis.', '3.2. Bad Pixel Masking and Quality Control': 'We consider a number of additional checks for quality control. Firstly, all (1D and 2D) spectra are visually inspected by two authors (G.R.B. and T.T.) and pixels corresponding to clear contamination or artifacts are masked - examples of these include detector gaps, bad background subtraction, and/or contamination from nearby objects, among others. Furthermore, we take steps to ensure the spectra are not affected by any lingering instrumental artifacts not completely accounted for by the pipeline. The main ones considered here (where relevant) are (i) the scaling of the spectrum to a rest-frame UV flux in order to account for the overall slit loss, and (ii) determining whether wavelengthdependent slit loss effects are present or significant in the prism observations. One limitation we do not account for is any potential difference between the photometry and spectroscopy attributable to spatially varying galaxy properties. For (i), we utilize the extracted HST or NIRCam photometry described in Section 2.2 to account for the overall slit loss by scaling the spectrum to the closest filter (preferentially NIRCam but also HST if the former is lacking) corresponding to a rest-frame wavelength of 1750 ˚ A, where the photometry probes the rest-frame UV but is unaffected by Lymanα emission or the Lyman break. For (ii), we use the available NIRCam photometry longward of the Lyman break to evaluate whether wavelength-dependent slit loss is a prevalent effect in prism observations. A detailed discussion of this residual slit loss and the derivation of an average correction function is presented in the Appendix, from which we find minor losses of 1-2% in < 2 µ m NIRCam bands and up to 10-20% losses (on average) in redder bands. We derive a correction function from these average slit losses and correct all our composite spectra with this function prior to any analysis. As with our composite spectra, we render the correction function publicly available for community use.', '3.3. Magnification Factors of Lensed Galaxies': 'Asignificant fraction ( ≃ 37%) of our galaxy sample derives from lensing cluster fields (see Table 1) and as such require correction for lensing magnification. For objects behind the A2744 cluster, we use the magnification maps from Bergamini et al. (2023), which make use of the recently confirmed z spec . = 9 . 75 triply imaged galaxy by  \nRoberts-Borsani et al. (2023b) and 149 multiple images for which a large number include spectroscopic redshifts. For objects behind the RXJ2129 cluster we adopt the magnification maps of Caminha et al. (2019), while for objects behind the MACS 0647 cluster we adopt the magnification maps of Meena et al. (2023, private communication). The spectra and photometry of all lensed sources in our sample of confirmed galaxies are corrected for their adopted magnification factors prior to any analysis, and the M UV values are derived using the measured flux and standard deviation of the resulting spectra in a wavelength window λ rest = 1500 -2000 ˚ A. The distribution of our fiducial sample in the z -M UV plane is shown in Figure 1.', '3.4. Identification of Active Galactic Nuclei': 'One of the most surprising findings from early JWST results has been the revelation of significant numbers of AGN in the early Universe (e.g., Harikane et al. 2023b; Larson et al. 2023; Greene et al. 2023; Goulding et al. 2023; Maiolino et al. 2023; Kocevski et al. 2023). While optimal identification methods and the nature of some objects has yet to reach consensus, we err on the side of caution and flag all AGN objects identified in the literature (i.e., objects identified as AGN in Harikane et al. 2023b, Larson et al. 2023, Greene et al. 2023, Goulding et al. 2023, Maiolino et al. 2023, Kocevski et al. 2023). Furthermore, for all objects in our sample not included in those analyses, we evaluate via visual inspection whether a two-component (i.e., broad plus narrow) model is clearly required to reproduce the H α profile when the [O III ] λλ 4960,5008 ˚ A profile is well fit by a single, narrow component (characteristic of Type 1 AGN). We opt not to extend such an analysis to measurements of other Balmer lines (e.g., H β , H γ , and H δ ) since they often lack sufficient SNR in individual spectra for robust fits. This amounts to 27 AGN candidates (identified here and in the literature), which are flagged (also listed in the bottom half of Table 3) and not utilized in the subsequent analyses. Of course, such a procedure does not preclude some degree of AGN (broad- or narrowline) contamination remaining in our sample, particularly at redshifts z > 7 where H α moves beyond the NIRSpec coverage, however these are likely to be weak AGN whose identifications are beyond the typical (per pixel) SNR ratios of the individual spectra used here.', '3.5. A Representative Sample of z ⩾ 5 Galaxies': "With our fiducial sample of 427 confirmed starforming galaxies in hand, we assess the degree of our targets as a representative sample of high-redshift sources by evaluating whether they adequately populate wellknown planes of global galaxy properties. One of the  \nFigure 1. The sample of spectroscopically confirmed sources at z ⩾ 5 compiled in this work (listed in Tables 3-4 in the Appendix), with continuum measurements available for spectrophotometric scaling in the rest-frame UV. Sources identified in blank fields are shown in blue, while lensed sources (approximately 37% of the parent sample) are plotted in purple. Sources considered to be AGN in this study are indicated with black squares.  \n<!-- image -->  \nmost ubiquitous correlations in galaxy evolution is the so-called 'main sequence' (Noeske et al. 2007; Speagle et al. 2014), i.e., that between galaxy SFRs and stellar masses ( M ∗ ) which is expected to reflect the statistical accumulation of stellar mass from smooth gas accretion and minor mergers. While individual sources may temporarily lie above or below the main sequence through stochastic levels of star formation (e.g., from mergers), generally speaking a source's position on the SFRM ∗ diagram dictates whether it represents a 'typical' source at a given epoch.  \nTo illustrate the regions of parameter space populated by our sample, we derive global galaxy properties following the procedure outlined by Morishita et al. (2023b). In short, we utilize a best-fit continuum model (see Section 3.7 below) at the redshifts of the photometrically scaled and magnification-corrected NIRSpec spectra to obtain stellar masses and SFRs, the latter derived using the rest-frame UV luminosity ( L UV ) measured at ∼ 1600 ˚ A from the best-fit SED. To account for dust attenuation, each L UV is corrected using the measured UV continuum slope and the A 1600 -β relation from Meurer et al. (1999) - in principle, such a relation is valid only down to slopes of approximately -2 . 2, whereas a large number of our sample display even bluer slopes. In such  \ncases, we assume negligible dust attenuation and do not correct the UV luminosity. The attenuation-corrected luminosity is then converted to an SFR using the Kennicutt (1998) relation, scaled to a Chabrier (2003) initial mass function (IMF):  \nSFR[M ⊙ yr -1 ] = 8 . 8 × 10 -29 L UV [erg s -1 Hz -1 ] (1)  \nWe plot the resulting values in four redshift intervals (5 < z < 7, 7 < z < 9, 9 < z < 11, and z > 11) in Figure 2, along with a best-fit straight line to the data in the lowest bin where the numbers dominate. We find the sources in each of the four redshift intervals populate expected regions of the parameter space, following a clear trend of increasing SFR with stellar mass. While the number of available data points declines rapidly with redshift, the consistency of the data points in each bin relative to the measured z ≃ 5 -7 main-sequence line (log SFR [ M ⊙ / yr]=0.66 log M ∗ [ M ⊙ ] -5 . 50) suggests little evolution in the main sequence between the redshifts probed here. We also compare our data points to the Speagle et al. (2014) z ∼ 6 main-sequence relation, finding good consistency between our points and the slope of their relation, with minor differences in the normalization (the Speagle et al. 2014 relation resides  \nFigure 2. The star-forming main sequence as defined by our sample of confirmed sources, in four redshift bins. Stellar masses are derived from best-fit gsf SED models of the stellar continuum, while SFRs are derived from attenuationcorrected UV luminosities at ∼ 1600 ˚ A. A best fit straight line to the 5 < z < 7 data is shown in each panel, along with the Speagle et al. (2014) relation at z ∼ 6 in the top left panel.  \n<!-- image -->  \napproximately +0.35 dex above ours, possibly due to potential differences in SFR calibrations, magnification corrections, etc.) which we note are both small and effectively constant with stellar mass, indicating that our points scatter symmetrically (across SFR) around the shifted version of the relation at fixed stellar mass, as expected for a representative sample of star-forming galaxies. The consistency of our inferred galaxy properties with the highz main sequence, as well as the large range of absolute magnitudes probed in Figure 1, supports the labeling of this compilation as a representative galaxy sample.", '3.6. A Note Regarding Selection Effects': 'The goal of the present paper is to construct representative composite spectra of high-redshift galaxies. Whether a composite spectrum is considered representative of the general galaxy population depends heavily on the variety of objects used in a stack. It is thus important to note that at the time of writing the selection functions for the NIRSpec programs used are extremely heterogenous, often unpublished, and therefore virtually impossible to characterize. This caveat means that one cannot define precisely in what sense our stacks (e.g., especially those at z > 10) are representative of the general population. They are certainly not flux-limited  \nsamples, nor simple color-selected samples. However, a posteriori characterization (Figures 1 and 2) shows that our individual sources lie on the main sequence, and span more or less uniformly a range in M UV between -21 mag and -17 mag. Furthermore, the properties of individual sources are in excellent agreement with independently derived properties from the literature. Therefore, the stacks presented here should be considered broadly representative of typical star-forming galaxies in the appropriate M UV and redshift ranges, although this caveat should be kept in mind when using them.  \nFurthermore, we also acknowledge that not all prism spectra have sufficient SNR and spectral resolution to separate, e.g., AGN and Ly α emitters (LAEs; see Sections 3.4 and 5.3) from their star formation-dominated and/or Ly α -attenuated counterparts. Such separations are limited by the data and some contamination may still be present in our stacks (although these are unlikely to be strong AGN or LAEs). Our efforts represent a first step with JWST/NIRSpec, and our stacks will continue to improve in SNR, representation, and diversity of galaxy spectral type as additional data sets become public and selection methods are better understood.', '3.7. Stacking Method and Spectral Modeling': 'We create composite galaxy spectra at high redshift through a stacking of prism data based on a number of galaxy properties. To construct a stack, each galaxy spectrum (already photometrically scaled and magnification corrected) is first de-redshifted then interpolated over a common wavelength grid (with 10 ˚ A intervals) and added to the stack after normalization to unity at its rest-frame UV flux (between λ rest = 1500 -2000 ˚ A). Composite fluxes and associated uncertainties are derived from the nonzero median and the semidifference of the 16th and 84th percentiles of the distribution in each wavelength bin (the latter divided by a square root of N i term, where N i is the number of nonzero fluxes in the stack at a given wavelength element, i ), respectively, after which the resulting spectrum and uncertainties are shifted and scaled to the median redshift and rest-frame UV magnitude of the stack. We do not apply any weighting to the fluxes in a given stack, in order not to bias toward apparently bright objects or those with longer integration times. An example composite spectrum (1D and 2D) derived using the full sample of z ⩾ 5 star-forming galaxies (and left unscaled) is shown in Figure 3, for illustration.  \nFigure 3. The composite (unscaled) NIRSpec prism spectrum (2D above, 1D below) of all z ⩾ 5 star-forming galaxies from this study, in the rest frame. Several notable features are seen, from a clear Lymanα break at λ rest ≃ 1216 ˚ A highlighting the neutral intergalactic medium (IGM) within which most of the objects reside (with some excess UV flux emanating from z < 6 sources where the IGM becomes ionized), to weak UV lines (the blended He II λ 1640 ˚ A and [O III ] λλ 1661,1666 ˚ A doublet, as well as C III ] λ 1909 ˚ A emission) and both weak and strong optical lines ([O II ] λλ 3727,3730 ˚ A, [Ne III ] λλ 3870,3969 ˚ A, H δλ 4103 ˚ A, blended H γλ 4342 ˚ A and [O III ] λ 4364 ˚ A, He I λ 4472 ˚ A, He II λ 4687 ˚ A, H βλ 4863 ˚ A, [O III ] λλ 4960,5008 ˚ A, He I λ 5877 ˚ A, [O I ] λ 6302 ˚ A, H αλ 6564 ˚ A, [S II ] λλ 6718,6732 ˚ A, and He I λ 7067 ˚ A), affording a plethora of emission-line indicators for metallicity, dust, radiation field, and star-formation diagnostics. Furthermore, despite the low resolution of the prism, strong hints of Fe II absorption can be seen in the rest-frame UV. A clear continuum is also seen at all wavelengths longward of the Lyman and Balmer breaks, indicating the presence of predominantly young, massive stars and longer-lived, lower-mass stars, respectively. The positions of all detected lines are marked by dashed lines (grey for emission lines, black for absorption lines).  \n<!-- image -->  \nEach composite spectrum is then fit with the gsf 5 SED-fitting code (Morishita et al. 2019) to obtain a continuum model. The code fits an assortment of galaxy templates derived from the simple stellar population (SSP) models of Conroy et al. (2009), which are convolved to the resolution of the observed prism spectrum. To ensure a representative set of templates able to fit a large variety of spectral types, we generate the templates using a Chabrier (2003) IMF, a metallicity range of log Z/Z ⊙ = 0 . 01 -1, and stellar ages of 0.001-1.0 Gyrs. For the fit itself, we adopt a nonparametric star formation history and an SMC-like dust attenuation law (Gordon et al. 2003; see also Reddy et al. 2015, 2020), leaving A v as a free parameter with  \nA v = 0 -3. Wavelengths corresponding to emission lines strong enough to potentially influence the fit (i.e., [O II ] λλ 3727,3730 ˚ A, [Ne III ] λλ 3969,3870 ˚ A, H δλ 4103 ˚ A, [O III ] λ 4364 ˚ A, H βλ 4863 ˚ A, [O III ] λλ 4960,5008 ˚ A, He I λ 5877 ˚ A, H αλ 6565 ˚ A, and [S II ] λλ 6718,6733 ˚ A) are masked prior to the fit and are measured separately. IGM absorption is modeled via the prescription in Dijkstra (2014), however we mask wavelengths between Ly α and (rest-frame) 1500 ˚ Agiven the presence of strong line emission and damping wing effects identified by a number of recent studies for sources in our sample (e.g., Tang et al. 2023; Saxena et al. 2023a; Heintz et al. 2023). The continuum spectrum and all gsf -derived parameters are taken from the 50th percentile of the posterior distribution, with the 16th to 84th percentile ranges used for uncertainties.  \nEmission lines are fit using the emcee code (ForemanMackey et al. 2013) after subtraction of the best-fit gsf continuum model, the latter of which also accounts for stellar absorption. Given the large wavelength (and thus resolution) range over which each of the aforementioned lines span, we perform the (stellar and emission line) fits in separate wavelength groups and model each of the lines within a group simultaneously e.g., [O II ] λλ 3727,3730 + Ne III λλ 3969,3870 + H δλ 4103 + [O III ] λ 4364 ˚ A, H βλ 4863 + [O III ] λλ 4960,5008, and H αλ 6565 + [S II ] λλ 6718,6733 ˚ A. For each emission line in a set, we adopt a simple Gaussian profile with a line width (0 . 001 < σ < 100 ˚ A) and derive a posterior distribution for each of the free parameters (central wavelength, amplitude, and line width). All fluxes are derived using the best-fit parameters and their uncertainties, taken as the 50th percentile and semidifference of the 16th-84th percentiles of their posterior distributions, respectively. EW measurements are derived using those fluxes as well as the underlying gsf continuum model, unless stated otherwise.', '4.1. Composite Spectra in Bins of Redshift, Absolute Magnitude, and Stellar Mass': 'We have already provided a first look at the median spectrum of z ⩾ 5 sources from the composite spectrum in Figure 3. The spectrum displays a rich variety of features hinting at a variety of environmental, stellar and nebular effects: a clear Lyman break arising from the absorption of UV photons by an increasingly neutral IGM, a blue rest-frame UV continuum ( β = -2 . 42 ± 0 . 08) dominated by stellar and nebular light from young OB stars, rest-frame optical emission governed primarily by a mix of A- and OB-type stars, and a host of clear line emission from He II λ 1640 ˚ A and [O III ] λλ 1661,1666 ˚ A (all blended), C III ] λλ 1907,1909 ˚ A, [O II ] λλ 3727,3729 ˚ A, [Ne III ] λλ 3869,3968 ˚ A, H γλ 4102 ˚ A, H δλ 4341 ˚ A and [O III ] λ 4264 ˚ A(blended), H βλ 4862 ˚ A, [O III ] λλ 4960,5008 ˚ A, He I λ 5877 ˚ A, [O I ] λ 6302 ˚ A, H αλ 6564 ˚ A and [N II ] λλ 6549,6585 ˚ A (blended), [S II ] λλ 6718,6732 ˚ A, and He I λ 7067 ˚ A (as well as a number of additional marginal lines from, e.g., He I λ 4472 ˚ A and He II λ 4687 ˚ A). Additionally, we tentatively detect absorption features from resonant Fe II λλλλ 2374,2382,2586,2600 ˚ A lines in the rest-frame UV, which in principle allow for determinations of gas outflows along with estimates of gas covering fraction and column density, among other properties (e.g., Martin et al. 2012).  \nWhile remarkably clear, such a spectrum combines galaxies at different epochs, luminosities, and stellar masses, which together yield unphysical results. As  \nsuch, we construct three sets of more physically motivated stacks, one in bins of redshift only, a second in bins of redshift and M UV , and a third in bins of redshift and stellar mass. The redshift ranges for each set are z =[5-6, 6-7, 7-8, 8-9, 9-10, 10-13], while the absolute magnitude bins are -22 . 00 < M UV < -19 . 25 or -18 . 75 < M UV < -16 . 00, and the stellar mass bins are log M ⊙ / M ∗ < 8 . 4 and log M ⊙ / M ∗ > 8 . 6. The limits of M UV and stellar mass are chosen to cover comparable parameter ranges and completeness levels either side of the divide, while leaving a gap to ensure sufficient contrast in average galaxy properties at fixed redshift (the exact size of the gap is limited by the sample size at the time of writing). The resulting composite spectra for the first set of stacks (i.e., redshift only) are shown in their rest frame in Figure 4 and we adopt these as our fiducial set of templates. The median redshifts, absolute magnitudes, and stellar masses for those spectra, as well as the number of objects in each stack, are listed in Table 2.', '4.2. Blue Ultraviolet Continuum Slopes': "Similarly to Figure 3, we find overall blue continuum slopes in each of our high-redshift stacks. Adopting a wavelength window of 1600-2800 ˚ A, masking the region of the C III ] λλ 1907,1909 ˚ A doublet, and binning the continuum by a factor of × 10 (making it less sensitive to flux outliers), we use emcee with a simple powerlaw parameterization (i.e., F λ ∝ λ β ) and measure increasingly blue UV slopes of -2.29 ± 0.01, -2.50 ± 0.10, -2.50 ± 0.04, -2.62 ± 0.02, -2.60 ± 0.01, -2.66 ± 0.02 at redshifts z ≃ 5 . 5, z ≃ 6 . 3, z ≃ 7 . 5, z ≃ 8 . 6, z ≃ 9 . 6, and z ≃ 10 . 7, respectively (from the redshift-binned stacks). The results and their uncertainties (from the 50th and semidifference of the 16-84th percentiles of the posterior distribution, respectively) are shown as a gray line and shading in Figure 5, where a clear and steep decline in UV slope at fixed absolute magnitude is seen for z ≃ 5 -7 sources ( β ≃ -2.1 to -2.5) before a much more gradual evolution at redshifts beyond ( β ≃ -2.5 to -2.6). The trend is suggestive of evolving dust contents and stellar populations, where the lowestredshift sources are characterized by the most dust- and metal-rich interstellar media and aged stellar populations, while higher-redshift sources reflect more pristine ISM conditions and younger stellar populations.  \nA secondary dependence of β on luminosity is well known (see, e.g., Bouwens et al. 2014 and Cullen et al. 2023), and we find that our z -M UV stacks display consistent but systematically offset values compared to the redshift-binned stacks: the most luminous galaxies harbor redder slopes compared to their fainter counterparts  \nTable 2. Spectral Properties of the Composite Spectra, for Star-forming Galaxies Stacked by Redshift, M UV , and Stellar Mass.  \nNote -All quoted uncertainties represent 1 σ values. Those quoted for the median redshift refer to the standard deviation of the values in the stack, those for M UV represent the standard deviation of the spectral fluxes within the selected wavelength window, while all others derive from the semidifference of the 16th and 84th percentiles of the posterior distribution from emcee parameter fitting. All values are corrected for magnification when necessary.  \nFigure 4. Composite NIRSpec spectra (solid lines and filled) and their 1 σ standard deviations (i.e., dashed line) as a function of spectroscopic redshift. Each composite spectrum is shifted to the rest frame using the median spectroscopic redshift (indicated) of the stack and normalized to its λ ≃ 1750 ˚ A flux. All spectra are vertically offset for clarity, and the positions of several rest-frame optical emission lines are indicated by dashed gray lines.  \n<!-- image -->  \n(approximately ∆ β =0.1-0.4 redder, with smaller differences at higher redshifts), likely reflecting more evolved stellar populations and metal- and dust-rich gas reservoirs for the most luminous systems.  \nThe UV-luminous and UV-faint stacks are characterized by trends mirroring the evolution seen in the redshift-binned stacks, while a similar trend is also found with stellar mass, particularly for the most massive systems, which display a virtually identical set of continuum slopes at fixed redshift compared to the UVluminous stacks. The trend is less pronounced for lowermass galaxies compared to their UV-faint counterparts, however, suggestive of stellar mass playing a more dominant role in regulating β with redshift, through the chemical enrichment of the ISM and aging of stellar populations in the most massive galaxies.  \nThe normalization and evolution of the UV slopes with redshift can first be compared to HST-derived results, which generally display redder slopes at comparable redshifts likely attributable to the modest sampling of near-ultraviolet wavelengths compared to the longer baseline afforded by JWST. For instance, for sources at z phot ∼ 5 -7 Wilkins et al. (2011) report median slopes of β ≃ -1 . 9 to -2 . 3, Finkelstein et al. (2012) report median slopes of β ≃ -1 . 9 to -1 . 95 at z phot ∼ 5 -8, Dunlop et al. (2013) report median slopes of β ≃ = -2 . 1 to -1 . 8 at z phot ∼ 7 -9, and Bouwens et al. (2014) report median values of β ≃ -2 . 2 to -1 . 8 at z phot ∼ 5 -9 (for sources at comparable luminosities to here). The z < 8 values are consistent only with the most luminous of our stacks and, interestingly, also point to a steepening with redshift from z ≃ 5 to z ≃ 7 (Wilkins et al. 2011; Finkelstein et al. 2012; Bouwens et al. 2014). However, beyond z ≃ 7 where WFC3 photometry samples only up to λ rest ≃ 1600 ˚ A, the HST photometry lacks the wavelength coverage necessary for accurate constraints, as is evident from the markedly shallower reported slopes (∆ β ≈ 0 . 2 -0 . 7) compared to our spectroscopic results.  \nThe extended baseline afforded by the NIRCam imaging alleviates some of the tension, yielding trends and a range of UV slopes similar to those reported here. For instance, Nanayakkara et al. (2023) report median slopes of β ≃ -2 . 1 to -2 . 4 at z phot ∼ 5 . 1 -6 . 5, while Topping et al. (2023) report median values of β ≃ -2 . 3 to -2 . 5 at z phot ∼ 6 -12, and Cullen et al. (2023) report a median value of β -2 . 32 for sources at z phot ∼ 6 -12 (with comparable luminosities to our stacks).  \nThose NIRCam studies all report trends of declining UV slopes with redshift that are similar to what we report here, with some differences in slope and normalization likely reflecting selection effects and possible sample contamination. For instance, a straight-line fit to the  \nredshift trend produced by our fiducial stacks in Figure 5 yields a slope of dβ/dz = -0 . 06 ± 0 . 01, which is moderately steeper than (but consistent with) a fit to the Topping et al. (2023) median points ( dβ/dz = -0 . 03 ± 0 . 01) over the same redshift range. Conversely, a similar comparison at redshifts of z ≳ 8 with median values from the Cullen et al. (2023) NIRCam points yields slopes of dβ/dz = -0 . 11 ± 0 . 04 and dβ/dz = -0 . 02 ± 0 . 02 for our stacks and their photometric data points, respectively. For the comparisons above, we also find y-axis intercepts of -2 . 04 ± 0 . 12 (these stacks) versus -2 . 06 ± 0 . 05 (fit to the Topping et al. 2023 points) and -2 . 45 ± 0 . 20 (these stacks, z > 8) and -1 . 36 ± 0 . 41 (fit to the derived median Cullen et al. 2023 points).  \nOf particular interest are redshifts beyond which galaxies are deemed to become 'dust-free,' facilitating the transmission of UV photons from massive stars and enhancing their observed UV luminosities perhaps beyond first-order expectations (e.g., Castellano et al. 2022; Naidu et al. 2022a; Ferrara et al. 2023; Mason et al. 2023; Casey et al. 2023). While this analysis and the studies discussed above all associate the bluest slopes to sources at the highest redshifts, none consistently measure dust-free slopes of β ≃ -3 or bluer. A small minority of individual sources within these samples are reported to show slopes bluer than the canonical β ≃ -2 . 7 limit set by nebular continuum contributions (Topping et al. 2023; Cullen et al. 2023), however singlesource measurements from photometry or low SNR spectroscopy are prone to large scatter and uncertainties. Median values consistently measure values redder than or close to this limit and the z > 10 values presented here confirm the lack of β < -3 slopes, at least on average.  \nWhile a direct comparison with photometric results is challenging, the secure redshifts and high SNR measurements of our NIRSpec composite spectra confirm a clear steepening in UV slopes (i.e., toward bluer values) from lower redshifts to z ∼ 7 across the general galaxy population, and a probable flattening out beyond. The measurements point primarily to a declining dust content, possibly arriving close to relatively pristine systems (on average) by redshifts z ≃ 8 -9. While this does not preclude some exceptionally luminous sources to have accelerated their dust buildup compared to the general population (e.g., Ferrara et al. 2022), such examples are rarely seen at redshifts beyond z ≃ 8 -9 and likely represent the exception rather than the rule.  \n4.3. Continuum Breaks Indicative of Evolving Stellar Populations  \n<!-- image -->  \nFigure 5. The UV continuum slopes of z ⩾ 5 galaxies from individual sources (points) and composite JWST/NIRSpec prism spectra, the latter binned as a function of redshift only (gray line and gray shaded uncertainties in both panels), z -M UV (navy and orange points for sources with -22 . 00 < M UV < -19 . 25 mag and -18 . 75 < M UV < -16 . 00 mag, respectively, in the left panel), and z -M ∗ (navy and orange points for sources with log M ∗ /M ⊙ > 8 . 6 and log M ∗ /M ⊙ < 8 . 4, respectively, in the right panel). Individual points (color coded to their M UV and log M ∗ values, with UV-luminous and/or higher-mass sources in bluer colors and the opposite in redder colors) display a large scatter of values (generally between β ∼ -1 . 8 to β ∼ -2 . 8, with a mean 1 σ uncertainty of ≃ 0 . 22) and an overall declining slope with redshift. Composite points confirm the trend, with clear differences between faint/low-mass systems and their bright/high-mass counterparts and a particularly strong redshift evolution with absolute magnitude but shallower evolution with stellar mass.  \n<!-- image -->  \nOne of the primary discontinuities in the rest-frame optical spectra of galaxies are the so-called Balmer and 4000 ˚ A breaks, both of which serve as indirect tracers for the ages of the underlying stellar populations. While both derive from absorption within the atmospheres of stars and are often treated as a single feature, their origins differ - the Balmer break at 3646 ˚ A represents the limit of the Balmer series and occurs in stellar populations younger than 1 Gyr, while the 4000 ˚ A break derives from the absorption of stellar continuum by ionized metals in the atmospheres of stars (including Ca II H and K) and as such increases monotonically with age (Kauffmann et al. 2003), with some dependence also on metallicity.  \nUntil recently, for z ⩾ 5 sources these diagnostics were accessible only through Spitzer photometry (e.g., Hashimoto et al. 2018; Roberts-Borsani et al. 2022; Laporte et al. 2021). However, in near-identical fashion to wide JWST/NIRCam bands, the wide coverage of photometric bands can lead to significant degeneracies (see, e.g., Roberts-Borsani et al. 2022; Desprez et al. 2023). Thanks to the extensive wavelength coverage of the NIRSpec prism, the Balmer break is covered in each of our stacks, allowing for a direct and consistent measurement of the continuum across ∼ 600 Myr of cosmic time. A number of definitions have been utilized for sources  \nin the local Universe, differing primarily in the central wavelengths and widths of the pseudobandpasses where the red and blue portions of the continuum are measured. Given the smearing of emission-line fluxes across the stellar continuum, traditional definitions using 40504250 ˚ A and 3750-3950 ˚ A (Bruzual A. 1983) or narrower 3850-3950 ˚ A and 4000-4100 ˚ A (Balogh et al. 1999) windows would result in emission-line-contaminated measurements. For a consistent and clean set of measurements, therefore, we take care to define spectral windows available across our entire redshift range that are free of line emission. We adopt wavelength ranges of λ rest = 3500 -3630 ˚ A and 4160 -4290 ˚ A and, to mitigate the impact of noise or lingering features, determine the average fluxes from straight-line fits to the fluxes in each of those windows. We refer to the measurements as Balmer 'indices' and plot the results for our redshift-binned stacks in Figure 6, where uncertainties are defined as the 1 σ error on the best-fit normalization parameter.  \nWe find a steadily declining Balmer index with redshift. The largest index is found in our z ∼ 5 . 5 stack, with F ν, 4225 /F ν, 3565 = 1 . 05 ± 0 . 06, which rapidly shifts to F ν, 4225 /F ν, 3565 < 1 by z ∼ 6 . 5 and declines further with redshift to 0.75 ± 0.26 by z ≃ 9 . 6. For our highest-redshift stack at z ≃ 10 . 7, the red continuum is  \ntoo noisy to determine a reliable break, while for our z -M UV and z -M ∗ stacks the SNR is insufficient to verify whether differences between subsamples exist (however the general decline toward reduced values with redshift remains).  \nThe small values we observe here are unlike those found for the bulk of the general population at z ∼ 0, which show indices of ∼ 1 -2 and include effects from both the Balmer and 4000 ˚ A breaks (Kauffmann et al. 2003). For those sources, even the youngest galaxies are limited to values of ∼ 1. Local analogs, however, typically associated with extremely metal-poor ISM, young stellar ages, and enhanced nebular continuum (all of which are conducive to Balmer jumps; Mingozzi et al. 2022; Inoue 2011), show greater consistency with the results shown here (e.g., Izotov et al. 2021). At z > 5, only a handful of sources have been studied with spectroscopic measurements, and these show significant scatter in the reported indices (e.g., Looser et al. 2023, Cameron et al. 2023a, Vikaeus et al. 2024, Witten et al. 2024), likely reflective of stochastic star formation histories.  \nThe trend observed here is reminiscent of the declining UV continuum slopes found in the previous section, and is also matched by decreasing mass-weighted stellar ages inferred from our gsf continuum fitting. Those ages begin at values of ∼ 95 Myr for our z ∼ 5 . 5 stack, decrease to ∼ 18 Myr at z ∼ 6 . 5, and steady out to ages of < 10 Myr at redshifts beyond. The inferred breaks and ages are in line with photometric expectations (e.g., Stefanon et al. 2023; Whitler et al. 2023a) and recent spectroscopic measurements of some individual sources (Vikaeus et al. 2024; Cameron et al. 2023a), which all require young stellar ages, low quantities of metals and dust, and significant contributions from the nebular continuum to explain the observations.  \nTo verify whether our observations agree with expectations from theory and simulations, we overplot in Figure 6 the expected break index ranges measured in identical fashion on the spectral models of Schaerer (2003), for a variety of metallicities (0 . 01 -0 . 2 Z ⊙ ) and stellar ages (5-20 Myr), all of which include contributions from the nebular continuum. The models are in excellent agreement with the values measured here. We also overplot (as navy lines) the median values from the SPHINX 20 (Katz et al. 2023) and FLARES (Lovell et al. 2020; Vijayan et al. 2020) hydrodynamical simulations, limiting the comparison to sources with metallicities of < 0 . 01 Z ⊙ , stellar masses of < 10 8 . 0 M ⊙ , and stellar ages comparable to those inferred by our gsf modeling at each redshift interval. For the former, we measure the indices in identical fashion using the simulated spectra, while the later are derived directly from the FLARES  \ncode output using the formalism of Wilkins et al. (2024). Both simulations are able to match the range of our measured values, and moreover show excellent agreement with the trend toward reduced values at higher redshift (see Vikaeus et al. 2024 and their comparison to the FLARES and DELPHI simulations with a smaller number of confirmed JWST sources). Similar conclusions can be derived when comparing with the model of Ferrara et al. (2023), which also shows agreement with the range of values we probe and a decline in the Balmer break index with redshift.  \nThe strength of the break is sensitive to a number of physical properties, including stellar age, dust attenuation, metallicity, the escape fraction of ionizing Lyman continuum photons, and the shape of the star formation history (Wilkins et al. 2024). Determining the precise contributions of each to the breaks measured here is beyond the scope of the present work, however we note that the overall blue UV slopes seen in Section 4.2 disfavor significant dust contributions and the consistency between the gsf ages and the age-limited simulations suggest stellar age is the main driver.", '4.4. Prominent C III ] Emission in Early Galaxies': "C III ] λλ 1907,1909 ˚ A emission is a useful indicator for the metal enrichment and electron density of the ISM (e.g., Jaskot & Ravindranath 2016; D'Eugenio et al. 2023), the hardness of the ionizing spectrum (e.g., Stark et al. 2017), and emission-line diagnostics to distinguish between star formation and AGN activity (e.g., Laporte et al. 2017a; Nakajima et al. 2018). Characterizing the presence and evolution of C III ] thus represents an important step toward understanding the ISM conditions of galaxies out to the earliest times. Our composite spectra reveal significant C III ] emission at all redshifts, with measured rest-frame EWs (EW 0 ) of 4.6 ± 0.7 ˚ A, 7.4 ± 1.6 ˚ A, 7.8 ± 0.9 ˚ A, 11.2 ± 1.2 ˚ A, 12.8 ± 1.0 ˚ A, and 13.7 ± 0.9 ˚ A for the redshift-binned stacks at z ≃ 5 . 5, z ≃ 6 . 3, z ≃ 7 . 5, z ≃ 8 . 6, z ≃ 9 . 6, and z ≃ 10 . 7, respectively, marking an increase of a factor of ∼ 3 from z ∼ 5 to z ∼ 11 and a median value of 9.5 ± 0.7 across the full sample. The values are plotted in Figure 7 and showcase a very clear trend of increasing emission with redshift, with an evolution best described by a straight-line fit of EW 0 (C III ])=1.84( ± 0.19) z -5.38( ± 1.53) (indicated by a dashed line and 1 σ uncertainty fill). Furthermore, we find no clear evidence for varying C III ] strength in our z -M UV stacks, but do find tentative evidence for some enhanced emission in the least massive sources of our z -M ∗ stacks compared to their higher-mass counterparts.  \nFigure 6. The Balmer break as measured directly from our fiducial composite spectra ( F ν 4225 /F ν, 3565 ), as a function of redshift (red squares and uncertainties). For comparison, the expected range of values for low-metallicity (0.01-0.2 Z ⊙ ) and young (5-20 Myr) sources Schaerer (2003) are shown as the gray region. Median values for comparable sources in the SPHINX 20 (Katz et al. 2023) and FLARES (Wilkins et al. 2024) simulations are also shown as solid and dashed navy lines, respectively, and the model of Ferrara et al. (2023) as a dotted line. All comparisons show agreement with the range of values measured here, as well as the observed trend with redshift.  \n<!-- image -->  \nThe implied prevalence of C III ] and the observed correlation are consistent with extrapolations of measurements from lower redshift. High resolution spectroscopy over statistical samples of z ≃ 2 -4 sources find a high incidence of C III ] emission in both individual sources  \n(e.g., ∼ 24-29%; Le F'evre et al. 2019; Llerena et al. 2022; Stark et al. 2014) and stacked spectra (e.g., Shapley et al. 2003; Amor'ın et al. 2017; Le F'evre et al. 2019; Llerena et al. 2022). Reported EW 0 s for typical sources generally range between a few to 10 ˚ A although most studies also report a small tail in the distribution toward larger values, generally attributed to extreme radiation fields from very metal-poor galaxies or AGN (Jaskot & Ravindranath 2016; Nakajima et al. 2018); for instance, for a statistical sample of 3899 (Le F'evre et al. 2019, see also Erb et al. 2010; Stark et al. 2014; Amor'ın et al. 2017; Llerena et al. 2022) report EWs of 3-10 ˚ A for 20% of their sample and find higher values of > 10 ˚ A for only ∼ 4% of the sample. The same authors find median values of 2.0 ˚ A and 2.2 ˚ A for sources at 2 < z < 3 and 3 < z < 4, respectively, while Shapley et al. (2003) find a marginally smaller but consistent value of 1.7 ˚ A in stacked spectra of UV-selected z ∼ 3 LBGs. Each of these are in good agreement with the trend we observe in Figure 7, which if extrapolated suggests EWs of 0.12.0 ˚ A at z = 3 -4. At even lower redshifts of z ∼ 1, Du et al. (2017) find a detection rate of ∼ 22% and a median EW of 1.3 ˚ A, consistent with the decline seen in Figure 7. The ubiquity of our (stacked) detections suggests C III ] is likely to be at least as prominent in z ⩾ 5 sources as in z ∼ 1 -4 sources, while a well-established link between C III ] EW and metallicity (Jaskot & Ravindranath 2016; Nakajima et al. 2018; Llerena et al. 2022) suggests the latter is the main driver of our trend.  \nFor comparisons at high redshift, in Figure 7 we also plot a number of reported values for known C III ] emitters in the Epoch of Reionization (Stark et al. 2015, 2017; Bunker et al. 2023a; Larson et al. 2023; D'Eugenio et al. 2023; Castellano et al. 2024). We find a large scatter of individual objects around the relation (albeit with significant uncertainties), where a number of objects (A383-5.2, EGS-zs8-1, GS-z12, and GHZ2) fall above the relation and by this metric alone could indicate extreme ISM conditions - e.g., EGS-zs8-1 is a known z = 7 . 7 LAE with extreme [O III ]+H β line emission (Oesch et al. 2015; Roberts-Borsani et al. 2016), while GHZ2 displays remarkably luminous rest-frame UV lines possibly indicative of AGN activity (Castellano et al. 2024). For the sources instead lying close to the range probed by the general highz population, these include apparently extreme sources such as GN-z11 and CEERS 1019, two Ly α -emitting sources with atypical UV line emission and nitrogen enhancement (Bunker et al. 2023a; Marques-Chaves et al. 2024; Senchyna et al. 2024). In the case of GN-z11, we note the source exhibits a measured EW 0 (C III ]) of 12.5 ± 1.1 ˚ A (Bunker et al. 2023a), compared to a predicted value of 14.1 ± 2.5  \n˚ A. Our results imply the classification of 'extreme' systems for these particular objects are not reflected in their C III ] emission, based on a comparison to the range of values characteristic of the general galaxy population.  \nAt even higher redshifts of z ∼ 14, JADES-GS-z140 appears well below the extrapolation of the trend, suggestive of an exceptionally low metallicity (Carniani et al. 2024). Should other sources at similarly high redshifts reveal comparable C III ] EWs, this might imply a turnover in the trend somewhere around z ∼ 13 which is not captured by the extrapolation of the evolution from z < 11, where sources begin to reveal virtually pristine ISM conditions.  \nLastly, although the presence of C III ] in our stacks is unambiguous, the presence or otherwise of additional rest-frame UV lines is less clear. For instance, although we find detections of He II +[O III ] rest-frame UV line emission in Figure 3 and Figure 4, those detections are lower SNR and do not translate to the smaller samples in any of our z -M UV or z -M ∗ stacks. Similar conclusions can be derived for other rest-frame UV lines such as N V λ 1243 ˚ A, N IV λ 1486 ˚ A, and C IV λλ 1548,1550 ˚ A, all of which are seen in individual Lyman break or C III ]-emitting galaxies at z ∼ 2 -3 (e.g., Shapley et al. 2003; Stark et al. 2014) and out to higher redshifts of z > 6 (e.g., Bunker et al. 2023a, Topping et al. 2024), but not in our stacks.", '4.5. Rest-frame Optical Line Emission': "By far the most common features found in our stacked spectra are signatures of strong rest-frame optical line emission at virtually all redshifts, UV magnitudes, and/or stellar masses, highlighting the nature by which most of our z > 5 sample were selected. These include clear signatures of line emission from [O II ] λλ 3727,3729 ˚ A, [Ne III ] λ 3869 ˚ A, [Ne III ] λ 3968 ˚ A, H δλ 4102 ˚ A, H γλ 4341 ˚ A+[O III ] λ 4364 ˚ A, H βλ 4862 ˚ A, [O III ] λλ 4960,5008 ˚ A, He I λ 5877 ˚ A, H αλ 6564 ˚ A+[N II ] λ 6584 ˚ A, and [S II ] λλ 6718,6732 ˚ A. Due to the wavelength resolution and limits of the NIRSpec prism, most line doublets and closely spaced pairs (e.g., H γ +[O III ]) are blended, while emission from He I , H α , and [S II ] are detectable only out to z ∼ 7 -7 . 5 and strong lines of [O III ]+H β out to z ≃ 9 . 6.  \nThe weakest lines across the λ rest ≃ 3750-4500 ˚ A wavelength range are generally the [Ne III ] λ 3968 ˚ A and H δ lines, with rest-frame equivalent widths in the range EW 0 ≃ 11-30 ˚ A for the former and EW 0 ≃ 16-22 ˚ A for the latter. Next in strength are the [O II ] λλ 3727,3729 ˚ A doublet and H γ +[O III ], with EW 0 ≃ 14-39 ˚ A and EW 0 ≃ 48-76 ˚ A, respectively.  \nFigure 7. The rest-frame equivalent width of C III ] λλ 1907,1909 ˚ A line emission as a function of spectroscopic redshift, as measured from our redshift-binned composite spectra. The trend of increasing EW with redshift is well parametrized by a straight line fit (dashed purple line with pink uncertainties). Measurements from the literature (Stark et al. 2015, 2017; Larson et al. 2023; Bunker et al. 2023a; D'Eugenio et al. 2023; Castellano et al. 2024; Carniani et al. 2024) are added as grey points, for comparison.  \n<!-- image -->  \nProminent [O III ]+H β lines are found across all redshifts (except the highest redshift of z ∼ 10 . 66 where the lines fall beyond the red limit of the prism) with combined rest-frame EW 0 s of 752 ± 10 ˚ A, 1031 ± 31 ˚ A, 1292 ± 68 ˚ A, 1214 ± 54 ˚ A, and 910 ± 18 ˚ A for redshifts z ≃ 5 . 5, z ≃ 6 . 3, z ≃ 7 . 5, z ≃ 8 . 6, and z ≃ 9 . 6. An increase of ∼ × 1 . 7 is found between the lowest and highest values, which reflect some possible redshift evolution at least between z ∼ 5 . 5 -8 . 5. The range of values reported here illustrate the strong line nature of the underlying sample.  \nThe median value of the sample is 1031 ± 31 ˚ A which is in excellent agreement with values reported in the literature from JWST photometry at comparable redshifts (e.g., Endsley et al. 2021b). Similarly to observations of the C III ] line, we find no significant evidence for stronger line emission in UV-faint galaxies compared to their UV-bright counterparts, but do find evidence for such a trend with stellar mass. This is comparable to the trend found by Reddy et al. (2018), who show a clear decrease in virtually all rest-frame optical line EWs at higher stellar masses. Prominent He I λ 5877 ˚ A and H α emission are also seen in the lowest-redshift stacks (at redshifts of z ≃ 5 . 5 -7 . 5 for the former and z ≃ 5 . 5 -6 . 3 for the latter) before moving beyond the prism cover-  \nFigure 8. The distribution of EW 0 s for the rest-frame optical emission lines observed in our redshift-stacked spectra, as measured from the Gaussian fits described in Section 3.7.  \n<!-- image -->  \nr EWs are reported to be 16-32 ˚ A for He I at redshifts z ≃ 5 . 5 -7 . 5, and 559-922 ˚ A for H α at z ≃ 5 . 5 -6 . 3. The EWs of each of the rest-optical lines are shown in Figure 8, and we list these along with their H β -normalized fluxes in Table 5 in the Appendix.", '5.1. Dust Attenuation from Spectral Energy Distribution Fitting and Balmer Decrements': "Constraints on the prevalence and evolution of dust at high redshift have yielded mixed results. Until recently, constraints in Reionization Era galaxies have come from SED fitting of (primarily) rest-frame UV photometry (e.g., Strait et al. 2020; Roberts-Borsani et al. 2022), or NIR continuum measurements from Atacama Large Millimeter/submillimeter Array (ALMA; e.g., Laporte et al. 2017b; Tamura et al. 2019; Bakx et al. 2021; Bouwens et al. 2021; Ferrara et al. 2022). Although the vast majority of objects have yielded apparently low to negligible dust contents, some marked exceptions from ALMA have challenged expectations with large dust quantities by z ∼ 7 -8 (Laporte et al. 2017b; Ferrara et al. 2022). More direct constraints of nebular dust attenuation with JWST spectroscopy of the Balmer decrement (i.e., the ratio of Balmer lines to H β ) have indicated ISM approaching 'dust-free' states by z ≃ 6, with a clear correlation of enhanced color excesses ( E ( B -V ) neb ) with stellar mass (above log M ∗ /M ⊙ ≳ 9 . 5) that appears invariant at z ≃ 2 -6 (Shapley et al. 2022, 2023). Importantly, however, the scatter and uncertainties on both  \nindividual measurements and small sample sizes remain large, making general conclusions challenging.  \nWe evaluate the potential contributions of dust to the attenuation of our composite spectra, via both continuum and nebular measurements. The evolution of the UV continuum slope presented in Section 4.2 already provides strong indication of declining dust obscuration toward higher redshift. This interpretation is mirrored by our SED-fitting results, which, assuming the SMClike dust curve of Gordon et al. (2003) and assuming R V = 2 . 74, yield negligible E ( B -V ) stel values of 0.030.10 and a clear decline in values from lower to higher redshift.  \nA consistent assessment of the nebular color excess with redshift is more challenging: H α is covered only in our two lowest-redshift stacks and blended with [N II ] λ 6585 ˚ A, H γ is blended with [O III ] λ 4364 ˚ A, and the faintness of H δ makes its measurements extremely sensitive to small fluctuations in the underlying continuum fit. Assuming the Milky Way extinction law of Cardelli et al. (1989, see Reddy et al. 2020 for its appropriateness to highz nebular spectra), we derive upper limits on the H α /H β decrement (intrinsic ratio H α /H β =2.860) of 2.65 and 3.38 for our z ≃ 5 . 5 and z ≃ 6 . 3 spectra, which are consistent with no and low ( E ( B -V ) neb ≃ 0 . 17) dust scenarios, respectively. Given those values are derived from our lowest-redshift stacks, they can be considered upper limits on the dust attenuation for our higher-redshift stacks as well. The low E ( B -V ) stel values inferred from our SED fitting and lowest-redshift spectra suggest little to no dust in our stacks, and thus we do not make additional corrections to our spectra in forthcoming sections.", '5.2. Metallicities and Ionization States': "The plethora of emission lines in our stacks allows us to constrain the chemical enrichment of galaxies to high significance and out to the earliest times. In theory, the presence of the [O III ] λ 4364 ˚ A and [O III ] λ 5008 ˚ A lines in our spectra also enables us to determine direct T e -based oxygen abundances out to z ∼ 9 . 7. However, the prism's spectral resolution blends the [O III ] λ 4364 ˚ A line with H γ , preventing us from modeling it with high confidence (see Hu et al. 2024). We thus derive gas-phase oxygen abundances using the JWST-derived calibrations in Sanders et al. (2023b), specifically using the so-called 'O32', 'R23', and 'Ne3O2' ratios which yield metallicities of 12+log (O/H) ≈ 7.3-7.9, 7.1-7.5, and 7.3-7.6, respectively (approximately 0.03-0.16 solar, assuming 12+log(O/H) ⊙ =8.69; Asplund et al. 2009). The most obvious trend comes from the O 32 calibration, which at face value appears to show a well-defined but small evo-  \nreasing metallicity with redshift that is in agreement with single-object and average derivations in the literature (Sanders et al. 2023b; Curti et al. 2023; Nakajima et al. 2023; Hu et al. 2024). We verify the significance of the evolution by comparing the difference between stacked points and the uncertainties on those differences. We find the decrease between points to be larger than the associated uncertainties by 2-5 σ (except for the decrease between the z ∼ 7 . 5 and z ∼ 8 . 6 points), adding statistical robustness to the trend. Despite the agreement, suggestive of some evolution in the massmetallicity relation within the redshift range probed, we consider whether moderate stellar mass differences between our stacks could induce some artificial evolution if one assumes an invariant mass-metallicity relation across those redshifts. In such a case, based on the stellar mass differences between our stacks alone, we would infer a metallicity difference of only 0.09 dex (Nakajima et al. 2023). This is significantly smaller than the metallicity difference we observe between our z ∼ 5 . 5 and z ∼ 9 . 5 stacks, which amounts to 0.58 dex. Furthermore, were the metallicity evolution we observe driven by differences in the stellar masses sampled in each redshift stack, the redshift-metallicity trend would then have to mirror the redshiftM ∗ evolution between our stacks, which it does not - i.e., the metallicity would have to scatter in the same fashion as the stellar mass between our stacks. As such, while the uncertainties mean we cannot rule out some modest influence by stellar mass variations, the overall trend of decreasing metallicities with redshift appears most consistent with some small evolution in the mass-metallicity relation between z ∼ 5 and 10.  \nIonization versus excitation diagrams also provide a means of interpreting the underlying conditions of the ISM. We show two variations of such a diagram in Figure 10, in similar fashion to Figure 8 of Bunker et al. (2023a). In the left panel we show the R 23 -O 32 line ratios, where R 23 is used as a proxy for gas-phase metallicity while O 32 is invoked as a proxy for ionization parameter (e.g., Kewley et al. 2019). Low-redshift ( z < 0 . 3) points from the Sloan Digital Sky Survey (SDSS) MPAJHU catalogs 6 (Aihara et al. 2011) are added to the plots (in black), for reference. A clear trend appears from these diagnostics, whereby the lowest R 23 and O 32 values are found in local galaxies and the highest values in the highest-redshift sources, indicating an evolution towards more extreme and ionized ISM with decreased gas-phase metallicities at higher redshift, i.e., the harder  \nFigure 9. Gas-phase oxygen metallicities as a function of redshift for the stacks derived in this work (blue squares), as estimated from the strong line O 32 prescription of Sanders et al. (2023b). Literature points from Sanders et al. (2023b), Curti et al. (2023), and Nakajima et al. (2023) are plotted as gray symbols. The dashed lines mark levels of constant metallicity in solar units, assuming 12+log (O/H) ⊙ =8.69 (Asplund et al. 2009). A clear trend of decreasing metallicity toward higher redshifts is evident, likely reflecting larger H I gas fractions of highz sources.  \n<!-- image -->  \nionizing spectrum (traced by an ionization parameter, U ) associated with highz sources drags the R 23 ratio to higher values at fixed O 32 (Steidel et al. 2016). Our points, shown as color-coded stars, are in excellent agreement with the single-source measurements at comparable redshifts by Mascia et al. (2023a), Cameron et al. (2023b), Sanders et al. (2023a), and Nakajima et al. (2023).  \nThe evolution is seen not only between the locus of SDSS sources and locus of stacked (and literature) highz points, but possibly between the stacked points themselves, where the spectrum at z ∼ 5 . 5 displays the lowest values of R 23 = 6 . 3 ± 0 . 2 and O 32 = 6 . 8 ± 0 . 1 and increases to R 23 = 8 . 4 ± 1 . 1 and O 32 = 13 . 0 ± 1 . 3 by z ∼ 8 . 6 (our redshift z ∼ 9 . 5 point appears to be somewhat scattered off the main locus of points with R 23 = 5 . 5 ± 0 . 9 and O 32 = 31 . 1 ± 5 . 3, possibly due to the comparatively small numbers in the stack). Despite the limited range in parameter space probed by the stacks, the z ∼ 5 . 5 -8 . 6 spectra hint at nonnegligible evolution in the ISM conditions of star-forming galaxies within the first ∼ 1 Gyr of the Universe (or ∼ 500 Myr within the redshift range). For comparison, we plot the updated line ratios for star-forming sources  \nfrom Gutkin et al. (2016), as a function of metallicity and ionization parameter, log U . Our points are consistent, as expected, with subsolar metallicity tracks and are characterized by large log U values ( -2 . 5 to -1 . 5) compared to z ∼ 0 sources. Lastly, the two locuses of low- and high-redshift sources are also well connected by measurements at intermediate redshifts, namely those of Sanders et al. (2023a) at z ∼ 3 . 0 -7 . 5 and Nakajima et al. (2023) at z ∼ 4 -10, which sit firmly between the SDSS points and ours.  \nAnalogous to the R 23 -O 32 diagram is the Ne 3 O 2 line ratio paired with ([Ne III ] λ 3869 ˚ A+[O II ] λλ 3726,3729 ˚ A) / H δ , neither of which make use of strong [O III ] and H β lines and thus allow for an inspection of ISM conditions with NIRSpec in principle out to z ∼ 13 . 2. In this case the ([Ne III ]+[O II ])/H δ ratio traces the gas-phase metallicity of galaxies while the Ne 3 O 2 ratio traces their ionization strengths. We plot our values and a number of literature values (Sanders et al. 2023a; Nakajima et al. 2023; Cameron et al. 2023b) in the right panel of Figure 10. Similarly to above, we also plot the MPA-JHU values for z < 0 . 3 SDSS star-forming galaxies, as well as the photoionization models of Gutkin et al. (2016). For reference, we also show the reported values for GNz11 (Bunker et al. 2023a). Again, we find a clear separation between sources at low redshift, which dominate regions of parameter space associated with lowionization and high-metallicity systems, and high redshift, which shift toward high values of Ne 3 O 2 and low values of ([Ne III ]+[O II ])/H δ where high-ionization and low-metallicity sources dominate. Interestingly, the ratios presented by Bunker et al. (2023a) for GN-z11 are in excellent agreement with those found for our z ∼ 9 . 5 stack, reflecting ISM conditions comparable to those of 'typical' sources at higher redshift and a conclusion those authors found when comparing to the individual measurements in z ∼ 5 . 5 -9 . 5 sources by Cameron et al. (2023b). A similar conclusion can be derived for GHZ2, whose Ne 3 O 2 value (see Calabro et al. 2024) places it at the extreme end of the distribution but within the locus of high-redshift sources clearly distinct from the ISM conditions of z ∼ 0 galaxies. Given its proximity to the measurements of other individual sources at high redshift (e.g., Sanders et al. 2023b; Nakajima et al. 2023) and our stacked spectra, the result suggests GHZ2 may too be reflective of somewhat typical ISM conditions above z > 5.  \nTo conclude, the modest change in the O 32 and R 23 values from the stacked and literature results shown in the left panel of Figure 10 point to some possible (albeit modest) evolution in the ionization state and gas-phase metallicities from redshifts z ≃ 9 -10 to z ≃ 3 -5.  \nThe combination of both panels clearly indicate a far more drastic change in the ionization state and excitation conditions of galaxies from redshifts z ≃ 3 -5 to z ∼ 0, suggesting the bulk of metal and dust buildup in galaxies occurs between the latter two populations.", '5.3. Do Extreme Properties Lead To Extreme Ly α ?': "A significant number of individual sources in our starforming sample display strong ( ≳ 30 ˚ A) Ly α emission. Characterizing the properties of such populations and determining whether they exhibit enhanced intrinsic ionizing capabilities compared to their Ly α -attenuated counterparts represents a key goal of JWST. Moreover, a thorough understanding of their intrinsic properties is a prerequisite toward using them as tracers of Cosmic Reionization (e.g., through the opacity of Ly α ; Mason et al. 2019a).  \nWe identify a sample of 54 LAEs (12% of our fiducial sample) via a combination of previous works (specifically Tang et al. 2023, Saxena et al. 2023a, Saxena et al. 2023b, and Napolitano et al. 2024) and visual inspection of our compiled spectra. For reference, 10/54 of the LAEs ( ∼ 17%) fall in the Abell 2744 field, 16/54 ( ∼ 32%) fall in the CANDELS-EGS field, 1/54 ( ∼ 2%) fall in the GOODS-North field, 23/54 ( ∼ 42%) fall in the GOODS-South field, 3/54 ( ∼ 5%) fall in the MACS 0647 field, and 1/54 ( ∼ 2%) fall in the RXJ2129 field. While the main driver of Ly α visibility in each of these fields remains subject to debate, the variability is likely linked to a mixture of intrinsic galaxy properties, large-scale ionized bubbles, and cosmic variance. While the sample of LAEs used here is not sufficiently large to perform stacks in bins of global properties, it is large enough for a comparison to a stacked sample of redshift- and M UV -matched non-LAEs.  \nTo separate the two populations, we first derive Ly α EWs for each of our LAEs using a Gaussian profile and the median prism continuum flux within a wavelength range λ rest ≃ 1250 -1750 ˚ A. A straight-line fit is included in the Gaussian fit to account for underlying continuum emission (when present). The resulting range of values span EW 0 , Ly α ≃ 11-919 ˚ A (with a median value EW 0 , Ly α ≃ 155 ˚ A), consistent with ranges derived by previous works. For all other objects, where the line remains undetected in the prism spectra, we derive 2 σ upper limits using instead the flux measured within a 200 ˚ A window centered at the redshifted position of Ly α . The resulting distribution of values, along with the redshifts of their respective hosts, is shown in Figure 11 where we highlight in maroon our LAE sample and in gray the upper limits, which unsurprisingly span a lower range of values, EW 0 , non -LAE ≃ 1-86 ˚ A.  \n<!-- image -->  \nFigure 10. The R 23 -O 32 (left) and Ne 3 O 2 -([Ne III ]+[O II ])/H δ (right) diagrams tracing the ionization and excitation states of galaxy ISM across cosmic time. In each plot, the values and uncertainties (often smaller than the symbols) for our stacked sources are shown as colored stars (color-coded according to Figure 4), while values from the literature at comparable redshifts (Sanders et al. 2023a; Mascia et al. 2023a; Cameron et al. 2023b; Nakajima et al. 2023) are shown as gray points with varying marker symbols. GN-z11 is shown as a green diamond in the right panel using values from Bunker et al. (2023a), while GHZ2 is shown as a green square using values from Calabro et al. (2024). For reference, the locus of z < 0 . 3 points from SDSS are shown as a density plot. Tracks from the photoionization models of Gutkin et al. (2016) are shown in blue for a variety of metallicities and ionization parameters, where the shade of blue denotes the metallicity (from 0.01 to 0.4 Z ⊙ going from light shaded to dark shaded) and the size of the connecting points describes the ionization parameter (from log U = -3 to -1 in -0 . 5 intervals, with the largest values represented by the largest circles).  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 11. Left: the distribution of Ly α rest-frame equivalent widths as a function of spectroscopic redshift for the objects in our fiducial sample. Objects with confirmed Ly α emission from prism observations are highlighted as maroon circles, while undetected sources are represented by 2 σ gray upper limits. Dashed lines represent the EW limits used for the construction of the composite spectra. Middle: the same EW measurements from the left panel for Ly α -detected sources only, as a function of absolute magnitude. Similar measurements from Oesch et al. (2015), Tang et al. (2023), and Saxena et al. (2023a) are plotted for comparison. Right: the same as the middle plot but as a function of O 32 parameter. Results for both low- and high-redshift studies using measurements from Flury et al. (2022), Tang et al. (2021, 2023), and Saxena et al. (2023a) are plotted for comparison.  \n<!-- image -->  \nFor comparison, we also include in the middle and right plots of Figure 11 the Ly α EWs for the linedetected z ≃ 6 -9 sources in Oesch et al. (2015), Tang et al. (2023), and Saxena et al. (2023a) and examine potential trends. Firstly, in the middle plot of Figure 11 we see a clear trend of increasing EW with decreasing luminosity. This is expected for UV-selected sources (see Schaerer et al. 2011; Schenker et al. 2014; De Barros et al. 2017; Stark et al. 2017) where Ly α is (in general) more easily attenuated in luminous sources dominated by enhanced H I column densities and dust obscuration. However, important selection effects apply: our measurements naturally probe the strongest Ly α EWs (which correspond to moderately and extremely faint sources with M UV > -20 mag) and R ≳ 1000 spectroscopy from the literature probe the weaker emission from more luminous sources, however weak Ly α (EW 0 , Ly α ≲ 10 ˚ A) from the intrinsically faint populations ( M UV > -18 mag) probed by flux-limited samples are beyond the sensitivity of the data sets considered here.  \nSecondly, we verify whether at fixed O 32 ratio LAEs at higher redshifts display lower EWs. Although the dynamical range where such comparisons are viable is limited, we do not find any significant evidence for decreasing Ly α EW at fixed O 32 ratio across the combined samples plotted in the figure (i.e., this work; Tang et al. 2021, Tang et al. 2023, Saxena et al. 2023a, and Flury et al. 2022), which together span a redshift range z ∼ 0 . 3 -9 . 0. We caution, however, that such an interpretation is sensitive to selection effects, where our prism sample probes only the strongest LAEs and the higher resolution spectroscopy (sensitive to lower EW measurements) is limited in sample size. An enlarged sample of R > 1000 measurements at z > 6 is needed for a more complete comparison with which to derive conclusive results.  \nFor our composite spectra we perform two stacks, one of all sources with EW 0 , Ly α > 30 ˚ A (which corresponds to the approximate lower limit seen in the LAEs in Figure 11) and z spec < 8, and one with EW 0 , Ly α < 10 ˚ A and z spec < 8. In total, this results in 53 sources for the LAE stack and 249 sources for the non-LAE stack. The composites are an excellent match in redshift ( z spec , LAE = 5 . 92 ± 0 . 60 and z spec , non -LAE = 5 . 92 ± 0 . 74, respectively) and M UV ( M UV , LAE = -18 . 50 ± 0 . 08 and M UV , non -LAE = -18 . 70 ± 0 . 03, respectively). Furthermore, while approximately 75% of compiled LAEs fall in either the GOODS-North or CANDELS-EGS field (where the presence of large-scale ionized bubbles are suspected; Tilvi et al. 2020; Larson et al. 2022; Endsley &Stark 2022), the sampling of sources from six indepen-  \nlds should, to a degree, mitigate cosmic variance biases (see Roberts-Borsani et al. 2024 for examples of spectroscopy over independent sight lines). The combination of the matched redshifts and M UV 's, as well as the sampling from independent fields, allows for a close comparison of the populations' intrinsic properties independently of clustering effects or IGM opacity (see Roberts-Borsani et al. 2023a).  \nWe note modest differences in the UV continuum slopes of the two stacks, with β = -2 . 59 ± 0 . 05 for the LAE spectrum and β = -2 . 40 ± 0 . 02 for the nonLAE spectrum, both within the range of values probed by the stacks in Figure 5 and suggestive of dust-poor systems. Remarkably, we find striking differences in the presence and intensity of rest-frame UV line emission. Our LAE spectrum displays a host of line emission from C IV λλ 1548,1550 ˚ A (EW 0 = 2 . 4 ± 0 . 1 ˚ A), He II λ 1640 + [O III ] λλ 1660,1666 ˚ A (EW 0 = 3 . 8 ± 0 . 1 ˚ A), and C III ] λλ 1907,1909 ˚ A (EW 0 = 6 . 9 ± 0 . 1 ˚ A), while the non-LAE spectrum displays only reduced C III ] λλ 1907,1909 ˚ A emission (EW 0 = 4 . 5 ± 0 . 1 ˚ A) and a hint of He II λ 1640+[O III ] λλ 1660,1666 ˚ A(EW 0 < 1 . 2 ˚ A at 2 σ ), without any obvious C IV emission (EW 0 < 1 . 8 ˚ A at 2 σ ). The C III ] EWs of the two stacks are consistent with the upper and lower limits of the relation presented in Figure 7 for their respective redshifts, while the enhanced C III ] emission in the LAE composite relative to the non-LAE composite is not unexpected given results from individual and stacked spectra at z ∼ 2 -4 almost always show larger C III ] EWs with larger Ly α EWs (Shapley et al. 2003; Amor'ın et al. 2017; Stark et al. 2017; Marchi et al. 2018).  \nThe presence of C IV in our LAE spectrum alone is also suggestive of particularly extreme ionizing properties, requiring 48 eV to convert C III ] to C IV and most typically seen in extremely metal-poor systems or AGN (Feltre et al. 2016; Senchyna et al. 2019). C IV has already been identified in a number of objects in the Reionization Era (e.g., Bunker et al. 2023a; Topping et al. 2024) and their lower-redshift analogs (e.g., Erb et al. 2010; Stark et al. 2014; Rigby et al. 2015; Berg et al. 2019; Du et al. 2020), with a significant fraction of those also displaying Ly α emission (e.g., Schenker et al. 2012; Schmidt et al. 2017; Amor'ın et al. 2017; Naidu et al. 2022b; Bunker et al. 2023a) which hints at a possible link between C IV and high Ly α escape. This latter consideration has been indirectly proposed in a number of recent studies, which link high C IV /C III ] ratios to a nonzero escape fraction of Lyman continuum photons (Berg et al. 2019; Schaerer et al. 2022; Mascia et al. 2023b; Kramarenko et al. 2024) and thus by extension Ly α . Whether such systems are representative of the  \nFigure 12. A comparison between the median spectrum of sources with EW 0 , Ly α > 30 ˚ A (LAEs; blue) and those with EW 0 , Ly α < 10 ˚ A (non-LAEs; black). The two spectra are well matched in median redshift and absolute magnitude, allowing for a close comparison of intrinsic properties: the two spectra reveal modest differences in UV slope, but significant differences in line intensity and ratios (e.g., the presence or otherwise of C IV , C III ], and He II +[O III ] in the rest-frame UV), as well as in the rest-frame optical continuum. Overall, the LAE stack appears to represent sources with more extreme ionizing fields, lower metallicities and stellar masses, and young stellar populations, suggesting intrinsic properties - in addition to environmental effects - play some role in regulating the visibility of strong Ly α .  \n<!-- image -->  \nglobal population at z > 6 remains unclear, although our stacked results suggest they likely are, at least in Ly α -emitting systems.  \nWe also find notable difference in the strengths of the rest-frame optical line emission, which we list (fluxes and EWs) in Table 6 in the Appendix. The LAE spectrum displays line strengths × 1.4-1.7 larger than its non-LAE counterpart for [Ne III ] λλ 3969,3870 ˚ A, H δ , and H γ +[O III ], and moderately higher factors of ∼ × 2 in its He I λ 4473 ˚ A, [O III ]+H β , and H α emission. The largest difference is seen in the [O III ]+H β line strength, with EWs of 1573 ± 34 ˚ A and 734 ± 11 ˚ A for the LAE and non-LAE stacks, respectively. The two measurements are in excellent agreement with the differences observed by Castellano et al. (2017) for LAEs and nonLAEs through stacked Spitzer/IRAC photometry, as well as more general measurements inferred from both SED fitting (e.g., Roberts-Borsani et al. 2016; Endsley et al. 2021b, 2023; Laporte et al. 2023) or NIRSpec spectroscopy (e.g., Williams et al. 2022; Tang et al. 2023; Saxena et al. 2023a). The more extreme environments inferred from enhanced [O III ]+H β EWs are reflected in the R 23O 32 ratios of the two stacks, with  \nR 23 LAE = 6 . 12 ± 0 . 23 and O 32 LAE = 21 . 00 ± 0 . 57 compared to R 23 non -LAE = 6 . 54 ± 0 . 22 and O 32 non -LAE = 6 . 69 ± 0 . 30. While the R 23 ratios of the two stacks are consistent with each other, the LAE stack displays a O 32 ratio more than × 3 that of its non-LAE counterpart and ∼ × 1 . 6 -3 . 0 larger than our redshift-binned stacks, clearly indicating a higher ionization parameter likely paired with lower metallicity. The inferred metallicities for the two stacks using the O 32 prescription of Sanders et al. (2023b) are log (O/H) LAE = 7 . 48 ± 0 . 01 and log (O/H) non -LAE = 7 . 91 ± 0 . 02, supporting the interpretation (see Cullen et al. 2020).  \nFinally, we also observe significant differences in the rest-frame optical continuum. Our LAE spectrum displays a Balmer index of F red /F blue = 0 . 79 ± 0 . 07, similar to the values observed at z ∼ 8 . 5 in our redshift-binned stacks (see Figure 6) which suggests young stellar ages and enhanced contributions from the nebular continuum. Our non-LAE stack instead displays a larger Balmer index of F red /F blue = 1 . 03 ± 0 . 04, more typical of aged stars and enhanced stellar contributions to the spectrum. The interpretation is corroborated by contrasting stellar mass and mass-  \nhted ages, where gsf -derived estimates point to log M ∗ , LAE /M ⊙ ∼ 7 . 54 and ∼ 4 Myr for the LAE composite, and log M ∗ , non -LAE /M ⊙ ∼ 8 . 36 with ∼ 110 Myr for the non-LAE composite. The significant contrast in line and continuum features between the two stacks is demonstrative of intrinsic properties playing some role in regulating the visibility of Ly α , and possibly the recent star formation histories of the hosts. This is likely in combination with the presence of overdense regions and/or large ionized bubbles in sources at z ⩾ 6 (c.f., Roberts-Borsani et al. 2023a), and/or with effects from a largely neutral medium for sources at z ≲ 5 . 3 -5 . 8 (see Eilers et al. 2018; Bosman et al. 2022; Zhu et al. 2023 for discussions on the end of reionization).  \nThe SED fitting of the LAE composite also indicates an effectively dust-free system ( E ( B -V ) ≃ 0 . 07), therefore as a final inspection we can use the measured H β , H α , and rest-frame UV continuum fluxes to infer a number of fundamental quantities related to the Ly α visibility. One of the most informative is the production rate of ionizing continuum photons, ξ ion . Given ξ ion is related to H α via:  \nξ ion = N ( H 0 ) L 1500 = L (H α ) (1 -f esc ) 1 L 1500 [7 . 4 × 10 11 s -1 erg -1 s -1 Hz -1 ] (2)  \nwhere N ( H 0 ) is the production rate of ionizing continuum photons, L (H α ) is the luminosity of H α , L 1500 is the continuum luminosity at approximately λ rest ≃ 1500 ˚ A, and f esc is the escape fraction of ionizing continuum photons (assumed to be zero here). From the measured H α flux ( f (H α ) = 1 . 1 × 10 -18 erg s -1 cm -2 ), we derive log ξ ion , LAE = 25 . 49 ± 0 . 01 Hzerg -1 , which is consistent with canonical values derived from SED fitting of high-redshift sources and/or low-redshift populations (e.g., Bouwens et al. 2016; Matthee et al. 2017; Shivaei et al. 2018; Simmonds et al. 2024). Applying the same procedure to the non-LAE stack, we derive a similarly consistent value of log ξ ion , non -LAE = 25 . 30 ± 0 . 01 Hzerg -1 . We note those values increase by only 0.04 dex for modest escape fractions of 10%.  \nAdditionally, we evaluate the Ly α escape fraction, f Ly α esc of our composite LAE spectrum. From Case B recombination theory, we adopt an Ly α /H β ratio of 24.273 (assuming Ly α /H α =8.7 and H α /H β =2.790) and, multiplied by our measured H β flux, f (H β ) = 4 . 0 × 10 -19 erg s -1 cm -2 , derive an intrinsic Ly α flux of f Ly α int = 9 . 8 × 10 -18 erg s -1 cm -2 . From a simple Gaussian fit to the continuum-subtracted Ly α profile, we derive an  \nobserved Ly α flux of f Ly α obs = 3 . 5 × 10 -18 erg s -1 cm -2 which equates to an escape fraction of approximately f Ly α esc ≃ 36%, on the high end of but consistent with measurements of individual objects at z ⩾ 5 reported in the literature (see Saxena et al. 2023a; Napolitano et al. 2024).  \nA more in-depth interpretation of such values and their implications for the contributions of galaxies to the reionization process requires a thorough investigation into the underlying properties of the sample going into our stack (which in this case probes the strongest LAEs), which we leave to future papers. Nevertheless, our comparison of the average spectra of LAEs and nonLAEs very clearly illustrates important intrinsic differences between the two populations. Considering the identical median redshifts of the stacks ( z ∼ 5 . 9), those differences are likely related to the physics within the galaxies themselves and their immediate surroundings as opposed to the neutrality of the intervening IGM. These differences will have to be taken into account when using, e.g., galaxies and their Ly α emission as a probe of the neutral density of the Universe. For example, the high Ly α detection fraction in z > 7 samples selected for their extreme optical emission lines (RobertsBorsani et al. 2016) is easier to reproduce if they have intrinsically higher Ly α emission than UV-continuumselected galaxies (Roberts-Borsani et al. 2023a). In order to further characterize these differences, higher spectral resolution spectra are needed, so as to increase the sensitivity to Ly α and determine its velocity profile in comparison to the rest-frame optical Balmer lines. Our stacks, and average line ratios, are therefore intended as an initial benchmark for reference in future studies.", '6. SUMMARY AND CONCLUSIONS': "We have performed the largest compilation to date of z ⩾ 5 sources from public JWST/NIRSpec prism spectra. We identify 482 sources with secure redshifts primarily from the [O III ]+H β emission-line complex at z ≃ 5 . 0 -9 . 7 and the Lymanα break at higher redshifts, with 356 of those sources representing novel confirmations. After removal of candidate unobscured AGN sources, we verify that our sample is commensurate with the star-forming main sequence at comparable redshifts, and use our statistical sample to construct representative and high SNR composite spectra as a function of redshift, absolute magnitude, stellar mass, and Ly α EW. We report on their average properties while making those templates available to the community. Our main findings can be summarized as follows.  \n- · We provide the tightest constraints yet on the UV continuum slopes of high-redshift galaxies and  \ntheir evolution between z ≃ 5 -11, finding blue slopes of β ≃ -2 . 3 to -2 . 7, indicative of overall dust-poor ISM and young stellar populations. Sources at z ∼ 5 . 5 display the reddest slopes, followed by a steep decline to bluer slopes for sources at z ∼ 6 -7, and a flattening off in the evolution between z ∼ 7 -11. UV-luminous and higher-mass sources display systematically redder UVslopes compared to their UV-fainter and lowermass counterparts, likely reflecting more evolved ISM conditions. On average, sources at z > 10 do not display slopes of β ≲ -2 . 7 indicative of dust-free systems.  \n- · We report a steadily declining rest-frame optical continuum break with redshift, indicative of declining stellar ages and metallicities. The largest breaks ( F ν, 4200 /F ν, 3500 > 1) are attributed to our lowest-redshift ( z ∼ 5 . 5) stack, while our higherredshift composites are characterized by optical jumps of F ν, 4200 /F ν, 3500 < 1 attributable to extremely young and metal-poor stellar populations, and a possible enhancement of the nebular continuum.\n- · All of our composite spectra display unambiguous detections of C III ]1907,1909 ˚ A emission, with equivalent widths of EW 0 ≃ 5 -14 ˚ A and a clear increase of the line strength with redshift. The evolving EW of the line is consistent with an extrapolation of results at z ∼ 0 -3, while its prevalence and strength in our stacks suggests a high cadence among individual spectra. A comparison to confirmed objects in the literature suggests a large number classified as 'extreme' in fact reside close to or on the relation, reflecting C III ] strengths comparable to the general star-forming population. Prominent He II +[O III ] line emission is also consistently detected, but at comparatively low strengths.\n- · Ratios of strong rest-optical line emission ( O 32, R 23, Ne 3 O 2, and (Ne III +[O III ])/H δ ) indicate subsolar gas-phase metallicities (log (O/H) ≃ 7 . 2 -7 . 9; ∼ 0.03-0.16 Z ⊙ ) which decline with redshift. Comparisons to photoionization models further reveal that z ⩾ 5 galaxies populate regions of ionization versus excitation diagrams typically associated with especially strong ionization parameters, indicative of extreme ionizing radiation fields not observed at z ∼ 0.\n- · Lastly, a direct comparison between the average (redshift- and M UV -matched) spectra of LAEs and  \nnon-LAEs reveals important differences in their intrinsic properties. The spectrum of the former displays enhanced UV emission (from C III ], C IV , and He II +[O III ]), bluer UV continuum slopes, an inverse Balmer jump, and stronger line emission with enhanced ratios compared to the latter. Inferences from H α and Ly α emission in the LAE stack suggest a production rate of ionizing continuum photons of log ξ ion = 25 . 49 ± 0 . 01 Hzerg -1 and Ly α escape fraction of f Ly α esc ≃ 35%, consistent with canonical values and measurements of individual highz sources in the literature. The observed intrinsic differences suggest LAEs are governed by more extreme ionizing fields, young stellar populations, dust-poor ISM, and low metallicities, the combination of which - combined with environmental effects - may be more conducive to higher escape fractions of Ly α photons compared to non-LAEs regardless of environmental effects.  \nThe combination of JWST/NIRSpec's unprecedented sensitivity and wavelength coverage is opening a new window with which to characterize galaxies over the first ∼ 1 . 2 Gyr of the Universe. The plethora of available spectroscopic data sets from only ∼ 1 . 5 years of operations now allow us to move from single-object characterizations to statistical analyses generally attributed to photometric studies, placing novel and direct constraints on the evolution of galaxy properties across entire populations and out to the earliest times. This study represents one of the first (and largest) examples to do so, mirroring spectroscopic studies at z ∼ 2 -4 and extending them to star-forming sources representative of the populations dominating the aftermath, the conclusion, and the thick of the reionization process. The findings and templates derived in this study represent a new and observational benchmark for future studies of highredshift galaxies, allowing for a more direct comparison of galaxy evolution across cosmic time and acting as a spectroscopic reference point for the interpretation of individual sources out to the highest redshifts.", 'ACKNOWLEDGMENTS': "We are grateful and extend our thanks to a number of people. First, we thank Adi Zitrin and Asheesh Meena for providing the magnification factors of confirmed sources behind the MACS 0647 cluster. We also thank Anna Feltre, St'ephane Charlot, and Alex Cameron for providing and discussing the photoionization models used in this study. G.R.B. also thanks Gabe Brammer for helpful discussions on the reduction of the JWSTspectra, and for his extensive efforts in developing and maintaining the msaexp code from which the highz community has greatly benefited. Last but not least, we also thank Daniel Schaerer, Rui Marques-Chaves, and Emma Giovinazzo for additional models used in this study and for useful discussions on the science described here. Some/all of the data presented in this article were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via DOIs 10.17909/j4rg-3f57 and 10.17909/8tdj-8n28.", 'REFERENCES': "- Aihara, H., Allende Prieto, C., An, D., et al. 2011, ApJS, 193, 29, doi: 10.1088/0067-0049/193/2/29  \nAmor'ın, R., Fontana, A., P'erez-Montero, E., et al. 2017,  \nNature Astronomy, 1, 0052, doi: 10.1038/s41550-017-0052  \nArrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al.  \n2023a, arXiv e-prints, arXiv:2304.05378,  \ndoi: 10.48550/arXiv.2304.05378  \n-. 2023b, arXiv e-prints, arXiv:2303.15431,  \ndoi: 10.48550/arXiv.2303.15431  \nAsplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481,  \ndoi: 10.1146/annurev.astro.46.060407.145222  \nBagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2023, ApJL, 946, L12, doi: 10.3847/2041-8213/acbb08  \nBakx, T. J. L. C., Sommovigo, L., Carniani, S., et al. 2021, MNRAS, 508, L58, doi: 10.1093/mnrasl/slab104  \nBalogh, M. L., Morris, S. L., Yee, H. K. C., Carlberg, R. G., & Ellingson, E. 1999, ApJ, 527, 54, doi: 10.1086/308056 Barbary, K. 2016, Journal of Open Source Software, 1, 58, doi: 10.21105/joss.00058  \nBerg, D. A., Chisholm, J., Erb, D. K., et al. 2019, ApJL, 878, L3, doi: 10.3847/2041-8213/ab21dc  \nBergamini, P., Acebron, A., Grillo, C., et al. 2023, arXiv e-prints, arXiv:2303.10210,  \ndoi: 10.48550/arXiv.2303.10210  \nBertin, E., & Arnouts, S. 1996, A&AS, 117, 393,  \ndoi: 10.1051/aas:1996164  \nBezanson, R., Labbe, I., Whitaker, K. E., et al. 2022, arXiv e-prints, arXiv:2212.04026,  \ndoi: 10.48550/arXiv.2212.04026  \nBosman, S. E. I., Davies, F. B., Becker, G. D., et al. 2022, MNRAS, 514, 55, doi: 10.1093/mnras/stac1046  \nBouwens, R. J., Smit, R., Labb'e, I., et al. 2016, ApJ, 831, 176, doi: 10.3847/0004-637X/831/2/176  \nBouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2014, ApJ, 793, 115, doi: 10.1088/0004-637X/793/2/115 -. 2015, ApJ, 803, 34, doi: 10.1088/0004-637X/803/1/34  \nBouwens, R. J., Smit, R., Schouws, S., et al. 2021, arXiv e-prints, arXiv:2106.13719.  \nhttps://arxiv.org/abs/2106.13719  \nBrammer, G. 2021, eazy-py, 0.5.2,  \ndoi: 10.5281/zenodo.5012704  \n-. 2023, grizli, doi: 10.5281/zenodo.1146904  \nBruzual A., G. 1983, ApJ, 273, 105, doi: 10.1086/161352  \nBunker, A. J., Saxena, A., Cameron, A. J., et al. 2023a, arXiv e-prints, arXiv:2302.07256,  \ndoi: 10.48550/arXiv.2302.07256  \nBunker, A. J., Cameron, A. J., Curtis-Lake, E., et al. 2023b, arXiv e-prints, arXiv:2306.02467,  \ndoi: 10.48550/arXiv.2306.02467  \nCalabro, A., Castellano, M., Zavala, J. A., et al. 2024, arXiv e-prints, arXiv:2403.12683, doi: 10.48550/arXiv.2403.12683", 'A. REDSHIFT CONFIRMATIONS OF z ⩾ 5 SOURCES': 'Here we tabulate the spectroscopic redshifts and a number of spectrophotometric properties associated with our parent sample of confirmed z ⩾ 5 sources. In Table 3 we present all sources (star-forming sources and those flagged as potential AGN, see Section 3) with secure redshifts and where sufficient photometric and spectroscopic coverage existed for a scaling of the spectrum in the rest-frame UV. In Table 4, we adopt the same format as Table 3 but list sources where either the redshift was uncertain (e.g., where the msaexp -derived redshift was z ⩾ 5 but the spectrum lacked clearly identifiable features) and/or where UV coverage necessary for a scaling of the spectrum was lacking.  \nTable 3 . The List of Spectroscopically Confirmed Sources Reduced and Compiled in This Study, along with Several Spectrophotometric Properties.  \nNote -Star-forming sources are listed in the top half of the table, while AGN candidates (identified here and in other works) are listed in the bottom half. Only the first 10 sources from each classification are shown here, while the table is available in its entirety online. All apparent magnitudes remain uncorrected for gravitational lensing, however absolute magnitudes are corrected for such effects. Objects for which coordinates are absent represent multiply-imaged sources where the apparent and absolute magnitudes are derived by de-magnifying the photometry of the individual images and taking the median value.  \nTable 4. The Same as Table 3 But for Sources with Uncertain Redshifts and/or Where Rest-frame Ultraviolet Photometry or Prism fluxes Were Not Available for Normalization Purposes.', 'B. AVERAGE SLIT SLOSS CORRECTION FUNCTION': 'The correction of wavelength-dependent slit losses (when present) in galaxy spectra represents a crucial step toward obtaining accurate measurements. As mentioned in Section 3.2, residual wavelength-dependent slit losses are found to be present in a number of our reduced spectra that are not fully accounted for in the path loss step of the data reduction pipeline. Here we expand on the determination of those slit losses and the derivation of an average correction function.  \nFor each photometric band redward of the adopted UV-scaling filter, we evaluate the percentage difference between the extracted photometry and prism pseudophotometry generated directly from the spectrum with the relevant filter response curves (i.e., NIRCam photometry minus prism photometry, divided by NIRCam photometry times 100). The results are plotted in Figure 13, where for each source and photometric band the percentage of slit loss is plotted as a gray point. The distribution of points for each filter are approximately Gaussian, and thus to highlight their centroids and spread, we fit each density profile with a Gaussian function and plot the best fit as colored lines (out to their 5 σ width limits), along with a marker (shifted in wavelength for clarity) indicating the centroid and its 1 σ uncertainty. We note a nonnegligible amount of wavelength-dependent slit loss, starting at approximately 1-2% (absolute values) on average in the blue ( < 2 µ m) bands (albeit with some scatter in the F115W and F150W bands) and increasing to 10-20% on average in the red bands, undoubtedly reflecting loss of light from an increasing point-spread function width with wavelength. We note that the apparently larger deviations in the F115W and F150W centroids are due to poorly constrained density profiles, which begin to deviate from a well-defined Gaussian. The average wavelength-dependent slit loss is well fit by an exponential function (plotted as a turquoise line), which begins at 10-20% slit loss in the reddest bands and progressively decreases to subpercent levels by ∼ 3 µ m and shorter-wavelength bands. We use this parameterization to correct all composite spectra derived in this analysis (although we do not propagate uncertainties on the slit loss correction into our spectral measurements) and render the correction function publicly available for community use.  \nWe also note there is nonnegligible scatter among individual sources, some of which is not uniquely attributable to loss of light due to the NIRSpes MSA shutter. These include, in a few cases, some possible oversubtraction of the sky background or overcorrection of the path loss by the data reduction pipeline, in which case the spectrum would lie below and above the photometric data points, respectively, and result in values beyond ± 100%. Accurate corrections on an object-by-object basis are not trivial, requiring knowledge of the object size and position in the NIRSpec MSA shutter as well as a robust photometric baseline, which is not always the case in low SNR regimes where even HST and NIRCam photometry can display significant scatter. However, an average correction for stacked spectra is both possible and required: the high SNR of composite spectra ensures a more robust parameterization of the trend which is less sensitive to outliers (e.g., apparently faint galaxies or low SNR spectra) and encapsulating of various degrees of slit loss deriving from, e.g., size and morphology and/or position in the slit.  \nFigure 13. The percentage slit loss of sources observed with NIRSpec prism compared to NIRCam photometry. The slit loss is defined as the NIRCam photometry minus the NIRSpec pseudophotometry, divided by the NIRCam photometry, and multiplied by 100. Measurements for individual sources are plotted as grey points, while Gaussian fits to the density distribution in each filter are plotted as solid lines with normalized amplitudes. The centroid of each Gaussian is shown as a colored square, shifted by +0.2 µ mfor clarity. A best-fit exponential function to the squares is shown as a turquoise line. We note that the distributions of single-object points for virtually all bands are Gaussian in shape, except the F115W and F150W bands, which are less well defined and thus appear slightly up- or down-scattered compared to their redder counterparts.  \n<!-- image -->', 'C. LINE FLUXES AND EQUIVALENT WIDTHS': 'We list here the H β -normalized integrated line fluxes and EW 0 s, calculated as described in Section 3.7 and discussed throughout the various sections of the paper.  \nTable 5. Emission line Ratios Relative to H β and Rest-frame Equivalent Widths for the z ∼ 5 . 5 -9 . 5 Redshift-binned and Ly α -Emitter Composite Spectra listed in Table 2.  \nNote -The H β fluxes are 3 . 3( ± 0 . 1) × 10 -19 , 3 . 1( ± 0 . 1) × 10 -19 , 2 . 9( ± 0 . 1) × 10 -19 , and 2 . 7( ± 0 . 3) × 10 -19 , and 2 . 8( ± 0 . 3) × 10 -19 erg s -1 cm -2 for the z5 SF, z6 SF, z7 SF, z8 SF, and z9 SF stacks, respectively. All equivalent widths are quoted in the rest-frame and have units of ˚ A.  \nTable 6. The Same as Table 5 But for the Ly α Emitter and non-Ly α Emitter Stacks from Table 2.  \nNote -The H β fluxes are 4 . 0( ± 0 . 1) × 10 -19 erg s -1 cm -2 for the LAE stack and 3 . 0( ± 6 . 0) × 10 -19 erg s -1 cm -2 for the non-LAE stack. All equivalent widths are quoted in the rest-frame and have units of ˚ A.'}