One of the surprising early findings with JWST has been the discovery of a strong rollover or a softening of the absorption edge of Lyalpha in a large number of galaxies at zgtrsim 6 in addition to systematic offsets from photometric redshift estimates and fundamental galaxy scaling relations. This has been interpreted as damped Lyalpha absorption DLA wings from high column densities of neutral atomic hydrogen HI signifying major gas accretion events in the formation of these galaxies. To explore this new phenomenon systematically we assemble the JWSTNIRSpec PRImordial gas Mass AssembLy PRIMAL legacy survey of 494 galaxies at z5.513.4. We characterize this benchmark sample in full and spectroscopically derive the galaxy redshifts metallicities starformation rates and ultraviolet slopes. We define a new diagnostic the Lyalpha damping parameter Drm Lyalpha to measure and quantify the Lyalpha emission strength HI fraction in the IGM or local HI column density for each source. The JWSTPRIMAL survey is based on the spectroscopic DAWN JWST Archive DJASpec. All the software reduced spectra and spectroscopically derived quantities and catalogs are made publicly available in dedicated repositories. The fraction of strong galaxy DLAs are found to be in the range 6595 at zgt5.5. The fraction of strong Lyalpha emitters LAEs is found to increase with decreasing redshift in qualitative agreement with previous observational results and are predominantly associated with lowmetallicity and UV faint galaxies. By contrast strong DLAs are observed in galaxies with a variety of intrinsic physical properties. Our results indicate that strong DLAs likely reflect a particular early assembly phase of reionizationera galaxies at which point they are largely dominated by pristine HI gas accretion. abridged
['Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
The JWSTPRIMAL Legacy Survey. A JWSTNIRSpec reference sample for the physical properties and Lymanalpha absorption and emission of sim 500 galaxies at z5.513.4
{'A JWST/NIRSpec reference sample for the physical properties and Lymanα absorption and emission of ∼ 500 galaxies at z = 5 . 5 -13 . 4': "K. E. Heintz, 1 , 2 , 3 G. B. Brammer, 1 , 2 D. Watson, 1 , 2 P. A. Oesch, 3 , 1 , 2 L. C. Keating, 4 M. J. Hayes, 5 Abdurro'uf, 6 , 7 K. Z. Arellano-Córdova, 4 A. C. Carnall, 4 C. R. Christiansen, 1 , 2 F. Cullen, 4 R. Davé, 4 P. Dayal, 8 A. Ferrara, 9 K. Finlator, 10 , 1 , 2 J. P. U. Fynbo, 1 , 2 S. R. Flury, 11 V. Gelli, 1 , 2 S. Gillman, 1 , 12 R. Gottumukkala, 1 , 2 K. Gould, 1 , 2 T. R. Greve, 1 , 12 S. E. Hardin, 13 T. Y.-Y Hsiao, 6 , 7 A. Hutter, 1 , 2 P. Jakobsson, 14 M. Killi, 15 N. Khosravaninezhad, 16 , 1 , 2 P. Laursen, 1 , 2 M. M. Lee, 1 , 12 G. E. Magdis, 1 , 12 , 2 J. Matthee, 17 R. P. Naidu, 18 ⋆ D. Narayanan, 19 , 1 , 2 C. Pollock, 4 M. Prescott, 10 V. Rusakov, 1 , 2 M. Shuntov, 1 , 2 A. Sneppen, 1 , 2 R. Smit, 20 21 C. Terp, 1 , 2 S. Toft, 1 , 2 F. Valentino, 22 , 1 A. P. Vijayan, 1 , 12 J. R. Weaver, 11 J. H. Wise, 13 J. Witstok, 23 , 24  \nN. R. Tanvir,  \n(A ffi liations can be found after the references)  \nSubmitted -, accepted -, published -", 'ABSTRACT': "Context. One of the surprising early findings with JWST has been the discovery of a strong 'roll-over' or a softening of the absorption edge of Ly α in a large number of galaxies at z ≳ 6, in addition to systematic o ff sets from photometric redshift estimates and fundamental galaxy scaling relations. This has been interpreted as strong cumulative damped Ly α absorption (DLA) wings from high column densities of neutral atomic hydrogen (H i ), signifying major gas accretion events in the formation of these galaxies.  \nAims. To explore this new phenomenon systematically, we assemble the JWST / NIRSpec PRImordial gas Mass AssembLy (PRIMAL) legacy survey of 494 galaxies at z = 5 . 5 -13 . 4, designed to study the physical properties and gas in and around galaxies during the reionization epoch.  \nMethods. Wecharacterize this benchmark sample in full and spectroscopically derive the galaxy redshifts, metallicities, star-formation rates, and ultraviolet slopes. We define a new diagnostic, the Ly α damping parameter D Ly α to measure and quantify the net e ff ect of Ly α emission strength, H i fraction in the IGM, or local H i column density for each source. The JWST -PRIMAL survey is based on the spectroscopic DAWN JWST Archive (DJA-Spec). We describe DJA-Spec in this paper, detailing the reduction methods, the postprocessing steps, and basic analysis tools. All the software, reduced spectra, and spectroscopically derived quantities and catalogs are made publicly available in dedicated repositories.  \nResults. We find that the fraction of galaxies showing strong integrated DLAs with N HI > 10 21 cm -2 only increases slightly from ≈ 60% at z ≈ 6 up to ≈ 65 -90% at z > 8. Similarly, the prevalence and prominence of Ly α emission is found to increase with decreasing redshift, in qualitative agreement with previous observational results. Strong Ly α emitters (LAEs) are predominantly found to be associated with low-metallicity and UV faint galaxies. By contrast, strong DLAs are observed in galaxies with a variety of intrinsic physical properties, but predominantly at high redshifts and low metallicities.  \nConclusions. Our results indicate that strong DLAs likely reflect a particular early assembly phase of reionization-era galaxies, at which point they are largely dominated by pristine H i gas accretion. At z = 8 -10, this gas gradually cools and forms into stars that ionize their local surroundings, forming large ionized bubbles and produce strong observed Ly α emission at z < 8.  \nKey words. Galaxies: formation, evolution, high-redshift - Reionization - Intergalactic medium - Interstellar medium - Spectroscopy", '1. Introduction': "The first epoch of galaxy formation is believed to have occurred at z ∼ 15 -20 (Hashimoto et al. 2018; Robertson 2022), some 100 -200 million years after the Big Bang. This process was primarily driven by the accretion and gravitational collapse of primordial, neutral gas onto dark matter halos (White & Rees 1978; Blumenthal et al. 1984; White & Frenk 1991; Kereš et al. 2005; Schaye et al. 2010; Dayal & Ferrara 2018), which eventually led to the formation of stars in galaxies. The strong UV radiation originating from the first stars and supermassive black holes gradually ionized their immediate and then large-scale sur-  \nndings, initiating the reionization epoch, currently estimated to have ended by z ≈ 5 -6 (e.g. Stark 2016; Dayal & Ferrara 2018; Keating et al. 2020; Robertson 2022; Bosman et al. 2022; Fan et al. 2023). When these first stars ended their lives as supernovae, the pristine gas was chemically enriched with elements formed through explosive nucleosynthesis (Hoyle 1954; Cameron 1957; Woosley & Weaver 1995). Charting these independent components, their interplay through the 'baryon cycle' (Tumlinson et al. 2017; Péroux & Howk 2020), and their evolution with cosmic time is at the heart of contemporary galaxy formation and evolution studies.  \nThe launch of the James Webb Space Telescope ( JWST ) has now enabled us to peer deep into this early cosmic epoch, prob-  \ning the rest-frame optical and ultraviolet (UV) emission of galaxies, potentially all the way back to redshifts z ≳ 15 (Castellano et al. 2022; Naidu et al. 2022; Harikane et al. 2023; Atek et al. 2023a; Donnan et al. 2023; Bouwens et al. 2023; Wang et al. 2023b; Austin et al. 2023). This leap in near-infrared (NIR) capability now enables detailed characterization of the physical properties and baryonic components of galaxies, reaching the first epoch of galaxy formation, the cosmic dawn. This epoch seems to contain a larger population of more luminous ( M UV < -20) and more massive galaxies than expected from theoretical models (Labbé et al. 2023; Finkelstein et al. 2023a; Franco et al. 2023; Harikane et al. 2023; Adams et al. 2023; Casey et al. 2023; Bouwens et al. 2023; Chemerynska et al. 2023; McLeod et al. 2024). Various possibilities have been considered to resolve this discrepancy, including a varying initial mass function (IMF; Haslbauer et al. 2022; Trinca et al. 2023, though see Rasmussen Cueto et al. 2023), radiation pressure pushing dust away from star-forming regions (Ferrara et al. 2023), bursty star-formation histories (Sun et al. 2023), or an unexpected overabundance of active galactic nuclei (AGN; Inayoshi et al. 2022; Pacucci et al. 2022; Dayal et al. 2024). The apparent discrepancy may be exaggerated by selection biases towards the youngest, most highly star-forming galaxies with more bursty star formation histories at z > 10 (Mason et al. 2023; Mirocha & Furlanetto 2023; Shen et al. 2023), making the galaxies appear temporarily brighter in the UV. As neutral gas accretion is the primary driver and fuel for star formation, a robust measure of this component is key to resolving the debate.  \nRecent JWST spectroscopic studies have pushed the limit of the most distant galaxies confirmed to z ≈ 11-13 (Curtis-Lake et al. 2023; Wang et al. 2023a; Fujimoto et al. 2023b; Zavala et al. 2024; Castellano et al. 2024). The metal content of galaxies during the reionization epoch can now also be readily measured based on direct Te -based methods for su ffi ciently bright galaxies (e.g., Schaerer et al. 2022; Arellano-Córdova et al. 2022; Taylor et al. 2022; Brinchmann 2023; Curti et al. 2023a; Katz et al. 2023; Rhoads et al. 2023; Trump et al. 2023; Heintz et al. 2023c; Nakajima et al. 2023; Sanders et al. 2024) or inferred through rest-frame optical strong-line calibrations (Heintz et al. 2023a; Langeroodi et al. 2023; Curti et al. 2023b; Matthee et al. 2023). These observations show a substantial o ff set from the fundamental-metallicity relation (FMR, Mannucci et al. 2010; Lara-López et al. 2010; Curti et al. 2020; Sanders et al. 2021) established at z ≈ 0 -3 by up to 0.5 dex at z > 7 (Heintz et al. 2023a), suggesting copious pristine H i gas inflows onto galaxies at this epoch before chemical enrichment from star formation can catch up to equilibrium (e.g., Torrey et al. 2018). The exact redshift for this transition is still debated (Curti et al. 2023b; Nakajima et al. 2023).  \nConstraining the inflow of neutral, pristine H i gas is key to understanding the assembly of the primordial matter and the build-up of stars and metals in the first galaxies. However, the circumstantial evidence from the FMR for excessive H i gas accretion needs to be corroborated by measuring directly the neutral, atomic hydrogen (H i ) - the key missing element in the baryonic matter budget of high-redshift galaxies. Due to the weakness of the hyperfine H i 21 cm line, H i has historically been measured at redshifts z ≳ 2 through Lymanα (Ly α ) absorption in bright background sources such as quasars (Wolfe et al. 1986, 2005; Prochaska & Wolfe 2009; Noterdaeme et al. 2012; Péroux & Howk 2020) or gamma-ray bursts (Jakobsson et al. 2006; Prochaska et al. 2007; Fynbo et al. 2009; Tanvir et al. 2019; Heintz et al. 2023b). However, these only probe H i in narrow, pencil-beam sightlines and often at large impact parameters. Ly α  \nabsorption has also been observed directly in integrated spectra of Lyman-break galaxies at z ≈ 3 (Pettini et al. 2000; Shapley et al. 2003; Steidel et al. 2010; Cooke & O'Meara 2015; Lin et al. 2023), but only in rare, extreme cases.  \nRecently, Heintz et al. (2023d) reported the discovery of extremely strong ( N HI ≳ 10 22 cm -2 ) damped Ly α absorption systems (DLAs) in star-forming galaxies at z > 8 (see also Umeda et al. 2023; D'Eugenio et al. 2023; Chen et al. 2023). These seem to represent a galaxy population showing more prominent and substantially more prevalent DLAs than observed in Lymanbreak galaxies at z ≈ 3 (e.g. Shapley et al. 2003). The shape of the Ly α damping wings of galaxies had until these discoveries mainly been thought to trace the bulk H i of the intergalactic medium (IGM) at z ≳ 6 (Miralda-Escudé 1998; McQuinn et al. 2008), and had been used to infer the reionization history of the Universe with JWST spectroscopy (Chen 2023; Keating et al. 2023a,b; Umeda et al. 2023). The discovery of these massive interstellar or circumgalactic H i gas reservoirs carry important implications for our study of the early Universe: they constitute the first direct evidence for the mass accretion of primordial intergalactic gas into ordinary galaxies, i.e. these DLAs represent the creation of the baryonic component of galaxies. They also strongly a ff ect observations of early galaxies because they change the shape of the Ly α damping wings, altering and complicating inferences on the state of the bulk IGM. Further, they hinder the escape of ionizing photons (Laursen et al. 2011; Verhamme et al. 2015; Steidel et al. 2018; Hu et al. 2023; Hayes & Scarlata 2023), and may also be responsible for the bias in photometric redshifts compared to spectroscopic measurements (Heintz et al. 2023d; Finkelstein et al. 2023b; Fujimoto et al. 2023a; Hainline et al. 2023) of galaxies during the reionization era. Their discovery, however, also presents a fortuitous new avenue to directly probe the build-up of pristine H i in galaxies at this critical epoch.  \nIn this work, we systematically characterize these e ff ects by measuring the prevalence and prominence of strong Ly α emission and damped Ly α absorption in a large spectroscopic sample of star-forming galaxies at z = 5 . 5- 13.4 observed with JWST / NIRSpec. The observational data for these galaxies are all presented here as part of the DAWN JWST Archive (DJA), which includes fully reduced and post-processed 1D and 2D spectra of all public JWST / NIRSpec data. The large JWST / NIRSpec archive forms the basis of the JWST / NIRSpec PRImordial gas Mass AssembLy (PRIMAL) legacy survey presented here. In this first part of the survey we aim to disentangle the Ly α damping wings produced from H i gas in the immediate surroundings of galaxies (ISM or CGM) from the e ff ects of an increasingly neutral IGM, and establish statistical correlations for these features and the presence of Ly α emission with the physical properties of the galaxies. We further make the reduced spectra, source catalogs, spectroscopic redshifts, line identifications and fluxes, and the physical properties of each of the galaxies in the JWST -PRIMAL sample publicly available on dedicated webpages.  \nWe have structured the paper as follows. In Sect. 2 we present the observations, describe the spectroscopic reduction and post-processing, and outline the JWST -PRIMAL sample compilation. In Sect. 3, we detail the spectroscopic analysis of Ly α , the nebular emission line fluxes and equivalent widths, and the inferred physical properties of the sample galaxies. In Sect. 4, we consider the cosmic evolution of galaxy DLAs and Ly α emitters and quantify the physical correlations driving or hampering the observed Ly α emission and excess damped Ly α absorption. In Sect. 5, we discuss and conclude our work and provide a future outlook. All the spectroscopic data products and the analy-  \nTable 1 Archival references for the JWST -PRIMAL sample.  \nNotes. Column (1): Right ascension in degrees (J 2000). Column (2): Declination in degrees (J 2000). Column (3): Program ID under which the object was observed. Column (4): Designated MSA source ID. Column (5): Rest-frame UV signal-to-noise ratio. Column (6): Spectroscopic redshift. Column (7): The original survey paper references. A full version of this table can be found online.  \nReferences. (B23) Bunker et al. (2023a); (AH23) Arrabal Haro et al. (2023); (ECL23) Curtis-Lake et al. (2023); (W23) Wang et al. (2023a).  \nsis and results presented in this work are made publicly available on dedicated webpages. Throughout the paper, we assume concordance flat Λ CDM cosmology, with H 0 = 67 . 4 km s -1 Mpc -1 , Ω m = 0 . 315, and ΩΛ = 0 . 685 (Planck Collaboration et al. 2020).", '2.1. The DAWN JWST Archive (DJA) - Spectroscopy': 'The data considered in this work are processed as part of the DAWN JWST Archive (DJA), an online repository containing reduced images, photometric catalogs, and spectroscopic data for public JWST data products 1 . Here, we detail the spectroscopic data reduction and post-processing only, but see Valentino et al. (2023) for details on the imaging and photometric data. The DJA spectroscopic archive (DJA-Spec) is, at the time of writing, comprised of observations taken from some of the large Early Release Science (ERS), General Observer (GO) and Guaranteed Time (GTO) Cycle 1 & 2 programs, such as CEERS (ERS-1345; Finkelstein et al. 2023a), GLASS-DDT (DDT-2756 Treu et al. 2022; Roberts-Borsani et al. 2023) (ERS-1324; Treu et al. 2022), JADES (GTO-1180, 1210, GO-3215; Bunker et al. 2023a; Eisenstein et al. 2023), and UNCOVER (GO-2561; Bezanson et al. 2022). Many of the spectra analyzed in this work have been presented in dedicated survey papers cited above and in others (e.g., Arrabal Haro et al. 2023; Atek et al. 2023b; Cameron et al. 2023b), which also describe the sample selection and observational design of the di ff erent programs. With DJA-Spec we reduce and extract all of the spectra from this diverse array of observational programs with a single standardized pipeline.', '2.2. DJA spectroscopic reduction and data processing': "All data processing is performed with the publicly available G ri -zli (Brammer 2023a) 2 and MSAE xp (Brammer 2023b) 3 soft-  \nware modules, with the procedure as follows. The spectroscopic analysis begins with reprocessing the individual uncalibrated ( uncal ) exposures retrieved from the Mikulski Archive for Space Telescopes (MAST) with the JWST Detector1Pipeline. 4 . The pipeline is executed with the default parameters, but with snowblind 5 run between the jump and ramp\\_fit steps for improved masking of the bright cosmic ray 'snowballs'(Rigby et al. 2023). We remove a column average from the count-rate ( rate ) files produced by the first step to remove the 1 / f noise. We also find that scaling the read noise extension of the calibrated count-rate images by a factor of order ≈ 1 . 4 derived separately for each exposure is necessary to explain the variance of pixels in un-illuminated portions of the NIRSpec detectors. The corrected rate files are then run through the calwebb\\_spec2 pipeline with its default parameters up through the photometric calibration ( photom ) step. The products at this stage are fluxand wavelength-calibrated 2D spectra for every source indicated with an open shutter in the attached MSA metadata, in a frame cut out from the original detector pixel grid with curved spectral traces. These 2D SlitModel cutouts are saved to individual multiextension FITS files, and it is here that the subsequent processing with MSAE xp departs significantly from the standard jwst pipeline.  \nThe sky background of each source spectrum is removed by taking di ff erences of the 2D cutout spectra taken at di ff erent 'nod' positions of the telescope. Typically there are three nod o ff sets that shift the source by 0 '' . 5 ≈ 5 pixels within the MSA slitlet (the so-called 3-Shutter-Slitlet pattern), though the pipeline calculates the groupings automatically to handle other configurations of the heterogeneous archival observations. Average 2D source and sky spectra (and associated variance images) are calculated from exposures at a given nod position and from those at all other o ff set positions, respectively. We fit a 2D Gaussian profile to the source -sky spectra where the width is the quadrature sum of a free parameter and the wavelengthdependent width of the point spread function and where a centroid shift is fit starting from an initial guess provided by the exposure MSA metadata and the catalog used to design the MSA  \nFig. 1 The absolute rest-frame UV magnitude ( M UV) as a function of spectroscopic redshift for the JWST / NIRSpec Prism sources at z > 2 from DJA with robust spectroscopic redshifts (grade 3) shown in grey. The JWST -PRIMAL sources are highlighted by the red squares, selected by requiring z spec > 5 . 5, an integrated S / N > 3 around the redshifted Ly α wavelength region, and full NIRSpec Prism spectroscopic wavelength coverage from 0.6-5.3 µ mas outlined in Sect. 2.3.  \n<!-- image -->  \nmask. We note that the latter centroid shift is usually of order 0.1 pix (10 mas), which, while consistent with the expected precision of the input catalog astrometry and telescope pointing, is non-trivial and measurable and required for robust extraction of the 1D spectra. The sky-subtracted 2D spectra cutout spectra at di ff erent o ff set positions are combined in a single rectified frame with orthogonal wavelength and cross-dispersion pixel axes by binning with a 2D histogram algorithm. Finally, the one-dimensional spectra and uncertainties are extracted with 'optimal' weighting (Horne 1986) using the fitted 2D Gaussian profile. We note that Roberts-Borsani et al. (2024) have recently used the same pipeline presented here and developed for this work, presenting a similar e ff ort to spectroscopically characterize galaxies at z > 5 observed with NIRSpec Prism.  \nCurrently (Apr 2, 2024) DJA-Spec includes 7,319 individual combined spectra taken with the NIRSpec Prism / CLEAR and 1,665 combined spectra with NIRSpec medium and highresolution gratings. A full overview can be found on the dedicated webpage 6 , where the following data are available: best-fit and inspected redshifts, the median signal-to-noise ratios (S / N) of the spectra, grades from visual inspection on the quality of the spectra, Hubble Space Telescope ( HST ) and JWST / NIRCam imaging thumbnails, and optimally reduced and extracted 2D and 1D spectra.  \nIn this work, we focus on the JWST / NIRSpec Prism spectroscopic data from DJA due to the higher continuum sensitivity  \nwhich is essential for quantifying the evidence for Ly α emission or absorption damping wings. JWST 's NIRSpec Prism configuration covers the entire near-infrared passband from 0.6 µ m to 5.3 µ m, with a varying spectral resolution from R ≈ 50 in the blue end to R ≈ 400 in the red end (Jakobsen et al. 2022). We have carefully reduced and processed the spectra in DJA in a homogeneous way, limiting any biases between reduction and calibration codes across surveys. Further, this optimizes the consistency and uniformity of the JWST -PRIMAL sample measurements and the DJA-Spec sample as a whole.", '2.3. JWST-PRIMAL sample selection': 'The main objective of this paper is to accurately model the Ly α damping wings and / or measure strong Ly α emission, and spectroscopically derive the physical properties of galaxies during the reionization epoch. For this, we require the following criteria for the sources to enter our archival JWST -PRIMAL sample:  \n- 1. A robust, spectroscopic redshift measurement at z > 5 . 5 from a minimum of one emission line detected at > 3 σ (in addition to the Ly α break).\n- 2. An integrated signal-to-noise over the redshifted Ly α and rest-frame UV (1550 Å) regions of S / N > 3.\n- 3. Flux density measurements over the entire NIRSpec Prism wavelength coverage from 0 . 6 -5 . 3 µ m.  \nThese criteria ensure that the redshift of the foreground H i gas absorber, from interstellar to intergalactic scales, and the shape of the Ly α line profile can be accurately measured. Further, the  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFig. 2 Schematic highlighting the intrinsic physical models a ff ecting the shape of the Ly α transmission from an example galaxy at z = 9 . 0, all convolved by the nominal JWST / NIRSpec Prism spectral resolution. The default model with x HI = 0 . 5 is shown at the top left, where the grey-shaded region represents the e ff ect of a partly ( x HI = 0 . 01) to fully neutral ( x HI = 1 . 0) IGM. For illustrative purposes the expected IGM density field at z = 9 is shown in the top right, extracted from the A straeus simulations. In the bottom panels are shown the combined e ff ects of the default IGM model and various sizes of the ionized UV bubble (left, blue: R b = 1 , 10 , 50 cMpc) and DLAs from local H i gas reservoirs (right, red: N HI = 10 21 , 10 22 , 10 23 cm -2 ). The integration region for D Ly α is marked by the top line in all the Ly α transmission curve figures and is defined to encapsulate all the predicted physical scenarios.  \n<!-- image -->  \nrequirement of the full spectral coverage enables a full characterization of each source in the JWST -PRIMAL sample, including robust spectroscopic modelling of the stellar continuum and rest-frame UV spectral slope, emission line fluxes and equivalent widths EWs, and direct constraints on the star-formation rates (SFRs), metallicities, and ionization parameters for all sources.  \nIn total, 494 sources from DJA-Spec meet these criteria. Fig. 1 show the absolute UV magnitude, M UV, as a function of redshift for the JWST -PRIMAL sources and compared to the underlying DJA-Spec sample (at z > 2). Table 1 summarizes the target list, including their coordinates, original program and source ID, and emission-line redshifts.', '3. Analysis and results': 'Here we detail the spectroscopic measurements derived for the full JWST -PRIMAL sample. We focus on the Ly α damping wings, but also detail the measurements and basic physical properties of the sample galaxies. All the spectroscopically-derived quantities are made available on a dedicated webpage 7 .', '3.1. The Lymanα damping parameter': "During the reionization epoch at z ≳ 6, several factors add to the line shape of Ly α on integrated galaxy spectra: They can show strong Ly α emission (LAE; Matthee et al. 2018; Mason et al. 2018; Witstok et al. 2023), damping wings from an increasingly neutral IGM (Miralda-Escudé 1998; Keating et al. 2023a; Chen 2023), excess continuum flux in the wings related to the size of the immediate ionized bubbles (McQuinn et al. 2008; Castellano  \net al. 2016, 2018; Fujimoto et al. 2023b; Umeda et al. 2023; Hayes & Scarlata 2023), strong damped Ly α absorption (DLA) from local H i gas (Heintz et al. 2023d; D'Eugenio et al. 2023), intrinsic variations due to a changing velocity and density distribution of gas and dust in the ISM and CGM (Dayal et al. 2011; Verhamme et al. 2015, 2017; Gronke et al. 2017; Hutter et al. 2023) or even possibly be dominated by two-photon emission processes (Steidel et al. 2016; Chisholm et al. 2019; Cameron et al. 2023a). Many of these e ff ects are degenerate, so we have to disentangle them statistically. For this, we define a new simple diagnostic, which we denote the Ly α damping parameter :  \nD Ly α ≡ Z λ Ly α, up λ Ly α, low (1 -F λ/ F cont) d λ / (1 + z spec) , (1)  \nwhich is the equivalent width (EW) of the transmitted flux density, F λ over the wavelength region covered by the instrumentally-broadened Ly α transition. The continuum flux, F cont, over the same region is estimated by extrapolating the restframe UV slope, β UV, which is derived directly in each spectra from rest-frame 1400 -2600 Å (see Sect. 3.4 for further details). This approximation of the stellar continuum is generally consistent with models of the continuum flux using galaxy templates (Heintz et al. 2023d). We define the integration limits from λ Ly α, low = 1180 Å to λ Ly α, up = 1350 Å (rest-frame) to capture most of the damping feature for the most extreme cases and to simultaneously limit the contamination from absorption or emission lines at longer wavelengths. We tested various lower limits and found this to be the bluest wavelength at which information is not lost at the resolution of NIRSpec Prism.  \nA set of physical models, shown as imprints on the Ly α transmission curves, are visualized in Fig. 2 including: varying ionized bubble sizes from R b = 1 -50 cMpc, neutral hydrogen  \nfractions of the IGM from x HI = 0 . 1 -1 . 0, and DLA H i column densities N HI = 10 21 -10 22 cm -2 , all a ff ecting the shape of the Ly α damping wings. The IGM models are based on the prescription by Miralda-Escudé (1998), using the approximation described in Totani et al. (2006), and the UV bubble sizes are included by varying the upper bound of the IGM integration region, such that z IGM , upper < z gal (e.g., McQuinn et al. 2008). The DLAs with varying H i column densities are derived from the Voigt-profile approximation by Tepper-García (2006). The corresponding Ly α damping parameter D Ly α is listed for each model. The UV bubble and DLA models all assume a benchmark IGM fraction of x H i = 0 . 5, and may therefore vary with the same scatter in D Ly α as observed for the IGM-only models.  \nWe derive damping parameters for the full sample in the range D Ly α = -400 to 150 Å. With this, we define four parameter spaces probing di ff erent physical regimes: (i) cases with excess continuum flux suggesting contributions from ionized bubbles or weak LAEs have D Ly α = 0-35 Å; (ii) the e ff ect on the Ly α damping wings of an increasing neutral IGM results in D Ly α = 3550 Å; (iii) excess damped Ly α absorption from local H i (with N HI > 10 21 cm -2 ) have D Ly α > 50 Å; and (iv) strong LAEs are typically observed with D Ly α < 0 Å. Fig. 3 shows three example spectra from the JWST -PRIMAL sample. One shows an example of a broad Ly α damping wing consistent with local H i absorption with D Ly α = 61 . 9 ± 3 . 6 Å, consistent with N HI ≈ 10 22 cm -2 . Another shows strong Ly α emission with a measured D Ly α = -60 . 1 ± 8 . 0 Å.  \nThis simple diagnostic now allows us to quantify the number of sources in each distinct population. We caution, however, that the low spectral resolution of the JWST / NIRSpec Prism configuration may conceal weak Ly α emission in the continuum, which requires higher resolution spectroscopy to resolve (see e.g. the case for GN-z11; Bunker et al. 2023b). This will likely cause a slight decrease in D Ly α . Further, as demonstrated in the bottom panel of Fig. 3, rare cases ( ≲ 1% of the sample) that show both a substantial damped Ly α broadening and simultaneously a strong Ly α emission line will not be accurately identified, here with D Ly α = 19 . 3 ± 0 . 7 Å, misleadingly suggesting mild excess Ly α continuum flux (the continuum flux of this particular source might be dominated by two-photon nebular emission, e.g., Cameron et al. 2023a). The JWST / NIRSpec shutter may also not cover the regions of the strongest Ly α emission, thereby underestimating the total Ly α photon output of the target source (Jiang et al. 2023). The D Ly α parameter is thus most powerful in identifying the most extreme cases, such as strong Ly α emitters or galaxy DLAs with N HI > 10 21 cm -2 , which uniquely have D Ly α < 0 Å and D Ly α > 50 Å, respectively.", '3.2. Line flux measurements': 'For each source in the JWST -PRIMAL sample, we measure the spectroscopic redshifts from the most prominent nebular emission lines, when detected, including: the [O ii ] λλ 3726 , 3729 and [O iii ] λλ 4960 , 5008 doublets, [Ne iii ] λ 3870, He i λ 3889, the [S ii ] λλ 6718 , 6732 doublet, and the Balmer lines H α , H β , H γ , and H δ . The [O ii ] λλ 3726 , 3729 and [S ii ] λλ 6718 , 6732 doublets are unresolved at all wavelengths in the NIRSpec Prism spectra, whereas the [O iii ] λλ 4960 , 5008 doublet is marginally resolved at z ≈ 6 and fully resolved at z ≳ 7 . 5 due to the increasing spectral resolution with wavelength in this particular configuration. H α and [N ii ] λ 6585 are unresolved at all redshifts. [Ne iii ] λ 3870 and He i λ 3889 are also generally unresolved, ex-  \nFig. 3 Three examples of normalized NIRSpec Prism spectra from the JWST -PRIMAL sample. The top panel shows an example of a strong DLA from local H i absorption, the middle panel a strong LAE, and the bottom panel a combination of both. The redshifts, D Ly α measurements, and D Ly α integration regions are highlighted for each case.  \n<!-- image -->  \nÅ  \ncept for the highest redshifts at z ≳ 10. Due to the broad wavelength coverage (0 . 6 -5 . 3 µ m), H α and [S ii ] λλ 6718 , 6732 can be detected up to z ≈ 7, H β and the [O iii ] λλ 4960 , 5008 doublet up to z ≈ 9 . 5, and Ly α can be detected at all redshifts from z ≳ 5 . 5. The auroral [O iii ] λ 4363 line is also in the majority of cases blended with H γ but becomes resolved at z ≳ 9 . 5, enabling robust measurements of this important transition (Heintz et al. 2023d; Hsiao et al. 2023). Unfortunately, the [O iii ] λλ 4960 , 5008 doublet is redshifted out of the NIRSpec Prism coverage at approximately the same redshifts, hindering direct Te -based metallicity measurements with the JWST / NIRSpec Prism configuration, except for a narrow redshift range ( z ≈ 9 . 3 -9 . 5).  \nTo determine the physical properties of the star-forming regions in each of the target galaxies, we measure or derive upper bounds on the line fluxes for each of these nebular emission lines when covered by the spectra. We fit the continuum around the lines with a simple polynomial and superimpose a  \nTable 2 Line flux measurements for the most prominent nebular emission lines JWST -PRIMAL sample.  \nNotes. All measurements are reported in units of 10 -20 erg s -1 cm -2 . Uncertainties on measurements are stated at 1 σ standard deviations and upper limits at 2 σ confidence limits. Column (1): Spectroscopic redshift. Column (2): Unresolved [O ii ] λ 3726 , 3729 doublet line flux. Column (3): Unresolved [Ne iii ] λ 3870 + He i λ 3889 line flux measurements, except for the highest redshift sources ( z > 10). Column (4): Balmer H δ . Column (5): Unresolved Balmer H γ + auroral [O iii ] λ 4363 line flux measurements, except for the highest redshift sources ( z > 9 . 5). Column (6): Balmer H β . Column (7): [O iii ] λ 5008 line flux, where [O iii ] λ 4960 is assumed to be 1 / 3 of this. Column (8): Unresolved Balmer H α + [N ii ] λ 6585. Column (9): Unresolved [S ii ] λλ 6718 , 6732 doublet line flux. A full version of this table, including an expanded set of emission line identifications and measurements, can be found online.  \nset of redshifted Gaussian line profiles at the rest-frame wavelength of each transition and fit for the redshift z spec, line equivalent widths (EWs) and fluxes. The line widths (full-width-athalf-maximum; FWHM) are in most cases directly proportional to the line-spread function of the NIRSpec Prism configuration (Jakobsen et al. 2022) and is modelled as such. For each case, we tie z spec and the FWHM corrected by the spectral resolution, e ff ectively assuming that the nebular emission lines all originate from the same H ii regions. An example of the line model fits is shown in Fig. 4. All the main identified lines, their derived line fluxes, or 2 σ upper limits, are summarized in Table 2 for the JWST -PRIMAL sample and provided in full online 8 , which also includes an expanded line list.', '3.3. [ O iii ]+H β equivalent widths': 'Asimple diagnostic of the specific star-formation rate (sSFR) of galaxies is the emission line [O iii ] + H β EWs. These have been found to increase with increasing redshifts (Khostovan et al. 2016; Endsley et al. 2021; Matthee et al. 2023), and toward lower stellar masses and metallicities (Malkan et al. 2017; Reddy et al. 2018b), due to the more intense ionization fields and active star formation in early galaxies. We derive the [O iii ] + H β EW for each of the sources in our sample that have these lines covered by the spectra (311 sources at z < 9 . 5). As demonstrated in Fig. 5, we integrate the flux density over the spectral region covering the H β and [O iii ] λλ 4960 , 5008 lines, normalized by the continuum flux, and correct the redshifted EWs to the rest frame. For each source, we extrapolate the continuum flux at the position of the lines by fitting a linear polynomial to the spectral regions on each side of the line complex.  \nIn Fig. 6, we show the [O iii ] + H β EWdistribution of the full JWST -PRIMAL sample and divided into redshift bins: z = 5 . 56.0, z = 6 . 0-7.0, z = 7 . 0-8.0, and z = 8 . 0-9.5. The median of the full distribution is 1045 Å, with 16th and 84th distribution percentiles from 356 Å to 2543 Å. We observe no appar-  \nevolution with redshift, but note that the number statistics is largely dominated by the lowest redshift bins. We also compare the full distribution to the much larger, albeit photometrically derived distribution from the CEERS survey of galaxies at z = 6 . 5-8.0 (Endsley et al. 2023). The two distributions are in good agreement, though the peak of the spectroscopic distribution is slightly shifted towards higher EWs. We note that the spectroscopically derived [O iii ] + H β EWs are much more accurate than photometric estimates, however, the latter potentially being underestimated by up to ≈ 30% (Duan et al. 2023). The fact that EW distribution of the JWST -PRIMAL sample is only slightly shifted, but otherwise follow the same distribution, indicates that there are no strong biases in our spectroscopic sample toward strong emission-line sources from the underlying galaxy population.', '3.4. Ultraviolet spectral slopes, magnitudes and luminosities': 'The UV spectral slope, β UV, where F λ ∝ λ β UV , of the stellar continuum at wavelengths λ rest = 1250-2600 Å encodes key information regarding dust attenuation (Calzetti et al. 1994; Meurer et al. 1999), stellar metallicity (Cullen et al. 2021), the average stellar population age (Zackrisson et al. 2011), and the escape fraction of ionizing photons (Chisholm et al. 2022) from galaxies. The minimum value expected for standard stellar populations with a Galactic IMF and with negligible dust attenuation is β UV ≈ -2 . 5 (Cullen et al. 2017; Reddy et al. 2018a). Steeper spectral power-law indices of β UV ≈ -3 . 0 would require a very recent onset of star formation ( ≲ 2 Myr) or near-pristine, neutral gas and negligible nebular processing. JWST has enabled measurements of β UV of galaxies at z > 10, hinting at even steeper spectral slopes, particularly for the faintest, less massive systems (Topping et al. 2023; Cullen et al. 2023b). These results have, however, been inferred from photometric data alone, which may be contaminated by varying Ly α or metal emission (or absorption) line strengths or even low-luminosity AGN.  \n<!-- image -->  \n<!-- image -->  \nFig. 4 Example of line flux modeling for one of the PRIMAL sources at z = 6 . 6345. The JWST / NIRSpec Prism spectrum is shown in black and the associated error spectrum in gray. The best fit continuum and Gaussian line model is shown in red.  \n<!-- image -->  \nHere, we measure β UV directly from the spectra for each of the JWST -PRIMAL sources at z = 5 . 5-13.4, carefully masking any identified emission lines or broad Ly α absorption troughs. We use the emcee minimizer within the LMFit framework and recover the parameter covariances to estimate the median and 16th to 84th percentiles of the posterior distributions on β UV and the normalization constant derived at UV 1550 Å rest-frame, m UV of the best-fit model, F λ = F UV , 1550 λ β . An example of one of the model fits is shown in Fig. 7. For our full sample, we recover UV continuum slopes ranging from β UV = -1 down to a few cases with β UV ≲ -3 . 0. This suggests that this set of rare galaxies are particularly dust- and metal-poor ( Z gas / Z ⊙ < 1%), have mean stellar populations with ages ≲ 2 Myr (e.g., Cullen et al. 2023a), and likely have ionizing photon escape fractions of ≳ 10% (Chisholm et al. 2022). In Fig. 8, we show the de-  \nFig. 5 Example of a [O iii ] + H β line equivalent width measurement. The blue-shaded region marks the lines integrated over, also represented by the top red bar. The continuum flux is measured by fitting a linear polynomial to the spectrum on each side of the lines.  \n<!-- image -->  \nFig. 6 [O iii ] + H β emission line equivalent width distribution (bottom) and normalized cumulative distribution function (top). The full sample up to z ≈ 9 . 5, where these lines can be detected, is shown by the black step function, and divided into redshift bins as indicated by the colors. The observed distribution of the spectroscopic JWST -PRIMAL sample is compared to the larger photometric CEERS survey results from Endsley et al. (2023).  \n<!-- image -->  \nal slopes as a function of redshift. We observe a marginal tendency for galaxies to have steeper UV slopes  \nFig. 7 Example of spectral fitting of the UV power-law index, β UV, where F λ ∝ λ β UV , of the stellar continuum. The grey curve shows the full NIRSpec Prism spectrum of a galaxy at z = 9 . 4383, with the error spectrum shown by the dotted line. The black part of the spectrum highlights the rest-frame fitting region from λ rest = 1250 -2600 Å. The top green bar marks the integration region for the M UV measurements, λ = 1400 -1700 Å in the rest frame.  \n<!-- image -->  \nat higher redshifts, though still consistent with a non-evolution within the scatter in the sample distributions. We further note that distribution of dust and H ii regions may account for the large scatter in the population around the expected slope for a standard stellar population with negligible dust and maximum nebular continuum contribution (e.g., Calzetti et al. 1994; Vijayan et al. 2024).  \nFig. 8 The galaxy UV spectral slope, β UV as a function of redshift. The red points mark the JWST -PRIMAL sample from this work, and the shown quantities are all derived directly from the NIRSpec Prism spectroscopy. The yellow hexagons show the mean β UV in bins denoted by the horizontal errorbars, and the vertical errorbars represent the 16th to 84th percentile of the distribution in the respective bins. The blue shaded region represents the minimum value β UV = -2 . 6 to -2 . 4 for a standard stellar population with negligible dust and maximum nebular continuum contribution (see Sect. 3.4 and Cullen et al. 2023a).  \n<!-- image -->  \nTo determine the UV luminosity, L UV, and absolute magnitude, M UV, for each source we derive the flux density and corresponding UV magnitude, m UV at rest-frame 1550 Å from the spectral power-law model. We convert the apparent m UV to the intrinsic, absolute brightness, M UV, via  \nM UV = m UV -µ ( z ) + 2 . 5 × log(1 + z ) (2)  \nwhere µ ( z ) is the distance modulus at a given redshift derived from the astropy.cosmology module. The derived M UV for each source are summarized in Table 3, not corrected for potential magnification factors for galaxies drawn from lensing cluster surveys. Our full sample spans absolute UV magnitudes of M UV = -22 to -16 mag, respectively (see Fig. 1).', '3.5. Star-formation rate': 'The Balmer recombination lines originating from star-forming H ii regions, H α in particular, are robust tracers of star formation on short ( ≲ 20 Myr) timescales. The star-formation rate (SFR) has been found to increase at a given stellar mass (Whitaker et al. 2012; Speagle et al. 2014; Salmon et al. 2015; Thorne et al. 2021; Sandles et al. 2022) and rest-frame UV size (van der Wel et al. 2014; Ward et al. 2023; Langeroodi & Hjorth 2023) for galaxies at increasing redshifts, following the evolution in the peak of the stellar mass halo mass relation with redshift (e.g., Behroozi et al. 2019). The SFR can now be measured directly from the Balmer recombination lines for galaxies at z > 6 with JWST , revealing intense star-forming galaxies during the reionization epoch (e.g., Heintz et al. 2023a; Shapley et al. 2023; Fujimoto et al. 2023a; Sanders et al. 2024; Nakajima et al. 2023; Curti et al. 2023b). Given that H α is only detected for the subset of the sample at z < 7, and may be contaminated by [N ii ], we instead derive the SFR based on H β for the majority of sources at z < 10 as:  \nSFRH ( M yr -1 ) = 5 5 × 10 -42 L (erg s -1 ) × f ( z 10) (3)  \nβ ⊙ . H β H α/ H β < following Kennicutt (1998), but here assuming the more topheavy initial mass function (IMF) from Kroupa (2001), which is likely more representative of the high-redshift galaxy population (e.g., Steinhardt et al. 2023). We derive the H β line luminosity via L H β = F H β × 4 π D 2 L ( z ) with DL ( z ) the luminosity distance at the given redshift and assume f H α/ H β = 2 . 86 from the Case B recombination scenario at Te = 10 4 K (Osterbrock & Ferland 2006). We propagate the uncertainty of ≈ 20-30% from the exact choice of the IMF in the relation, in addition to the dependence on the electron density and temperature of the H ii region which is only of the order ≲ 5%.  \nFor the sources at z > 10, where H β is redshifted out of the NIRSpec Prism spectral coverage, we instead rely on the [O ii ] λ 3727 line strength  \nSFR[OII]( M ⊙ yr -1 ) = 1 . 0 × 10 -41 L [OII](erg s -1 ) ( z > 10) (4)  \nagain following Kennicutt (1998), and assuming the Kroupa (2001) IMF. For the full sample, we derive SFRs in the range 0.1-100 M ⊙ yr -1 , not corrected for the relevant magnification factors or dust extinction, as summarized in Table 3. The SFRs can also be derived from overall UV luminosity of the sources following Madau & Dickinson (2014), where  \nSFRUV( M ⊙ yr -1 ) = 1 . 0 × 10 -28 L UV (erg s -1 Hz -1 ) (5)  \nassuming 10% solar metallicity which traces star formation on a longer timescale ( ∼ 100 Myr) compared to the Balmer line ( ∼ 100 Myr, see Calzetti 2013). We find that the SFRUV estimate is consistent with that derived from the Balmer recombination lines or [O ii ] for the full sample, when allowing for up to AV ≈ 0 . 5 mag of dust attenuation.  \nTable 3 Spectroscopically derived physical properties of the JWST -PRIMAL sample.  \nNotes. Column (1): Spectroscopic redshift. Column (2): Absolute UV (rest-frame ≈ 1500 Å) magnitude. Column (3): Rest-frame [O iii ] + H β equivalent width. Column (4): Rest-frame UV (1250 -2600 Å) spectral slope. Column (5): The Ly α damping parameter. Column (6): The star-formation rate derived from H β ( z < 10) or [O ii ] ( z > 10). Column (7): The [O iii ] λ 5008 / [O ii ] λ 3727 line flux ratio. Column (8): Gas-phase metallicity inferred from the joint PDF of the di ff erent strong-line calibrations from Sanders et al. (2024). A full version of this table can be found online.', '3.6. Gas-phase metallicity': 'One of the most fundamental characteristics of galaxies is their metallicity. This is notoriously di ffi cult to measure directly, in particular for the highest redshift galaxies, motivating the use of strong-line diagnostics of the most prominent nebular emission lines (see e.g. Kewley et al. 2019; Maiolino & Mannucci 2019, for recent reviews). In rare cases, the auroral [O iii ] λ 4363 emission line can be detected directly, which enables Te sensitive measurements of the gas-phase metallicity (known as the direct Te -method), typically quantified as the oxygen abundance, 12 + log (O / H) (e.g., Izotov et al. 2006). This particular feature has now been detected in star-forming galaxies at z > 6 with JWST (e.g., Schaerer et al. 2022; Curti et al. 2023a; Heintz et al. 2023c; Nakajima et al. 2023; Sanders et al. 2024), even out to z ≈ 10 (Hsiao et al. 2023; Heintz et al. 2023a,d). The auroral [O iii ] λ 4363 line feature, however, will only be resolved from H γ in the minority of sources due to the increasing spectral resolution towards the red end of the spectra.  \nWe therefore base our main metallicity measurements on strong-line calibrations. To preserve the homogeneity of the measurements, we use a single set of line diagnostics from Sanders et al. (2024) for the di ff erent line ratios, derived for galaxies at z = 2-9 through the direct Te -method. Specifically, we consider the strong-line ratios: O3 = [O iii ] λ 5008 / H β , O2 = [O ii ] λ 3727 / H β , R23 = ([O iii ] λλ 4960 , 5008 + [O ii ] λ 3727) / H β , Ne3O2 = [Ne iii ] λ 3870 / [O ii ] λ 3727, and O32 = [O iii ] λ 5008 / [O ii ] λ 3727, which are available at most of the targeted redshifts. Based on the line-flux measurements and errors (representing the likelihood distributions of fluxes), we construct the probability density function (PDF) for each diagnostic ratio where the width denotes the scatter in the relation. We then combine the PDFs from the available calibrations, inversely weighted by their scatter, for each individual source to determine the gas-phase metallicity through the median and 16th to 84th percentiles of the joint PDF, as illustrated in Fig. 9. This more conservative approach takes into account the added uncertainty from sources with potential high ionization parameters, typically quantified via O 32, thus minimizing the uncertainty from assuming a single line-ratio.  \nFig. 9 Example of joint gas-phase metallicity estimate for one of the sources. The blue-colored function shows the normalized joint probability density function (PDF) based on the normalized PDFs of the individual line-ratio diagnostics from Sanders et al. (2024). The blue dashed line marks the median of the joint PDF.  \n<!-- image -->  \nWe derive gas-phase metallicities ranging from 12 + log (O / H) = 6 . 5 to 8 . 2, corresponding to 0.6% to 30% of solar abundance, respectively (given 12 + log (O / H) ⊙ = 8 . 69; Asplund et al. 2009). The relevant line fluxes used for these measurements are provided in Table 2 and the derived metallicities are listed for each source in Table 3.', '4. Lymanα emission and absorption in galaxies during the reionization epoch': 'With the spectroscopically-derived physical properties for the full set of JWST -PRIMAL sample sources we can now inves-  \nFig. 10 Ly α damping parameter, D Ly α , for the full JWST -PRIMAL sample as a function of redshift. The corresponding age of the Universe assuming the cosmological parameters from Planck Collaboration et al. (2020) is given at the top. The colored backgrounds indicate the typical physical regimes probed for a given D Ly α : Strong Ly α emission ( D Ly α < 0 Å; blue), extended ionized bubbles or weak Ly α emission ( D Ly α = 0-35 Å; orange), IGM absorption from x HI = 0-1 ( D Ly α = 35-50 Å; green), and strong galaxy integrated DLAs with N HI > 10 21 cm -2 from local H i gas ( D Ly α > 50 Å; red).  \n<!-- image -->  \ntigate and chart the redshift evolution and the physical driver of strong galaxy DLAs produced by massive H i gas reservoirs to prominent ionized bubbles and Ly α emission during the reionization epoch.', '4.1. Charting D Ly α as a function of redshift': 'To obtain a first overview of the prevalence and prominence of DLAs and LAEs across time in the reionization epoch, we show D Ly α as a function of redshift for the full sample in Fig. 10. This represents the evolution from z = 13 . 4 to z = 5 . 5, corresponding to approx. 300 Myr to 1 Gyr after the Big Bang. The D Ly α parameter space has been color-coded to highlight the regimes driven by the physical models for the expected dominant contribution to the observed Ly α transmission curve. The number statistics of sources located in each region are summarized in Table 4, divided into redshift bins of z = 5 . 5-6.0 (159 sources), z = 6 . 0-7.0 (169), z = 7 . 0-8.0 (103), z = 8 . 0-10.0 (46), and z = 10 . 0-13.4 (17).  \nWe find that there are 38 strong LAEs in the full JWST -PRIMAL sample (i.e. 8%), based on D Ly α < 0 Å. None of these are found at z > 9. At z < 8, we find that 38 out of 494 of sources show strong LAEs (8%) in approximately equal fractions from z = 5 . 5-6.0, z = 6 . 0-7.0, and z = 7 . 0-8.0, albeit with the most prominent LAEs quantified by the lowest D Ly α at z < 7. This fraction of LAEs could be as high as ≈ 20% if we include all sources with excess Ly α flux as potential weakly-emitting LAEs. These results are slightly lower than found previous lit-  \ne studies, showing a typical fraction of LAEs of ≈ 30% at z ≈ 6-8 (e.g., Saxena et al. 2023b; Nakane et al. 2023; Simmonds et al. 2023; Witstok et al. 2023; Jung et al. 2023; Jones et al. 2023; Witten et al. 2024). However, they were also found to be much more sparsely populated beyond z ≳ 8 (though see Tang et al. 2023; Bunker et al. 2023b; Fujimoto et al. 2023b, Witstok et al. in prep.), consistent with our results, and as also recovered by simulations (Hutter et al. 2023). These findings would also naturally explain the null detections of LAEs at z = 8 . 8 in large photometric surveys using narrow-band imaging such as the UltraVISTA (e.g., Laursen et al. 2019).  \nSimilarly, we consider the redshift evolution of strong DLAs, here classified as sources with D Ly α > 50 Å, corresponding to H i column densities N HI ≳ 10 21 cm -2 . We find that the fraction of galaxy DLAs only slightly increases from ≈ 60% at z ≈ 6 up to ≈ 65-90% at z > 8. A similar conclusion is reached by Umeda et al. (2023), who find an ubiquitous presence of extremely strong DLAs ( N HI ≳ 10 22 cm -2 ) at z > 10. This evolution is more clearly highlighted in Fig. 11 where we show the histograms of D Ly α for the five di ff erent redshift ranges, demonstrating higher prevalence of galaxy-integrated DLAs at increasing redshifts. These results are in sharp contrast to previous studies of Lyman-break galaxies near the peak of cosmic star formation at z ∼ 2-3 (Pettini et al. 2000; Steidel et al. 2010), where only a small subset of the population ( < 25%) show integrated galaxy line-of-sight H i column densities of N HI ≈ 3 × 10 20 cm -2 (Shapley et al. 2003). This corroborates our findings here of an increasing probability of having high H i column densities or  \nTable 4 Ly α absorption vs. emission statistics.  \nFig. 11 Histogram of the Ly α damping parameter, D Ly α , distribution (bottom) and normalized cumulative distribution function (top). The full JWST -PRIMAL sample is represented by the black step function, and divided into sub-bins according to their redshift as indicated by the colors. Strong DLAs are more prevalent in galaxy spectra at increasing redshift.  \n<!-- image -->  \nlarger fraction of DLAs at higher redshifts, down to z ≈ 2-3. This may indicate that we are probing the epoch of first-galaxy assembly with excessive amounts of accreting or overdense pristine H i gas that are yet to be processed into stars or ionized by the first stellar sources. Indeed, the average gas mass density of galaxies have been inferred to increase by an order of magnitude through accretion from z = 10 to 6 (see ? Heintz et al. 2022, and Sect. 5 for further discussion). (see Sect. 5 for further discussion). Fig. 11 also demonstrates an excess of galaxies with strong LAEs ( D Ly α < 00 Å) at z = 5 . 5-6.0 from the full distribution, distinct for the epoch when the universe is fully ionized.  \nWe note that the fraction of ≳ 70% galaxy DLAs at z > 9 may be a lower bound, given that our selection criteria on the S / NLy α will preferentially be biased towards strong Ly α emission. On the other hand, fainter, less massive sources are likely to have stronger Ly α EWs (Saxena et al. 2023a) and contribute to most of the ionizing photons (Lin et al. 2024). To investigate this, we compare D Ly α to the spectroscopically-derived M UV in Fig. 12. We confirm the trend of the intrinsically faintest  \nFig. 12 The absolute UV brightness, M UV, of the JWST -PRIMAL sources as a function of D Ly α . The sources with the highest Ly α EWs (negative D Ly α ) are observed to represent the UV faintest galaxy population at z > 5 . 5.  \n<!-- image -->  \nsources showing the highest Ly α EWs, with the strongest LAEs with D Ly α < -100 Å showing luminosities predominantly at M UV > -19 mag. This further suggests that the large H i gas overdensities are preferentially located close to bright, massive sources. In the following sections, we will further investigate the underlying physical properties driving the prevalence and prominence of strong Ly α emission and absorption in these reionization era sources.', '4.2. The role of stellar ionization fields on Ly α emission and absorption': "The intense UV radiation fields produced by younger, more massive stars are likely to play a key role in the ionization of H i and the production and escape of Ly α photons. To test this and quantify it in terms of our observations, we compare D Ly α to two common tracers of active star formation and high ionization: The [O iii ] + H β EW, and the ionization parameter O 32 = [O iii ] λ 5008 / [O ii ] λ 3727. We find no apparent correlations between D Ly α and any of the two quantities. For example, strong galaxy DLAs are observed in galaxies spanning the entire dynamical range of [O iii ] + H β EWs and O 32 ionization. Similarly, strong LAEs are observed in galaxies spanning [O iii ] + H β EWs of 100 -10 4 Å, and are not directly associated with sources with high or low O 32 ratios. This is at odds with studies of local 'green pea' galaxies at z ≈ 0, believed to be robust analogs to highz star-forming galaxies, that find that larger Ly α emission EWs and escape fractions are commensurate with higher, more intense stellar ionization fields (Yang et al. 2017; Jaskot et al. 2019; Hayes et al.  \n<!-- image -->  \nFig. 13 Correlations of D Ly α with the galaxy properties that represent the average the specific SFR ([O iii ] + H β EW, left), and the intensity of the stellar ionization field ( O 32 = [O iii ] λ 5008 / [O ii ] λ 3727, right). The square symbols show the JWST -PRIMAL sources at z < 9 . 5, where [O iii ] can be detected, and a color-coded according to their redshift (identical to Fig. 12).  \n<!-- image -->  \n2023). We note, however, that the galaxies at z > 8 predominantly show high D Ly α absorption and O 32 ionization and instead typically represent the large-scale environment or H i gas covering fraction at these redshifts. Consequently, our results disfavor a scenario where the H i abundance (or the lack thereof) is driven by the ionization output of the stellar sources in the galaxy.", '4.3. Gas-phase metallicity and Ly α absorption and emission': 'The gas-phase metallicity of galaxies is one of the key regulators of the gas cooling e ffi ciency, the ionization state of the gas, and the escape fraction of Ly α and LyC photons. Lower metallicities are often associated with higher ionization and temperatures, effectively reducing the observed line-of-sight H i column density. By contrast, high metallicity and likely thereby also high dust content, will enable less e ffi cient ionization of H i and as a consequence more e ffi ciently scatter the output ionizing photons, effectively suppressing the intrinsic and observed escape fraction of Ly α photons (Dayal et al. 2011; Henry et al. 2015; Dijkstra et al. 2016; Verhamme et al. 2017; Yang et al. 2017; Izotov et al. 2018; Jaskot et al. 2019). We therefore expect a strong correlation between the gas-phase metallicity and D Ly α for our highredshift sample sources.  \nTo test this scenario, we show D Ly α as a function of gasphase metallicity, 12 + log(O / H), for each of the sources in the JWST -PRIMAL sample in Fig. 14. The strongest LAEs ( D Ly α < 0 Å) are predominantly observed in galaxies with low metallicities, and ubiquitously found at 12 + log (O / H) < 7 . 7, corresponding to < 10% solar metallicity. These results are consistent with the expectations that prominent Ly α emission can only escape from low-metallicity systems. By contrast, we find that galaxies with strong DLAs ( D Ly α > 50 Å) span the entire dynamical range of inferred metallicities, 12 + log (O / H) = 6 . 5-8 . 2 (i.e. 0.6-30% solar). The most prominent DLA systems are found at lower metallicities at any redshifts z ≳ 7, supporting the scenario where these mostly trace young, less chemically enriched galaxies.  \nFig. 14 The gas-phase metallicity, 12 + log(O / H), as a function of the Ly α damping parameter, D Ly α . The strength of the Ly α absorption is observed to increase with increasing metallicity, whereas Ly α emission is most e ffi ciently produced or allowed to escape in low-metallicity galaxies.  \n<!-- image -->', '4.4. The impact of galaxy DLAs on the reionization history': 'One of the major implications of high H i -cloumn density DLAs observed in galaxy-integrated spectra at z ≳ 6 is their e ff ect on the escape fraction of Ly α and LyC photons, key to constrain the reionization history and topology. Numerical modelling and reionization simulations have previously only taken IGM absorption into account when predicting the radiative transfer and transmission of Ly α (e.g., Laursen et al. 2019; Hutter et al. 2023). These overabundant local H i gas reservoirs are likely to be the main component in the transfer of ionizing photons, at least in the early evolutionary stages of galaxies, so quantifying this additional local ISM or CGM absorption is vital to draw ro-  \n<!-- image -->  \nFig. 15 Comparison of D Ly α to the spectral β UV slope (left) and the ionizing photon production e ffi ciency, ξ ion (right). The JWST -PRIMAL sources are shown by the red squares and color-coded according to redshift (identical to Fig. 12). The light- to dark-blue colored regions in the left panel indicates increasing escape fractions of LyC photons, f Lyc esc , based on the low-redshift empirical relation between β UV and f Lyc esc from Chisholm et al. (2022). The grey-colored band in the right panel shows the canonical reionization values log( ξ ion / Hzerg -1 ) = 25 . 1 -25 . 3 (Robertson et al. 2013). Strong LAEs are observed with a large variety of β UV slopes, but are predominantly associated with galaxies with high ionizing photon production e ffi ciencies. Galaxies with strong DLAs span the entire probed range in ξ ion.  \n<!-- image -->  \nconclusions and physically interpret the output and impact of galaxy-wide DLAs in simulations.  \nTo gauge this e ff ect, we first compare D Ly α to the β UV spectral slopes of each galaxy in Fig. 15. β UV has been found to be strongly correlated with the escape fraction of Lyman Continuum (LyC) photons, f LyC esc , in low-redshift galaxies sample (Chisholm et al. 2022; Flury et al. 2022). We color-grade the β UV parameter space with the empirical relation presented in their paper Chisholm et al. (2022, their Eq. 11), highlighting the inferred escape fractions from f LyC esc = 5% to 50%. We find no apparent correlations between β UV or the inferred f LyC esc with the strength of the Ly α emission or DLA absorption. This indicates that the hardness of the UV spectral shape is not the main driver of the ionization of H i , and that larger H i abundances do not translate directly into lower escape fractions of ionizing photons. This is somewhat at odds with theoretical predictions and local observations, since both Ly α and LyC photons are a ff ected by resonant scattered from abundant interstellar or circumgalactic H i gas columns (e.g., Flury et al. 2022; Xu et al. 2022; Begley et al. 2024). We caution, however, that the relations for f LyC esc derived from the local galaxy sample at z ≈ 0 . 3 studied by Chisholm et al. (2022) may not be universally representative for galaxies at z ≳ 6 (Choustikov et al. 2024; Pahl et al. 2024). Since the highz galaxy population show even more abundant H i gas reservoirs than galaxies at equivalent metallicities or O 32 ionization parameters (Heintz et al. 2023d), the escape fraction of Ly α and LyC photons is likely correspondingly lower. Indeed, comparing the sources from Chisholm et al. (2022) that are mostly absorbed by the neutral H i gas, the inferred LyC escape fraction could be a ≈ 10 × lower.  \nTo more directly relate the amount of photons capable of ionizing H i to the presence and strength of Ly α emission or DLAs, we also compare the ionizing photon production e ffi ciency, ξ ion, to D Ly α in Fig. 15. ξ ion can be derived through the Balmer re-  \nination line strengths, such as H β , since these accurately trace the ionizing photon production rate from massive stars on short ( ≲ 10 Myr) timescales (Bouwens et al. 2016). Following Matthee et al. (2023), we define  \nξ ion (Hz erg -1 ) = L H β (erg s -1 ) c H β (erg) L UV(erg s -1 Hz -1 ) (6)  \nwhere c H β = 4 . 86 × 10 -13 erg is the H β line-emission coe ffi cient, assuming a Case B recombination scenario with Te = 10 4 K (e.g., Schaerer 2003). We find from Fig. 15 that galaxies with strong DLAs ( D Ly α > 50 Å) occupy most of the probed ξ ion parameter space, though with a median of log( ξ ion / Hzerg -1 ) = 25 . 55, higher than the canonical values, log( ξ ion / Hzerg -1 ) = 25 . 1 -25 . 3 (Robertson et al. 2013). The large scatter in ξ ion at D Ly α > 50 Å could indicate bursty star formation, particularly since ξ ion traces massive O-stars, and may reflect typical younger stellar populations that are yet to ionize the neutral H i gas in their local surroundings. By contrast, galaxies that are strong LAEs tend to exhibit high ionizing photon production efficiencies, occupying the region with log( ξ ion / Hzerg -1 ) > 25 . 5, higher than the canonical value. Since these are predominantly found at z < 8, we argue that it is mainly this particular subset of galaxies with high ξ ion and late in the reionization epoch that provides the strongest contribution to this last phase transition.', '5. Discussion': 'To place our results into context of the large-scale e ff orts to quantify the galaxy mass assembly and reionization timeline, we describe in this section a physical scenario connecting the observed prevalence and prominence of Ly α emission and DLA absorption from abundant, pristine H i gas in galaxies through redshift and as a function of their intrinsic physical properties.  \nWe further quantify the e ff ect on D Ly α from the neutral hydrogen fraction based on more sophisticated IGM simulations and estimate the e ff ect of strong galaxy DLAs on photometric redshifts measurements at z > 8, which appear to be systematically overestimated.', '5.1. H i gas assembly and reionization across cosmic time': 'In this work, we quantified the strength of DLA absorption (Ly α emission) via the positive (negative) strength of D Ly α . We found that galaxies with high D Ly α , indicating substantial H i gas column densities in the ISM or CGM of galaxies, are ubiquitously observed at z > 9 and in galaxies spanning the entire dynamical range of physical properties probed here. The fact that there appears to be no correlation between the amount of local H i gas (or the lack thereof) and the ionizing output of the galaxies themselves may suggest that DLAs trace particularly young galaxies, at all redshifts, that are yet to ionize the surrounding H i gas or convert it into molecules and eventually stars. For galaxies with D Ly α > 50 Å, equivalent to N HI > 10 21 cm -2 , the absorbing H i gas is completely self-shielding of ionizing photons, even with substantial clumping with density variations down ∆ N HI ≈ 10 -3 supporting this scenario.  \nOn the contrary, strong LAEs are only observed at z ≲ 8. This indicates that only after z ≈ 8 -9, on average, is the surrounding large-scale IGM su ffi ciently ionized for a prominent fraction of Ly α photons to escape. Indeed, the observed redshift evolution of D Ly α is in good qualitative agreement with the evolution of the Ly α escape fraction, as compiled by Saxena et al. (2023a); Hayes & Scarlata (2023). Consequently, this favors more late and rapid reionization models (e.g., Naidu et al. 2020), over a smooth, early transition (e.g., Finkelstein et al. 2019). A similar rapid evolution seems to be suggested by damped Ly α wings in quasar spectra (Greig et al. 2017, 2019; Davies et al. 2018; Bañados et al. 2018; Yang et al. 2020; Wang et al. 2021; Ocvirk et al. 2021; ˇ Durovˇcíková et al. 2020, 2024) and Ly α emission statistics of galaxies (Ouchi et al. 2010; Sobacchi & Mesinger 2015; Mason et al. 2018, 2019), indicating a near-fully neutral IGM at z ≈ 8.  \nThe strong LAEs are predominantly associated with UV faint ( M UV > -19 mag) galaxies with high ionizing photon production e ffi ciencies, log( ξ ion / Hzerg -1 ) > 25 . 5. This could either indicate that these particular systems were the most e ffi -cient at ionizing their immediate surroundings, or that viewing angles play a larger e ff ect at low ( z < 8) redshifts compared to their higher redshift counterparts. In this scenario, galaxies at z > 8 are deeply embedded into dense, neutral gas overdensities with large covering fraction, showing strong absorption independent of the viewing angle. At lower redshifts, low-density N HI channels are carved by ionizing photons and we thus observe a broader span in D Ly α at z < 8. This scenario would be commensurate with the non-dependence of the strength of LAEs with the O 32 ionization parameter, [O iii ] + H β EW, the β UV spectral slope and the inferred escape fraction of LyC photons. In fact, strongly directionally-dependent escape of ionizing radiation has been proposed by simulations (Cen & Kimm 2015; Rosdahl et al. 2022; Yeh et al. 2023). Since only ≈ 10 -20% of the galaxies at z < 8 have strong LAEs, we argue that it is either this particular subset of faint, low mass galaxies that drive reionization or that ionizing photons from galaxies at this epoch will only be able to escape through the ≲ 20% low-density sightlines.  \nFig. 16 The impact of considering residual neutral hydrogen on IGM transmission curves for models with ionized bubble sizes of 1 cMpc (blue), 10 cMpc (orange) and 50 cMpc (red) and a volume-averaged IGM neutral fraction x HI = 0 . 1 at redshift 9. The vertical black dashed line marks the observed Ly α wavelength at this redshift. The left column shows IGM damping wings calculated using the Miralda-Escudé (1998) model, assuming that the bubbles are completely ionized. The right column shows the e ff ect of accounting for a small amount of residual neutral gas in the bubbles, which can completely saturate the absorption in the Ly α forest. The top panel shows the models before convolving with the NIRSpec PRISM instrumental profile and the bottom panel shows the models after convolution.  \n<!-- image -->', '5.2. Insights into the D Ly α parameter from simulations': 'So far, the IGM regime of D Ly α has been investigated assuming a homogeneously ionized IGM with an average value of the ionized fraction x HI everywhere. It is therefore interesting to investigate how the D Ly α parameter evolves in simulations of patchy reionization, to better understand the allowed range of D Ly α and what some limitations may be.  \nAttenuation by the intergalactic medium is expected to play a significant role in regulating visibility of Ly α emission until the Universe was fully reionized (McQuinn et al. 2007), and will impact the D Ly α parameter to some extent over all of the regimes identified in Fig. 10. A full comparison of observed and simulated D Ly α across the whole parameter space would therefore require detailed modelling of the redshift evolution of ionized bubbles, as well as the detailed shape and of the Ly α emission profile (e.g., Hayes & Scarlata 2023). We instead focus here on identifying the allowed range of D Ly α for which attenuation by the IGM dominates the measurement. We measure D Ly α from IGM transmission spectra normalised by the continuum, i.e., with no contribution from Ly α emission.  \nFirst, we show how the Ly α transmission curves change based on the residual fraction of H i inside an ionized bubble. As discussed in several works, even a small amount of neutral gas can completely saturate absorption of the Ly α forest (Mesinger & Haiman 2007; Mason et al. 2018; Keating et al. 2023a). The  \nFig. 17 Redshift evolution of the damping parameter for continuum-normalised spectra dominated by IGM absorption. The grey shaded region shows the range of damping parameters calculated from a step function at the observed Ly α wavelength (which sets the lower limit) and the IGM damping wing calculated using the Miralda-Escudé (1998) model in a completely neutral Universe (which sets the upper limit). The coloured violin plots represent the distribution of damping parameters measured from IGM damping wings generated from the SherwoodRelics simulation of an inhomogeneously reionized IGM as presented in Keating et al. (2023a). The volume-averaged H i fraction x HI for each simulation snapshot is indicated in the legend.  \n<!-- image -->  \ndi ff erence in accounting for this residual neutral gas is shown in Fig. 16, which shows IGM transmission curves calculated using the Miralda-Escudé (1998) analytic model. In the top left panel, the bubbles are assumed to be completely ionized, allowing for transmission on the blue side of Ly α . When this neutral gas is accounted for, a sharp cuto ff in IGM transmission blueward of the Ly α emission is observed. The bottom panel shows the same transmission profiles, but now convolved with the NIRSpec PRISM instrumental profile. It is di ffi cult to statistically disentangle the di ff erence between the models with di ff erent bubble sizes once the residual H i is included.  \nBased on these results, we estimate an allowed range of IGM-dominated D Ly α parameters, shown in Fig. 17. The lower value of D Ly α is given by a step function, where there is no transmission blueward of Ly α and complete transmission redward of Ly α ˙ This mimics the limit of large bubbles containing a small amount of residual H i . The upper value of D Ly α is calculated from the Miralda-Escudé (1998) analytic model assuming the Universe is completely neutral and assuming a bubble size R b = 0 cMpc. We see that the maximum value of D Ly α grows with redshift, due to the evolution of the mean density of the Universe (Keating et al. 2023b).  \nFor comparison, we also measure IGM-dominated D Ly α parameters directly from damping wings constructed from a simulation of inhomogeneous reionization. We use the spectra presented and described in Keating et al. (2023a), which were generated from sightlines through a simulation from the SherwoodRelics simulation suite (Puchwein et al. 2023). We analyse results from five snapshots from this simulation at z = 6 -10. The  \ndistributions of D Ly α parameters we measure are shown as the violin plots in Fig. 17. We find that the range of allowed IGM D Ly α parameters is well captured by the grey shaded region we defined. At z = 6, some sightlines have smaller D Ly α parameters as the Ly α forest starts to show transmission at this redshift. At z = 10, some sightlines have larger D Ly α parameters due to effects not captured in the analytic model for IGM damping wings, such as dense clumps of gas close to the host halo or the e ff ect of infalling material. In general, however, the range of D Ly α parameters for which the reionization of the IGM plays a dominant role are in good agreement with the range of D Ly α = 35-50 Å as defined in Section 3.1. However, we note that the IGM will also play an important role in regulating the fraction of flux transmitted by galaxies showing evidence of Ly α emission, and hence falling into the regime of smaller values of D Ly α .', '5.3. The effect of Ly α damping wings on highz photometric redshift estimates': "Some of the first e ff orts to photometrically identify high-redshift galaxy candidates seemed to generally overpredict their redshifts when compared to the subsequent spectroscopic redshift determination of the sources at z > 8 based on strong nebular emission lines (e.g., Finkelstein et al. 2023b; Fujimoto et al. 2023a; Hainline et al. 2023; Steinhardt et al. 2023; Serjeant & Bakx 2023). In Heintz et al. (2023d), we proposed that this may due to the strong damping wings observed in the large fraction fraction of galaxies at these redshifts, suppressing the derived flux density in the blue-most photometric bands and thereby mimicking a slightly larger Ly α 'break redshift' (e.g., Arrabal Haro et al. 2023).  \nTo test this, we show the derived photometric redshifts for the subset of sources where these are available in the DJA (Valentino et al. 2023) 9 relative to the spectroscopically derived redshifts in Fig. 18. We find that the photometric redshifts are systematically overestimated for galaxies at z > 8, with median z phot -z spec = 0 . 11. Considering only the galaxies with strong DLAs ( D Ly α > 50 Å) the discrepancy become slightly larger, with median z phot -z spec = 0 . 15. To optimize future photometric redshift estimates, the e ff ects of strong Ly α damping wings needs to be incorporated into spectral energy distribution fitting algorithms like E azy (Brammer et al. 2008).", '6. Summary & Future Outlook': "In this work, we have compiled and systematically characterized 494 sources observed with JWST / NIRSpec during the reionization epoch, at z = 5 . 5-13.4, introducing the JWST -PRIMAL legacy survey. The main goal was to statistically quantify the presence and redshift evolution of strong damped Ly α absorption in galaxies, signifying abundant local H i gas reservoirs. To disentangle these systems from sources with damping parameters probing mainly a largely neutral IGM, or residing in extended ionized bubbles, we defined a new simple diagnostic, the Ly α damping parameter, D Ly α . This allowed us to statistically chart the H i gas assembly history, the onset and prevalence of strong Ly α -emitting galaxies and to explore the main underlying physical properties driving the emission and escape of ionizing and Ly α photons.  \nWe found that the majority ( ≈ 65 -90%) of galaxies at z > 8 are dominated by large H i gas reservoirs in their local  \nFig. 18 Photometric redshift, z phot, estimates relative to the spectroscopic redshifts, z spec, for the JWST -PRIMAL sample at z > 8. The sources are color-coded as a function of D Ly α , as indicated by the colorbar. The photometric redshifts systematically overpredicts the actual spectroscopically derived redshifts with median z phot -z spec = 0 . 11 ( z phot -z spec = 0 . 15 for D Ly α > 50 Å).  \n<!-- image -->  \nsurroundings (ISM or CGM) with galaxy-integrated H i column densities, N HI > 10 21 cm -2 . These strong galaxy DLAs appear to exist throughout the reionization epoch, even though the fraction was observed to decrease to ≈ 60% by z = 6 of the overall galaxy population. Similar H i column densities have been theoretically predicted from cosmological simulations (e.g., Finlator et al. 2018; D'Odorico et al. 2018; Pallottini et al. 2022). Further, we found that these strong DLA galaxies represent a large variety of intrinsic physical properties in terms of their UV luminosity, gas-phase metallicity 12 + log(O / H), UV spectral slope β UV, nebular [O iii ] + H β line EWs, O 32 ionization parameters, and ionizing photon production e ffi ciency, ξ ion. We thus surmised that strong galaxy DLAs are a sign of early galaxies in the process of formation, actively accreting and assembling large amount of primordial H i gas, that are yet to be ionized or processed into molecular gas and stars.  \nFrom the measured redshift evolution of D Ly α we were able to quantify the onset and increasing strength of Ly α , enabling us to decode the reionization history directly. We found that at z ≈ 8-10 galaxy spectra mainly probe a substantially neutral CGM, whereas evidence for low Ly α escape or significant ionized bubbles around these galaxies were only observed at later cosmic times ( z < 8). The fraction of LAEs was estimated to be ≈ 20% at z = 6. The strongest LAEs at z > 7 have been found to reside in particular overdense or ionized regions (Saxena et al. 2023b; Witstok et al. 2023; Witten et al. 2024). From theoretical considerations and cosmological simulations (Pallottini et al. 2022), we expect most of the galaxies studied here to be formed within the first 200-300 Myr of cosmic time at z ≈ 18-14  \nand rapidly build up their stellar mass and on-going star formation via substantial H i gas accretion and assembly. Over the next 100-300 Myr, down to z = 8 -100,  \nThe emergence of these sources showing strong 'roll-overs' near the Ly α edge in the first era of JWST / NIRSpec Prism spectra have important practical implications for future studies and observations: First, they complicate statistical inferences of the reionization history by tracing the increasing neutral fraction of the IGM from galaxy damping wing measurements (Keating et al. 2023a). Here, we provided the first step forward to statistically model the galaxy DLA population which will eventually be possible to remove to uncover the underlying IGM damping e ff ects. Second, the first e ff orts in photometrically identifying high-redshift galaxy candidates seemed to show a systematic positive o ff set in z phot -z spec towards higher redshifts (Finkelstein et al. 2023b; Fujimoto et al. 2023a; Hainline et al. 2023; Steinhardt et al. 2023; Serjeant & Bakx 2023). This may naturally be explained by the increasing fraction of galaxies with stronger Ly α damping wings at higher redshifts. This we quantified here and estimated a mean of 0.2 dex overestimations of the photometric redshifts of galaxies at z > 8. This e ff ect would need to be incorporated into photoz codes to optimize the accuracy on the z phot priors for the earliest galaxy population.  \nThe next natural avenues to explore the early galaxy DLAs in more depth are to examine the e ff ect of clustering, their correlations with physical sizes and SFR surface densities, and the H i gas content; a potential cause of the observed o ff set in the otherwise fundamental-metallicity relation (FMR). Similar to how strong LAEs are typically found to reside in high-density groups of galaxies, galaxy DLAs may trace overdense regions via their early assembly of H i . Moreover, the subset of galaxies showing strong DLAs integrated over their relative small UV sizes indicate exceptional dense gas surface densities (Kennicutt & Evans 2012). This, in combination with the inferred SFRs, may be key to understand the intense star formation and bright UV luminosities of early z > 9 galaxies (e.g., Finkelstein et al. 2023a; Franco et al. 2023; Harikane et al. 2023; Adams et al. 2023; Casey et al. 2023; Bouwens et al. 2023; Chemerynska et al. 2023; Mason et al. 2023; McLeod et al. 2024). Finally, the discovered o ff set from the FMR towards lower metallicities of galaxies at z = 710 have been interpreted as evidence for excessive pristine H i gas inflow (Heintz et al. 2023a), which have later been confirmed with larger galaxy samples, though the exact redshift for this transition is still debated (e.g., Curti et al. 2023b; Nakajima et al. 2023). This scenario will be possible to test directly by correlating the H i column densities of galaxies to their o ff set from the FMR, moving the 3D plane of galaxy properties to 4 dimensions. Obtaining higher-resolution spectroscopic data of these sources will further help quantifying the H i covering fraction and disentangle the e ff ects on the Ly α transmission curve from weak Ly α emission to extended ionized bubbles and expand the probed physical properties to include the electron densities and temperatures of the gas.  \nWhile previous e ff orts with ALMA to quantify the H i gas mass content using far-infrared tracers such as [C ii ] -158 µ mand [O i ] -63 , 145 µ m have been successfully applied to galaxies in the epoch of reionization (Heintz et al. 2021, 2022; Vizgan et al. 2022; Liang et al. 2024; Wilson et al. 2023), the results presented here provide the most direct and statistical census of the abundance and assembly of the primordial H i gas in early Universe. This naturally paves the way and provide an important benchmark for next-generation radio facilities, targeting H i of the first galaxies through the 21-cm transition, such as the square kilometre array (SKA; Dewdney et al. 2009).", 'Acknowledgements': 'We would like to thank Peter Jakobsen for his vision and heroic endeavour in optimally designing the JWST / NIRSpec instrument and some of its first on-sky observations and for enlightening discussions about the intricacies of the NIRSpec data. Further, we would like to thank John Chisholm for helpful clarifications and discussions related to the escape fraction of ionizing photons and Aayush Saxena for enlightening conversations on the escape and absorption of Lymanα photons.  \nThis work has received funding from the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number MB22.00072. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant DNRF140. The data products presented herein were retrieved from the DAWN JWST Archive (DJA). DJA is an initiative of the Cosmic Dawn Center, which is funded by the Danish National Research Foundation under grant DNRF140. P.D. acknowledge support from the NWO grant 016.VIDI.189.162 (\'ODIN") and warmly thanks the European Commission\'s and University of Groningen\'s CO-FUND Rosalind Franklin program. Support from the ERC Advanced Grant INTERSTELLAR H2020 / 740120 is kindly acknowledged (A.F.). S.G. acknowledges financial support from the Villum Young Investigator grants 37440 and 13160 and the Cosmic Dawn Center. M.K. was supported by the ANID BASAL project FB210003. G.E.M. acknowledges financial support from the Villum Young Investigator grants 37440 and 13160 and the Cosmic Dawn Center. J.W. acknowledges support from the Science and Technology Facilities Council (STFC), by the ERC through Advanced Grant 695671 \'QUENCH\', by the UKRI Frontier Research grant RISEandFALL. Support for this work was provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51515.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. F.C. acknowledges support from a UKRI Frontier Research Guarantee Grant (PI Cullen; grant reference EP / X021025 / 1). J.H.W. acknowledges support by NSF grant AST-2108020 and NASA grants 80NSSC20K0520 and 80NSSC21K1053. NRT acknowledges support through STFC consolidated grant ST / W000857 / 1. M.J.H. is supported by the Swedish Research Council, Vetenskapsrådet, and is fellow of the Knut & Alice Wallenberg foundation.  \nThis work is based in part on observations made with the NASA / ESA / CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST .', 'Data availability statement': 'All the data presented in this work are made publicly available through DJA: https://dawn-cph.github.io/dja/ , and are made easy to inspect and query in a dedicated webpage: https://s3.amazonaws.com/msaexp-nirspec/ extractions/nirspec\\_graded\\_v2.html . The spectroscopic data have been processed using the custom-made reduction pipeline M sa E xp , publicly available here: https://github. com/gbrammer/msaexp . All the data and catalogs related to JWST -PRIMAL are also available on a dedicated webpage: https://github.com/keheintz/jwst-primal .  \nSoftware: This work made use of and acknowledge the following software: NumPy (Harris et al. 2020), Matplotlib (Hunter 2007), LMfit (Newville et al. 2014), SciPy (Virtanen et al. 2020), grizli (Brammer 2023a), Astrodrizzle (Gonzaga et al. 2012), and MsaExp (v0.3; Brammer 2023b).', 'References': "Adams, N. J., Conselice, C. J., Austin, D., et al. 2023, arXiv e-prints, arXiv:2304.13721  \nArellano-Córdova, K. Z., Berg, D. A., Chisholm, J., et al. 2022, ApJ, 940, L23 Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, Nature, 622, 707 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 Atek, H., Chemerynska, I., Wang, B., et al. 2023a, MNRAS, 524, 5486 Atek, H., Labbé, I., Furtak, L. J., et al. 2023b, arXiv e-prints, arXiv:2308.08540 Austin, D., Adams, N., Conselice, C. J., et al. 2023, ApJ, 952, L7 Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473 Begley, R., Cullen, F., McLure, R. J., et al. 2024, MNRAS, 527, 4040 Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488,  \n- 3143\n- Bezanson, R., Labbe, I., Whitaker, K. E., et al. 2022, arXiv e-prints, arXiv:2212.04026\n- Blumenthal, G. R., Faber, S. M., Primack, J. R., & Rees, M. J. 1984, Nature, 311, 517\n- Bosman, S. E. I., Davies, F. B., Becker, G. D., et al. 2022, MNRAS, 514, 55 Bouwens, R. J., Smit, R., Labbé, I., et al. 2016, ApJ, 831, 176\n- Bouwens, R. J., Stefanon, M., Brammer, G., et al. 2023, MNRAS, 523, 1036\n- Brammer, G. 2023a, grizli, Zenodo\n- Brammer, G. 2023b, msaexp: NIRSpec analyis tools, Zenodo\n- Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503\n- Brinchmann, J. 2023, MNRAS, 525, 2087\n- Bunker, A. J., Cameron, A. J., Curtis-Lake, E., et al. 2023a, arXiv e-prints, arXiv:2306.02467\n- Bunker, A. J., Saxena, A., Cameron, A. J., et al. 2023b, A&A, 677, A88 Calzetti, D. 2013, in Secular Evolution of Galaxies, ed. J. Falcón-Barroso & J. H. Knapen, 419\n- Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1994, ApJ, 429, 582 Cameron, A. G. W. 1957, PASP, 69, 201\n- Cameron, A. J., Katz, H., Witten, C., et al. 2023a, arXiv e-prints, arXiv:2311.02051\n- Cameron, A. J., Saxena, A., Bunker, A. J., et al. 2023b, A&A, 677, A115\n- Casey, C. M., Akins, H. B., Shuntov, M., et al. 2023, arXiv e-prints, arXiv:2308.10932\n- Castellano, M., Dayal, P., Pentericci, L., et al. 2016, ApJ, 818, L3\n- Castellano, M., Fontana, A., Treu, T., et al. 2022, ApJ, 938, L15\n- Castellano, M., Napolitano, L., Fontana, A., et al. 2024, arXiv e-prints, arXiv:2403.10238\n- Castellano, M., Pentericci, L., Vanzella, E., et al. 2018, ApJ, 863, L3 Cen, R. & Kimm, T. 2015, ApJ, 801, L25\n- Chemerynska, I., Atek, H., Furtak, L. J., et al. 2023, arXiv e-prints, arXiv:2312.05030\n- Chen, H. 2023, arXiv e-prints, arXiv:2307.04797\n- Chen, Z., Stark, D. P., Mason, C., et al. 2023, arXiv e-prints, arXiv:2311.13683\n- Chisholm, J., Rigby, J. R., Bayliss, M., et al. 2019, ApJ, 882, 182\n- Chisholm, J., Saldana-Lopez, A., Flury, S., et al. 2022, MNRAS, 517, 5104\n- Choustikov, N., Katz, H., Saxena, A., et al. 2024, arXiv e-prints, arXiv:2401.09557\n- Cooke, J. & O'Meara, J. M. 2015, ApJ, 812, L27\n- Cullen, F., McLeod, D. J., McLure, R. J., et al. 2023a, arXiv e-prints, arXiv:2311.06209\n- Cullen, F., McLure, R. J., Khochfar, S., Dunlop, J. S., & Dalla Vecchia, C. 2017, MNRAS, 470, 3006\n- Cullen, F., McLure, R. J., McLeod, D. J., et al. 2023b, MNRAS, 520, 14\n- Cullen, F., Shapley, A. E., McLure, R. J., et al. 2021, MNRAS, 505, 903\n- Curti, M., D'Eugenio, F., Carniani, S., et al. 2023a, MNRAS, 518, 425\n- Curti, M., Maiolino, R., Curtis-Lake, E., et al. 2023b, arXiv e-prints, arXiv:2304.08516\n- Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2020, MNRAS, 491, 944 Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, Nature Astronomy, 7, 622\n- Davies, F. B., Hennawi, J. F., Bañados, E., et al. 2018, ApJ, 864, 142 Dayal, P. & Ferrara, A. 2018, Phys. Rep., 780, 1\n- Dayal, P., Maselli, A., & Ferrara, A. 2011, MNRAS, 410, 830\n- Dayal, P., Volonteri, M., Greene, J. E., et al. 2024, arXiv e-prints, arXiv:2401.11242\n- D'Eugenio, F., Maiolino, R., Carniani, S., et al. 2023, arXiv e-prints, arXiv:2311.09908  \nDewdney, P. E., Hall, P. J., Schilizzi, R. T., & Lazio, T. J. L. W. 2009, IEEE Proceedings, 97, 1482  \nDijkstra, M., Gronke, M., & Venkatesan, A. 2016, ApJ, 828, 71 D'Odorico, V., Feruglio, C., Ferrara, A., et al. 2018, ApJ, 863, L29 Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011 Duan, Q., Conselice, C. J., Li, Q., et al. 2023, arXiv e-prints, arXiv:2309.14961 Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023, arXiv e-prints, arXiv:2306.02465  \nEndsley, R., Stark, D. P., Chevallard, J., & Charlot, S. 2021, MNRAS, 500, 5229  \nEndsley, R., Stark, D. P., Whitler, L., et al. 2023, MNRAS, 524, 2312  \nFan, X., Bañados, E., & Simcoe, R. A. 2023, ARA&A, 61, 373  \nFerrara, A., Pallottini, A., & Dayal, P. 2023, MNRAS, 522, 3986  \nFinkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023a, ApJ, 946, L13  \nFinkelstein, S. L., D'Aloisio, A., Paardekooper, J.-P., et al. 2019, ApJ, 879, 36  \nFinkelstein, S. L., Leung, G. C. K., Bagley, M. B., et al. 2023b, arXiv e-prints,  \narXiv:2311.04279  \nFinlator, K., Keating, L., Oppenheimer, B. D., Davé, R., & Zackrisson, E. 2018, MNRAS, 480, 2628  \nFlury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022, ApJ, 930, 126  \nFranco,  \nM.,  \nAkins,  \nH. B., Casey, C. M., et al. 2023, arXiv e-prints,  \narXiv:2308.00751  \nFujimoto, S., Arrabal Haro, P., Dickinson, M., et al. 2023a, ApJ, 949, L25 Fujimoto, S., Wang, B., Weaver, J., et al. 2023b, arXiv e-prints, arXiv:2308.11609  \nFynbo, J. P. U., Jakobsson, P., Prochaska, J. X., et al. 2009, ApJS, 185, 526 Gonzaga, S., Hack, W., Fruchter, A., & Mack, J. 2012, The DrizzlePac Handbook Greig, B., Mesinger, A., & Bañados, E. 2019, MNRAS, 484, 5094  \nGreig, B., Mesinger, A., Haiman, Z., & Simcoe, R. A. 2017, MNRAS, 466, 4239  \nGronke, M., Dijkstra, M., McCourt, M., & Oh, S. P. 2017, A&A, 607, A71  \nHainline, K. N., Johnson, B. D., Robertson, B., et al. 2023, arXiv e-prints,  \narXiv:2306.02468  \nHarikane, Y., Ouchi, M., Oguri, M., et al. 2023, ApJS, 265, 5  \nHarris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357  \nHashimoto, T., Laporte, N., Mawatari, K., et al. 2018, Nature, 557, 392  \nHaslbauer, M., Kroupa, P., Zonoozi, A. H., & Haghi, H. 2022, ApJ, 939, L31  \nHayes, M. J., Runnholm, A., Scarlata, C., Gronke, M., & Rivera-Thorsen, T. E.  \n2023, MNRAS, 520, 5903  \nHayes, M. J. & Scarlata, C. 2023, ApJ, 954, L14  \n- Heintz, K. E., Brammer, G. B., Giménez-Arteaga, C., et al. 2023a, Nature Astronomy\n- Heintz, K. E., De Cia, A., Thöne, C. C., et al. 2023b, arXiv e-prints, arXiv:2308.14812\n- Heintz, K. E., Giménez-Arteaga, C., Fujimoto, S., et al. 2023c, ApJ, 944, L30 Heintz, K. E., Oesch, P. A., Aravena, M., et al. 2022, ApJ, 934, L27\n- Heintz, K. E., Watson, D., Brammer, G., et al. 2023d, arXiv e-prints, arXiv:2306.00647\n- Heintz, K. E., Watson, D., Oesch, P. A., Narayanan, D., & Madden, S. C. 2021, ApJ, 922, 147\n- Henry, A., Scarlata, C., Martin, C. L., & Erb, D. 2015, ApJ, 809, 19  \nHorne, K. 1986, PASP, 98, 609  \nHoyle, F. 1954, ApJS, 1, 121  \n- Hsiao, T. Y.-Y., Abdurro'uf, Coe, D., et al. 2023, arXiv e-prints, arXiv:2305.03042  \nHu, W., Martin, C. L., Gronke, M., et al. 2023, ApJ, 956, 39  \nHunter, J. D. 2007, Computing in Science and Engineering, 9, 90  \nHutter, A., Trebitsch, M., Dayal, P., et al. 2023, MNRAS, 524, 6124  \n- Inayoshi, K., Harikane, Y., Inoue, A. K., Li, W., & Ho, L. C. 2022, ApJ, 938, L10\n- Izotov, Y. I., Stasi'nska, G., Meynet, G., Guseva, N. G., & Thuan, T. X. 2006, A&A, 448, 955  \nIzotov, Y. I., Worseck, G., Schaerer, D., et al. 2018, MNRAS, 478, 4851 Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80 Jakobsson, P., Fynbo, J. P. U., Ledoux, C., et al. 2006, A&A, 460, L13 Jaskot, A. E., Dowd, T., Oey, M. S., Scarlata, C., & McKinney, J. 2019, ApJ, 885, 96  \nJiang, H., Wang, X., Cheng, C., et al. 2023, arXiv e-prints, arXiv:2312.04151  \nJones,  \nG. C., Bunker, A. J., Saxena, A., et al.  \n2023,  \narXiv  \ne-prints,  \narXiv:2306.02471  \n- Jung, I., Finkelstein, S. L., Arrabal Haro, P., et al. 2023, arXiv e-prints, arXiv:2304.05385  \nKatz, H., Saxena, A., Cameron, A. J., et al. 2023, MNRAS, 518, 592  \n- Keating, L. C., Bolton, J. S., Cullen, F., et al. 2023a, arXiv e-prints, arXiv:2308.05800  \nKeating, L. C., Puchwein, E., Bolton, J. S., Haehnelt, M. G., & Kulkarni, G. 2023b, arXiv e-prints, arXiv:2308.11709  \nKeating, L. C., Weinberger, L. H., Kulkarni, G., et al. 2020, MNRAS, 491, 1736 Kennicutt, Robert C., J. 1998, ARA&A, 36, 189  \nKennicutt, R. C. & Evans, N. J. 2012, ARA&A, 50, 531  \nKereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2  \nKewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511  \nKhostovan, A. A., Sobral, D., Mobasher, B., et al. 2016, MNRAS, 463, 2363 Kroupa, P. 2001, MNRAS, 322, 231  \n- Labbé, I., van Dokkum, P., Nelson, E., et al. 2023, Nature, 616, 266\n- Langeroodi, D. & Hjorth, J. 2023, arXiv e-prints, arXiv:2307.06336\n- Langeroodi, D., Hjorth, J., Chen, W., et al. 2023, ApJ, 957, 39  \nLara-López, M. A., Cepa, J., Bongiovanni, A., et al. 2010, A&A, 521, L53  \nLaursen, P., Sommer-Larsen, J., Milvang-Jensen, B., Fynbo, J. P. U., & Ra-  \nzoumov, A. O. 2019, A&A, 627, A84 Laursen, P., Sommer-Larsen, J., & Razoumov, A. O. 2011, ApJ, 728, 52 Liang, L., Feldmann, R., Murray, N., et al. 2024, MNRAS, 528, 499 Lin, X., Cai, Z., Zou, S., et al. 2023, ApJ, 944, L59 Lin, Y.-H., Scarlata, C., Williams, H., et al. 2024, MNRAS, 527, 4173 Madau, P. & Dickinson, M. 2014, ARA&A, 52, 415 Maiolino, R. & Mannucci, F. 2019, A&A Rev., 27, 3 Malkan, M. A., Cohen, D. P., Maruyama, M., et al. 2017, ApJ, 850, 5  \n- Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115  \nMason, C. A., Fontana, A., Treu, T., et al. 2019, MNRAS, 485, 3947 Mason, C. A., Trenti, M., & Treu, T. 2023, MNRAS, 521, 497 Mason, C. A., Treu, T., Dijkstra, M., et al. 2018, ApJ, 856, 2 Matthee, J., Mackenzie, R., Simcoe, R. A., et al. 2023, ApJ, 950, 67 Matthee, J., Sobral, D., Gronke, M., et al. 2018, A&A, 619, A136 McLeod, D. J., Donnan, C. T., McLure, R. J., et al. 2024, MNRAS, 527, 5004 McQuinn, M., Hernquist, L., Zaldarriaga, M., & Dutta, S. 2007, MNRAS, 381,  \n- 75  \nMcQuinn, M., Lidz, A., Zaldarriaga, M., Hernquist, L., & Dutta, S. 2008, MNRAS, 388, 1101  \n- Mesinger, A. & Haiman, Z. 2007, ApJ, 660, 923\n- Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64  \nMiralda-Escudé, J. 1998, ApJ, 501, 15  \n- Mirocha, J. & Furlanetto, S. R. 2023, MNRAS, 519, 843\n- Naidu, R. P., Oesch, P. A., Setton, D. J., et al. 2022, arXiv e-prints, arXiv:2208.02794  \nNaidu, R. P., Tacchella, S., Mason, C. A., et al. 2020, ApJ, 892, 109 Nakajima, K., Ouchi, M., Isobe, Y., et al. 2023, arXiv e-prints, arXiv:2301.12825 Nakane, M., Ouchi, M., Nakajima, K., et al. 2023, arXiv e-prints, arXiv:2312.06804  \n- Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT: NonLinear Least-Square Minimization and Curve-Fitting for Python  \nNoterdaeme, P., Petitjean, P., Carithers, W. C., et al. 2012, A&A, 547, L1  \nOcvirk, P., Lewis, J. S. W., Gillet, N., et al. 2021, MNRAS, 507, 6108  \n- Osterbrock, D. E. & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei\n- Ouchi, M., Shimasaku, K., Furusawa, H., et al. 2010, ApJ, 723, 869 Pacucci, F., Dayal, P., Harikane, Y., Inoue, A. K., & Loeb, A. 2022, MNRAS, 514, L6\n- Pahl, A. J., Shapley, A. E., Steidel, C. C., et al. 2024, arXiv e-prints, arXiv:2401.09526\n- Pallottini, A., Ferrara, A., Gallerani, S., et al. 2022, MNRAS, 513, 5621  \nPéroux, C. & Howk, J. C. 2020, ARA&A, 58, 363  \n- Pettini, M., Steidel, C. C., Adelberger, K. L., Dickinson, M., & Giavalisco, M. 2000, ApJ, 528, 96  \nPlanck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6 Prochaska, J. X., Chen, H.-W., Dessauges-Zavadsky, M., & Bloom, J. S. 2007, ApJ, 666, 267  \nProchaska, J. X. & Wolfe, A. M. 2009, ApJ, 696, 1543  \nPuchwein, E., Bolton, J. S., Keating, L. C., et al. 2023, MNRAS, 519, 6162  \nRasmussen Cueto, E., Hutter, A., Dayal, P., et al. 2023, arXiv e-prints,  \narXiv:2312.12109  \nReddy, N. A., Oesch, P. A., Bouwens, R. J., et al. 2018a, ApJ, 853, 56 Reddy, N. A., Shapley, A. E., Sanders, R. L., et al. 2018b, ApJ, 869, 92 Rhoads, J. E., Wold, I. G. B., Harish, S., et al. 2023, ApJ, 942, L14 Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001 Roberts-Borsani, G., Treu, T., Chen, W., et al. 2023, Nature, 618, 480 Roberts-Borsani, G., Treu, T., Shapley, A., et al. 2024, arXiv e-prints, arXiv:2403.07103  \nRobertson, B. E. 2022, ARA&A, 60, 121  \nRobertson, B. E., Furlanetto, S. R., Schneider, E., et al. 2013, ApJ, 768, 71  \nRosdahl, J., Blaizot, J., Katz, H., et al. 2022, MNRAS, 515, 2386  \n- Salmon, B., Papovich, C., Finkelstein, S. L., et al. 2015, ApJ, 799, 183  \nSanders, R. L., Shapley, A. E., Jones, T., et al. 2021, ApJ, 914, 19  \n- Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2024, ApJ, 962, 24\n- Sandles, L., Curtis-Lake, E., Charlot, S., Chevallard, J., & Maiolino, R. 2022, MNRAS, 515, 2951\n- Saxena, A., Bunker, A. J., Jones, G. C., et al. 2023a, arXiv e-prints, arXiv:2306.04536\n- Saxena, A., Robertson, B. E., Bunker, A. J., et al. 2023b, A&A, 678, A68 Schaerer, D. 2003, A&A, 397, 527\n- Schaerer, D., Marques-Chaves, R., Barrufet, L., et al. 2022, A&A, 665, L4  \nSchaye, J., Dalla Vecchia, C., Booth, C. M., et al. 2010, MNRAS, 402, 1536 Serjeant, S. & Bakx, T. J. L. C. 2023, Nature Astronomy, 7, 1143  \n- Shapley, A. E., Sanders, R. L., Reddy, N. A., Topping, M. W., & Brammer, G. B. 2023, ApJ, 954, 157\n- Shapley, A. E., Steidel, C. C., Pettini, M., & Adelberger, K. L. 2003, ApJ, 588, 65\n- Shen, X., Vogelsberger, M., Boylan-Kolchin, M., Tacchella, S., & Kannan, R. 2023, MNRAS, 525, 3254  \nSimmonds, C., Tacchella, S., Maseda, M., et al. 2023, MNRAS, 523, 5468 Sobacchi, E. & Mesinger, A. 2015, MNRAS, 453, 1843  \n- Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15\n- Stark, D. P. 2016, ARA&A, 54, 761\n- Steidel, C. C., Bogosavljevi'c, M., Shapley, A. E., et al. 2018, ApJ, 869, 123\n- Steidel, C. C., Erb, D. K., Shapley, A. E., et al. 2010, ApJ, 717, 289\n- Steidel, C. C., Strom, A. L., Pettini, M., et al. 2016, ApJ, 826, 159\n- Steinhardt, C. L., Kokorev, V., Rusakov, V., Garcia, E., & Sneppen, A. 2023, ApJ, 951, L40\n- Sun, G., Faucher-Giguère, C.-A., Hayward, C. C., et al. 2023, ApJ, 955, L35 Tang, M., Stark, D. P., Chen, Z., et al. 2023, MNRAS, 526, 1657\n- Tanvir, N. R., Fynbo, J. P. U., de Ugarte Postigo, A., et al. 2019, MNRAS, 483, 5380  \nTaylor, A. J., Barger, A. J., & Cowie, L. L. 2022, ApJ, 939, L3 Tepper-García, T. 2006, MNRAS, 369, 2025  \n- Thorne, J. E., Robotham, A. S. G., Davies, L. J. M., et al. 2021, MNRAS, 505, 540\n- Topping, M. W., Stark, D. P., Endsley, R., et al. 2023, arXiv e-prints, arXiv:2307.08835\n- Torrey, P., Vogelsberger, M., Hernquist, L., et al. 2018, MNRAS, 477, L16 Totani, T., Kawai, N., Kosugi, G., et al. 2006, PASJ, 58, 485\n- Treu, T., Roberts-Borsani, G., Bradac, M., et al. 2022, ApJ, 935, 110\n- Trinca, A., Schneider, R., Valiante, R., et al. 2023, arXiv e-prints, arXiv:2305.04944\n- Trump, J. R., Arrabal Haro, P., Simons, R. C., et al. 2023, ApJ, 945, 35 Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389 Umeda, H., Ouchi, M., Nakajima, K., et al. 2023, arXiv e-prints, arXiv:2306.00487  \nValentino, F., Brammer, G., Gould, K. M. L., et al. 2023, ApJ, 947, 20 van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28 ˇ Durovˇcíková, D., Eilers, A.-C., Chen, H., et al. 2024, arXiv e-prints, arXiv:2401.10328  \nˇ Durovˇcíková, D., Katz, H., Bosman, S. E. I., et al. 2020, MNRAS, 493, 4256 Verhamme, A., Orlitová, I., Schaerer, D., & Hayes, M. 2015, A&A, 578, A7 Verhamme, A., Orlitová, I., Schaerer, D., et al. 2017, A&A, 597, A13 Vijayan, A. P., Thomas, P. A., Lovell, C. C., et al. 2024, MNRAS, 527, 7337 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261 Vizgan, D., Heintz, K. E., Greve, T. R., et al. 2022, ApJ, 939, L1 Wang, B., Fujimoto, S., Labbe, I., et al. 2023a, arXiv e-prints, arXiv:2308.03745 Wang, B., Leja, J., Labbé, I., et al. 2023b, arXiv e-prints, arXiv:2310.01276 Wang, F., Yang, J., Fan, X., et al. 2021, ApJ, 907, L1  \n- Ward, E. M., de la Vega, A., Mobasher, B., et al. 2023, arXiv e-prints, arXiv:2311.02162\n- Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, L29  \nWhite, S. D. M. & Frenk, C. S. 1991, ApJ, 379, 52  \nWhite, S. D. M. & Rees, M. J. 1978, MNRAS, 183, 341  \n- Wilson, S. N., Heintz, K. E., Jakobsson, P., et al. 2023, arXiv e-prints, arXiv:2305.05213  \nWitstok, J., Smit, R., Saxena, A., et al. 2023, arXiv e-prints, arXiv:2306.04627 Witten, C., Laporte, N., Martin-Alvarez, S., et al. 2024, Nature Astronomy Wolfe, A. M., Gawiser, E., & Prochaska, J. X. 2005, ARA&A, 43, 861 Wolfe, A. M., Turnshek, D. A., Smith, H. E., & Cohen, R. D. 1986, ApJS, 61, 249  \nWoosley, S. E. & Weaver, T. A. 1995, ApJS, 101, 181 Xu, X., Henry, A., Heckman, T., et al. 2022, ApJ, 933, 202 Yang, H., Malhotra, S., Gronke, M., et al. 2017, ApJ, 844, 171 Yang, J., Wang, F., Fan, X., et al. 2020, ApJ, 904, 26 Yeh, J. Y. C., Smith, A., Kannan, R., et al. 2023, MNRAS, 520, 2757 Zackrisson, E., Rydberg, C.-E., Schaerer, D., Östlin, G., & Tuli, M. 2011, ApJ, 740, 13  \n- Zavala, J. A., Castellano, M., Akins, H. B., et al. 2024, arXiv e-prints, arXiv:2403.10491\n- 1 Cosmic Dawn Center (DAWN), Denmark\n- 2 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen N, Denmark\n- 3 Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland\n- 4 Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK\n- 5 Stockholm University, Department of Astronomy and Oskar Klein Centre for Cosmoparticle Physics, AlbaNova University Centre, SE10691, Stockholm, Sweden\n- 6 Department of Physics and Astronomy, The Johns Hopkins University, 3400 N Charles St. Baltimore, MD 21218, USA\n- 7 Space Telescope Science Institute (STScI), 3700 San Martin Drive, Baltimore, MD 21218, USA\n- 8 Kapteyn Astronomical Institute, University of Groningen, 9700 AV Groningen, The Netherlands\n- 9 Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa, Italy\n- 10 New Mexico State University, Las Cruces, 88003 NM, USA\n- 11 Department of Astronomy, University of Massachusetts Amherst, Amherst, MA 01002, United States\n- 12 DTU Space, Technical University of Denmark, Elektrovej, Building 328, 2800, Kgs. Lyngby, Denmark\n- 13 Center for Relativistic Astrophysics, School of Physics, Georgia Institute of Technology, 837 State Street, Atlanta, GA 30332, USA\n- 14 Centre for Astrophysics and Cosmology, Science Institute, University of Iceland, Dunhagi 5, 107, Reykjavik, Iceland\n- 15 Instituto de Estudios Astrofísicos, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército 441, Santiago 8370191, Chile\n- 16 Department of Physics and Astronomy, University of California, Riverside, 900 University Avenue, Riverside, CA 92521, USA\n- 17 Institute of Science and Technology Austria (ISTA), Am Campus 1, 3400 Klosterneuburg, Austria\n- 18 MIT Kavli Institute for Astrophysics and Space Research, 77 Massachusetts Ave., Cambridge, MA 02139, USA\n- 19 Department of Astronomy, University of Florida, 211 Bryant Space Sciences Center, Gainesville, FL 32611, USA\n- 20 Astrophysics Research Institute, Liverpool John Moores University, Liverpool, L35 UG, UK\n- 21 School of Physics and Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK\n- 22 European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching, Germany\n- 23 Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK\n- 24 Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge CB3 0HE, UK  \nArticle number, page 20 of 20"}