Primordial neutral atomic gas mostly composed of hydrogen is the raw material for star formation in galaxies. However there are few direct constraints on the amount of neutral atomic hydrogen H I in galaxies at early cosmic times. We analyzed James Webb Space Telescope JWST nearinfrared spectroscopy of distant galaxies at redshifts 8. From a sample of 12 galaxies we identified three that show strong damped Lyman absorption due to H I in their local surroundings. The galaxies are located at spectroscopic redshifts of 8.8 10.2 and 11.4 corresponding to 400 to 600 million years after the Big Bang. They have H I column densities 10SUP22SUP cmSUP2SUP which is an order of magnitude higher than expected for a fully neutral intergalactic medium and constitute a gasrich population of young starforming galaxies.
['Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
Strong damped Lyman absorption in young starforming galaxies at redshifts 9 to 11
{'Extreme damped Lymanα absorption in young star-forming galaxies at z = 9 -11': "1 , 2 † 1 , 2 1 , 2  \nKasper E. Heintz, Darach Watson, Gabriel Brammer, Simone Vejlgaard, 1 , 2 Anne Hutter, 1 , 2 Victoria B. Strait, 1 , 2 Jorryt Matthee, 3 Pascal A. Oesch, 4 , 1 , 2 P'all Jakobsson, 5 Nial R. Tanvir, 6 Peter Laursen, 1 , 2 Rohan P. Naidu, 7 , ‡ Charlotte A. Mason, 1 , 2 Meghana Killi, 1 , 2 Intae Jung, 8 Tiger Yu-Yang Hsiao, 9 Abdurro'uf, 9 , 8 Dan Coe, 8 , 9 , 10 Pablo Arrabal Haro, 11 Steven L. Finkelstein, 12 &Sune Toft 1 , 2  \n1 Cosmic Dawn Center (DAWN), Denmark  \n2 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N, Denmark 3 Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 27, Zurich, 8093, Switzerland 4 Observatoire de Gen'eve, Universit'e de Gen'eve, Chemin Pegasi 51, CH-1290 Versoix, Switzerland 5 Centre for Astrophysics and Cosmology, Science Institute, University of Iceland, Dunhagi 5, 107 Reykjav'ık, Iceland 6 School of Physics and Astronomy, University of Leicester, University Road, Leicester, LE1 7RH, UK 7 MIT Kavli Institute for Astrophysics and Space Research, 77 Massachusetts Avenue, Cambridge, 02139, Massachusetts, USA 8 Space Telescope Science Institute (STScI), 3700 San Martin Drive, Baltimore, MD 21218, USA 9 Center for Astrophysical Sciences, Department of Physics and Astronomy, The Johns Hopkins University, 3400 N Charles St. Baltimore, MD 21218, USA 10 Association of Universities for Research in Astronomy (AURA), Inc. for the European Space Agency (ESA) 11 NSF's National Optical-Infrared Astronomy Research Laboratory, 950 N. Cherry Ave., Tucson, AZ 85719, USA 12 Department of Astronomy, The University of Texas at Austin, Austin, TX, USA † Corresponding author; E-mail: keheintz@nbi.ku.dk ‡ NASA Hubble Fellow  \nThe onset of galaxy formation is thought to be initiated by the infall of neutral, pristine gas onto the first protogalactic halos. However, direct constraints  \non the abundance of neutral atomic hydrogen (H I ) in galaxies have been difficult to obtain at early cosmic times. Here we present spectroscopic observations with JWST of three galaxies at redshifts z = 8 . 8 -11 . 4 , about 400600 Myr after the Big Bang, that show strong damped Lymanα absorption ( N H I > 10 22 cm -2 ) from H I in their local surroundings, an order of magnitude in excess of the Lymanα absorption caused by the neutral intergalactic medium at these redshifts. Consequently, these early galaxies cannot be contributing significantly to reionization, at least at their current evolutionary stages. Simulations of galaxy formation show that such massive gas reservoirs surrounding young galaxies so early in the history of the universe is a signature of galaxy formation in progress.  \nIn the current paradigm of galaxy formation, primordial neutral atomic hydrogen (H I) accretes onto galaxy halos, cools, and forms stars ( 1, 2 ). Detecting the accretion and build-up of HI gas in the earliest galaxies is therefore central to understanding the first phases of galaxy formation. Constraining the amount of H I in galaxies at z ≳ 8 is also important to quantify how much ionizing radiation escapes from normal star-forming galaxies at this epoch as they are likely to initiate and serve as the main driver of the large-scale reionization ( 3 ). Unfortunately, direct detection of H I in emission from galaxies at cosmological distances is severely limited by the weakness of the H I 21 cm hyperfine transition ( 4, 5 ).  \nAlternatively, H I can be measured at high redshifts through the Lymanα (Ly α ) absorption feature imposed on the spectra of bright background point sources such as quasars. Ly α absorbers with column densities N H I ≥ 2 × 10 20 cm -2 , known as damped Ly α absorbers (DLAs), reveal the presence of extended, dense gas reservoirs surrounding star-forming galaxies ( 6 ). Similar, but far more centrally-located DLAs are also observed within the host galaxies of γ -ray bursts. However, extended DLAs covering the UV emission of entire galaxies are only  \nobserved in the most intensely ionized, metal-poor local galaxies at z ≈ 0 ( 7 ) or in Lymanbreak galaxies at z ≈ 3 with similar extreme physical properties ( 8 ), and always with column densities N H I ≤ 5 × 10 21 cm -2 . With the advent of JWST, we are now able to uncover and characterize in unprecedented detail galaxies during the critical period in the Universe's history, at redshifts ≳ 8 , when the Universe was being reionized and the first galaxies developed. In this paper, we present observations of three early universe galaxies, at redshifts z ≃ 9 -11 , which show extreme DLAs with N H I > 10 22 cm -2 . We argue that the discovery of this feature in several galaxies of the small sample that has so far been observed suggests that we are now beginning to observe the early development phase for galaxies: the build-up of pristine, neutral gas prior to the production of the bulk of the stars.  \nWe analyse the JWST observations of three select galaxies: two were observed and identified as high-redshift galaxy candidates in JWST/NIRCam imaging obtained through the Cosmic Evolution Early Release Science (CEERS, DD-ERS 1345, PI: S. Finkelstein) Survey ( 9-11 ) and later spectroscopically confirmed with dedicated JWST/NIRSpec prism spectroscopic followup observations (Prog. ID: DD-2750, PI: P. Arrabal Haro), CEERS-43833 with MSA ID 28 and CEERS-16943 with MSA ID 1 ( 12 ). The third source, MACS0647-JD, is gravitationally-lensed by the foreground galaxy cluster MACS J0647+7015 and was discovered as part of the Cluster Lensing And Supernova survey with Hubble (CLASH) ( 13 ). JWST/NIRCam and NIRSpec prism observations (observation program GO 1433, PI: D. Coe) show a system with multiple intrinsic components (A and B) which are then multiply-imaged ( 14 ), and the galaxy spectrum covering mainly component A shows a star-forming system at z = 10 . 17 based on the identification of seven emission lines ( 15 ).  \nThe spectroscopic and imaging data were obtained from the Mikulski Archive for Space Telescopes (MAST) and reduced using grizli and a custom-made pipeline (see Supplementary Materials). We scale and match the optimally-extracted 1D spectra to the available  \nJWST/NIRCam photometry using a wavelength-dependent polynomial function. This is to improve the absolute flux calibration and take into account potential slit-losses. The reduced and photometrically-calibrated 1D spectra of the three galaxies at z = 9 -11 are shown in Fig. 1. We determine the redshift of CEERS-43833 to be z = 8 . 7622 ± 0 . 0002 from the multitude of nebular emission lines detected in the spectrum (see Supplementary Materials, Fig. S2), consistent with the first analysis of this spectrum presented in ( 12 ). Similarly, we determine z = 11 . 409 ± 0 . 001 for CEERS-16943 based on the detection of [O II] λ 3727 , [Ne III] λ 3870 , and He I λ 3889 . This is the most distant detection of nebular emission lines to date. For MACS0647-JD we derive z = 10 . 170 ± 0 . 003 , consistent with the previous spectroscopic analysis of this galaxy ( 15 ). To select these galaxies, we began by examining the publicly available JWST spectra for galaxies at z > 8 , looking for any exhibiting a slow rollover at the wavelength region around the Ly α transition (Lyman-break), indicative of a strong damping wing absorption. Such signatures appear to be common at these redshifts, although not universal (e.g. ( 16, 17 ), also see Fig. S4). From about a dozen candidates, the galaxies selected here offer the three most robust examples, with clear emission lines pinning down the systemic redshift and high signal-to-noise (S/N) around Ly α . We therefore investigate these three galaxies as the best exemplars of this particular feature in more detail.  \nTo model the intrinsic continuum galaxy emission, we used the recent set of galaxy templates from ( 18 ), created from BPASS and CLOUDY models that are specifically designed to reproduce the blue rest-frame ultraviolet (UV) colors of galaxies at z > 8 . In order to account for absorption from the intergalactic medium (IGM), we model the optical depth of Ly α due to the Gunn-Peterson effect from an increasingly neutral IGM ( 19 ) as  \nτ GP ( z ) = 1 . 8 × 10 5 h -1 (Ω DM , 0 ) 1 / 2 x H I ( Ω m , 0 h 2 0 . 02 )( 1 + z 7 ) 3 / 2 , (1)  \nwhere h = H 0 / (100 km s -1 Mpc -1 ) is the present-day dimensionless Hubble parameter, Ω m , 0  \nand Ω DM , 0 the current baryonic and dark matter density parameters, and x H I is the average neutral to total hydrogen fraction of the IGM. This equation only holds outside the damping wing of the absorption line, while the damping wing itself will be affected by the velocity distribution of the H I atoms. This correction is given by ( 20 ) as  \nτ IGM ( λ obs ) = x H I R α τ GP ( z gal ) π ( 1 + z obs 1 + z gal ) 3 / 2 × [ I ( 1 + z IGM , u 1 + z obs ) -I ( 1 + z IGM , l 1 + z obs )] (2)  \nfollowing the formalism of ( 21 ). Here, z obs is the observed redshift of the neutral gas, z gal is the galaxy redshift, R α = Λ α λ α / (4 πc ) = 2 . 02 × 10 -8 is a dimensionless quantity that includes the damping constant of the Ly α resonance Λ α , and the Ly α wavelength λ α . This parameterization assumes that the IGM is uniformly distributed within the redshift range z IGM , l to z IGM , u . Following ( 21 ), we set the upper bound at the galaxy redshift z IGM , u = z gal and integrate the expression down to z IGM , l = 6 . The resulting spectral cut-off is not sensitive to the choice of z IGM , l at redshifts much below z gal ( 21 ). We further add a component describing the optical depth due to Ly α absorption from the interstellar medium (ISM), or gas in the immediate surroundings of the galaxy, in the shape of a Voigt profile  \nτ ISM = CaH ( a, x ) N H I (3)  \nwhere C is the photon absorption constant, a is the damping parameter, and H ( a, x ) is the VoigtHjerting function, following the approximation of ( 22 ). In the case of strong DLA features, where the Lorentzian wings dominate the shape of the line profile, the optical depth is only sensitive to the column density of neutral atomic hydrogen, N H I.  \nWe model the observed galaxy spectra by fixing the intrinsic galaxy template to the spectroscopic redshifts determined from the detected emission lines, and leave the dust extinction, A V , the neutral H I fraction of the IGM, x H I, and the column density of the Ly α absorber from the local galaxy environment, N H I, as free parameters. Specifically, we use the binc100z001age6 cloudy model following the recommendation of ( 18 ) but note that the  \nFigure 1: JWST/NIRSpec spectroscopic data. Main panels show the reduced and photometrically-calibrated NIRSpec/prism 1D spectra covering 1 µ m to 5 . 2 µ m (black) and the associated 1 σ error spectrum (grey). The galaxy ID and spectroscopic redshifts from the identified emission lines are stated in the top left. The MACS0647-JD spectrum has been demagnified assuming µ = 8 ± 1 ( 15 ). In the insets are shown the spectral regions covering the Ly α transition, with models overplotted showing varying neutral IGM H I fractions, x H I = [0 . 1 , 0 . 5 , 1 . 0] (blue) and the best-fit model with an additional DLA feature (red). An extremely strong DLA, log( N H I / cm 2 ) > 22 . 0 , is detected in all cases.  \n<!-- image -->  \nexact choice of the intrinsic template, considering for instance older stellar populations, does not have a significant impact on the result ( ≲ 0 . 1 dex on N H I). We sample the posterior distributions using the PyMultiNest -package ( 23 ), which uses a Python wrapper around the MultiNest algorithm ( 24 ) to run a multimodal Nested Sampling algorithm on the spectrum. We derive bestfit models with log( N H I / cm 2 ) = 22 . 10 ± 0 . 20 for CEERS-43833, log( N H I / cm 2 ) = 22 . 20 ± 0 . 40 for CEERS-16943, and log( N H I / cm 2 ) = 22 . 40 ± 0 . 15 for MACS0647-JD. The uncertainties represent the 16th and 84th percentiles of the posterior distributions and includes the systematic uncertainty of 0.1 dex due to variations in the intrinsic galaxy continuum emission templates. The models are largely insensitive to the exact value of x H I, since the line profiles of the DLAs dominate the optical depth at λ α . To highlight this, we show various models with x H I = 0 . 1 , 0 . 5 , 1 . 0 (i.e. 10%, 50%, and 100% neutral IGM) in Fig. 1. It is evident that none of these pure IGM models are able to capture the broad absorption features observed in the JWST spectra. While such average sightline models are clearly ruled out even in the case of a fully neutral IGM, we also consider a scenario where these extreme DLAs could be caused by the galaxies residing in overdensities of a largely neutral IGM. However, this seems unlikely since the overdensities required at this epoch would be a factor of 30 -50 larger than average, approaching galaxy ISM densities, a very rare occurrence. We do note, however, that there appears to be an overdensity in the CEERS field at z ∼ 8 . 7 ( 25, 26 ), at least × 5 over the field average for very bright, M UV < -20 galaxies, though still insufficient to explain the observed properties.  \nWe show in the Supplementary Materials, Fig. S6, the distributions derived from a cosmological hydrodynamical simulation using sightlines from galaxies with similar stellar masses and star-formation rates (SFRs) at z ≈ 8 . 8 to investigate the impact of density variability in the IGM. In those simulations ( 27 ), the 95% confidence limits are significantly offset from the observed damping wing, albeit assuming a low x H I fraction at this epoch. These models thus  \nprovide supporting evidence that the most likely interpretation of the DLAs is related to H I gas in the immediate surroundings of the galaxies. We also exclude scenarios involving nebular continuum emission and 2 γ processes as the cause of the damped Ly α profiles in the Supplementary Materials. We further compare the CEERS-43833 spectrum to a set of spectra of four galaxies at z ≈ 9 with similar quality data, i.e. high S/N around Ly α and with robust redshift measurements from nebular emission (Fig. S4). These comparison galaxies do not show similar broad absorption features, though are consistent with substantial IGM neutral fractions in the line of sight, substantiating the distinct detection of the extreme DLAs in the three z = 9 -11 galaxies in excess of Ly α absorption from the IGM only. In Fig. S4, we also compare the observed spectra to the intrinsic IGM and DLA models convolved with the JWST/NIRSpec prism wavelength-dependent spectral resolution, excluding the coarse spectral resolution as the cause of the observed broad Ly α line profiles.  \nTo examine the line emission and physical properties of the three z = 9 -11 galaxies in more detail, we measure the line fluxes of each transition detected in the spectra. For CEERS-43833, we detect at high significance the strong nebular emission lines [O II] λ 3727 , [O III] λλ 4960 , 5008 , He I λ 3889 and the Balmer lines H β , H γ , and H δ and recover the same emission lines features previously reported in MACS0647-JD ( 15 ). For CEERS-16943 we marginally detect [O II] λ 3727 , [Ne III] λ 3870 , and He I λ 3889 at ≈ 3 σ statistical significance. To constrain the line fluxes of each transition, we fit a set of superimposed Gaussian line profiles on the best-fit continuum galaxy model, tying the redshift and line widths of the emission lines, but keep the line fluxes for each transition as free parameters. We report the measured line fluxes in Table S1. We also detect absorption features from the Si IV λ 1397 and C IV 1550 transitions in CEERS-43833 at the same redshift as the nebular emission lines, further hinting at the large, abundant gas reservoir in the local environment of this galaxy. We derive SFRs in the range 1 -15 M ⊙ yr -1 based on the H β or [O II] λ 3727 line flux measurements and infer  \nTable 1: Physical properties of the three star-forming at z = 9 -11 with strong DLAs.  \na Measurements adopted from Hsiao et al. ( 15 ).  \nrelatively low metallicities of 12 + log(O / H) = 7 . 4 -7 . 7 , i.e. 5-10% solar abundance values, based on strong-line calibrations (Supplementary Materials). The three galaxies are all characterized by young stellar populations ( ≲ 100 Myr) and have stellar masses in the range M ⋆ = 10 8 -10 9 M ⊙ based on spectrophotometric fitting of their spectral energy distributions (SEDs), as summarized in Table 1. Overall, these three galaxies have physical properties consistent with the typical star-forming galaxy population at z = 7 -10 ( 17 ), see Fig. S3.  \nOur observations further enable a comparison between the properties of these three young, highz galaxies to green pea galaxies at redshifts z ≈ 0 , which are considered to be lower redshift analogs to highz galaxies ( 7 ). In Fig. 2, we show the H I column density of these extreme z ≈ 0 emission-line galaxies versus their gas-phase metallicity and ionization parameter O 32 , including also the most metal-poor local galaxies known: I Zw 18, SBS 0335 -052 , and SBS 1415 + 437 ( 28 ). The three galaxies at z = 9 -11 show substantially higher H I gas columns than any of these local extreme galaxies, even at similar low gas-phase metallicities. Wealso observe a hint of increasing H I column densities for galaxies with higher ionization parameters. This is somewhat at odds with predictions and previous empirical results ( 29,30 ) that  \n<!-- image -->  \nFigure 2: H I column density relations. ( Left: ) The H I abundances of the three young, z = 9 -11 star-forming galaxies (red stars) compared to local green pea (GP) galaxies (blue dots) ( 7 ) and the three most metal-poor galaxies at z ≈ 0 (grey squares) ( 28 ) as a function of gas-phase metallicity, 12 + log(O / H) . ( Right: ) H I abundances as a function of the O 32 ionization. The symbols denote the same galaxy samples. The z = 9 -11 young star-forming galaxies show the most abundant DLAs.  \n<!-- image -->  \nshow that larger Ly α emission equivalent widths and escape fractions are expected from galaxies with higher O 32 ratios. The complete absorption of Ly α in these three highz galaxies might thus be additional signatures of their young stellar populations, which have not yet had time to significantly ionize even their local surroundings, as a prelude to initiating large-scale IGM reionization. Alternatively, the strong DLA features in these particular galaxies could be related to viewing angle effects, indicating that we observe particular gas-rich sightlines, whereas other low column density channels (as indicated by the high O 32 ) might exist. Such strongly directionally-dependent escape of ionizing radiation has been proposed by simulations ( 31-33 ). However, due to the complete absorption of Ly α , the covering fraction must be larger than the entire rest-frame UV emission of the galaxies in the line of sight.  \nBased on the large H I column densities and the complete covering fraction of Ly α , we can place lower bounds on the H I gas masses of the three z = 9 -11 galaxies. For CEERS43833 and CEERS-16943, we derive effective areas of the projected ellipses from the measured  \nFWHM major and minor axes of the galaxies, yielding A e = 0 . 67 kpc 2 and A e = 0 . 28 kpc 2 , respectively. For MACS0647-JD, we take the distance from the edges of the A and B components to be the effective diameter of the galaxy, ∼ 0 . 5 kpc ( 14 ), resulting in an area of A e = 0 . 25 kpc 2 . This results in columnar H I gas masses, M H I = m H I N H I A e , of approx. 6 . 7 × 10 7 M ⊙ , 3 . 5 × 10 7 M ⊙ , and 5 . 0 × 10 7 M ⊙ for the foreground DLAs in the line of sight. For a spherical geometry and with gas covering sizes × 2 -3 greater than the visible rest-frame UVgalaxy radii ( 34 ), this implies total H I gas masses of approx M H I ≃ 10 8 . 5 -10 9 M ⊙ . These HI masses are approx. × 2 -7 larger than the inferred stellar masses. Further, they imply short gas depletion times of 1-50 Myr due to their current active star formation. The inferred column densities also correspond to H I gas surface densities of log(Σ gas /M ⊙ pc -2 ) = 2 . 0 (CEERS-43833), log(Σ gas /M ⊙ pc -2 ) = 2 . 1 (CEERS-16943), and log(Σ gas /M ⊙ pc -2 ) = 2 . 3 (MACS0647-JD), respectively, larger than the gas surface densities observed in typical galaxies at z ≈ 0 ( 35 ).  \nBased on the SFRs and physical sizes of the galaxies, we derive SFR surface densities of log(Σ SFR /M ⊙ yr -1 kpc -2 ) = 1 . 3 ± 0 . 3 (CEERS-43833) and log(Σ SFR /M ⊙ yr -1 kpc -2 ) = 0 . 9 ± 0 . 3 (CEERS-16943), and approximately log(Σ SFR /M ⊙ yr -1 kpc -2 ) ∼ 2 ( 15 ) for MACS0647JD. These measurements place the young, z = 9 -11 star-forming galaxies among the population with the highest SFR surface densities observed in the local Universe ( 35 ). However, their associated gas surface densities place them offset to the Kennicutt-Schmidt (KS) relation, Σ SFR ∝ Σ n gas with n ≈ 1 . 4 . To be consistent with the KS relation, the gas surface densities should approach log(Σ gas /M ⊙ pc -2 ) ∼ 3 -4 , more than an order of magnitude larger than inferred from the DLA fit. This suggests a potentially large inherent gas mass in the galaxy itself, which does not contribute strongly to the DLA feature. This is similar to what is seen in I Zw 18, where the 21 cm inferred gas surface density ( 36 ) is nearly an order of magnitude higher than the DLA-inferred column density ( 28 ). Hence, the gas masses roughly estimated  \nabove for the H I gas in a shell around the galaxies must be considered as lower limits to the galaxies total gas masses, which may well be greater by an order of magnitude and thus more consistent with the fundamental KS relation.  \nTo quantify this and make predictions for the overall gas abundances and accretion histories of these three young, z = 9 -11 star-forming galaxies, we extract the mean gas masses of galaxies at their derived stellar masses and redshifts from the Astraues simulations framework ( 37 ) (see Supplementary Materials for details). The gas mass build-up for this subset of galaxies from z = 25 to z = 5 , i.e. 130 Myr to 1.2 Gyr after the Big Bang, in addition to the gas mass accretion at a given redshift, is shown in Fig. 3. From this simulation, we predict gas masses of M gas = 3 × 10 9 M ⊙ (CEERS-43833), M gas = 8 × 10 7 M ⊙ (CEERS-16943), and M gas = 8 × 10 8 M ⊙ (MACS0647-JD). We further infer gas fractions M gas /M ⋆ ≈ 5 , or f gas = M gas / ( M gas + M ⋆ ) ≈ 0 . 85 for all three galaxies, consistent with the lower bound obtained from the line-of-sight H I column densities. The total gas masses and gas fractions predicted in this simulation overall show gradual increases as a function of decreasing redshift, from M gas ∼ 2 × 10 7 ( M gas /M ⋆ = 100 ) at z > 20 to M gas = 10 9 -10 10 ( M gas /M ⋆ = 1 -2 ) at z = 5 .  \nOur results provide the first direct detection of the abundant neutral atomic H I gas reservoirs in the local environment of galaxies at the onset of the epoch of reionization at z > 8 . Circumstantial evidence for these abundant, neutral gas reservoirs in early galaxies has also been presented in recent literature based on strong absorption from singly-ionized UV absorption features ( 38 ) or an apparent flux deficit from the expected damped Ly α profile from a neutral IGM ( 15 ). Such abundant H I gas reservoirs in early galaxies are also observed indirectly through far-infrared gas tracers ( 39, 40 ). Our results imply that extreme H I gas masses might be common in galaxies at z = 9 -11 , in stark contrast to the scarcity of similar strong Ly α absorption features in Lyman-break galaxies at z ≈ 3 ( 8 ). Based on the large abundances  \nFigure 3: Gas build-up and accretion histories. ( Top: ) The predicted evolutionary tracks from the Astraeus simulations ( 37 ) of the build-up of M gas in the three z = 9 -11 galaxies based on their stellar masses as a function of redshifts. ( Bottom: ) The predicted tracks for the gas fraction, M gas /M ⋆ , for the three z = 9 -11 galaxies. The curves are generally consistent and suggest gas fractions M gas /M ⋆ ≈ 5 for the highz galaxies at their respective redshifts.  \n<!-- image -->  \nof neutral, pristine gas and the fact that these galaxies are surrounded by massive shells or extended sheets of neutral gas we argue that these DLA galaxy systems may be unique to the  \nreionization epoch. Consequently, galaxies like these cannot be leaking significant fractions of ionizing photons due to their extensive gas covering fractions, at least at their current evolutionary stages. These findings have important implications for understanding the process of reionization, and future studies of the incidence and demographics of DLAs surrounding highredshift galaxies will help distinguish between late and rapid reionization models ( 41,42 ), over a more smooth and early transition ( 43 ). The discovery of these extreme DLAs during reionization also has important implications for how robustly the degree of ionization of the IGM on large scales can be inferred. Clearly the damping wings caused by the neutral gas in the IGM will be masked by the large H I abundance in the local environment of galaxies at this epoch, leading to potential overestimations of the degree of neutrality in the bulk of the IGM ( 44 ) and overall photometric distances of galaxies in the early Universe. The young stellar populations and strong ionization parameters in the galaxies studied here, in addition to the potential more numerous population of H I-rich star-forming galaxies at z > 8 , may well indicate that we are starting to uncover the first epoch of galaxy formation.", 'References': "- 1. D. Kereˇs, N. Katz, D. H. Weinberg, R. Dav'e, Monthly Notices of the Royal Astronomical Society 363 , 2 (2005).\n- 2. J. Schaye, et al. , Monthly Notices of the Royal Astronomical Society 402 , 1536 (2010).\n- 3. D. P. Stark, Annual Review of Astronomy and Astrophysics 54 , 761 (2016).\n- 4. X. Fern'andez, et al. , Astrophysical Journal Letters 824 , L1 (2016).\n- 5. N. Maddox, et al. , Astronomy & Astrophysics 646 , A35 (2021).  \n- 6. A. M. Wolfe, E. Gawiser, J. X. Prochaska, Annual Review of Astronomy and Astrophysics 43 , 861 (2005).\n- 7. J. H. McKinney, et al. , Astrophysical Journal 874 , 52 (2019).\n- 8. A. E. Shapley, C. C. Steidel, M. Pettini, K. L. Adelberger, Astrophysical Journal 588 , 65 (2003).\n- 9. S. L. Finkelstein, et al. , Astrophysical Journal Letters 940 , L55 (2022).\n- 10. S. L. Finkelstein, et al. , Astrophysical Journal Letters 946 , L13 (2023).\n- 11. M. B. Bagley, et al. , Astrophysical Journal Letters 946 , L12 (2023).\n- 12. P. Arrabal Haro, et al. , arXiv e-prints p. arXiv:2303.15431 (2023).\n- 13. D. Coe, et al. , Astrophysical Journal 762 , 32 (2013).\n- 14. T. Y.-Y. Hsiao, et al. , arXiv e-prints p. arXiv:2210.14123 (2022).\n- 15. T. Y.-Y. Hsiao, et al. , arXiv e-prints p. arXiv:2305.03042 (2023).\n- 16. M. Tang, et al. , arXiv e-prints p. arXiv:2301.07072 (2023).\n- 17. K. E. Heintz, et al. , arXiv e-prints p. arXiv:2212.02890 (2022).\n- 18. R. L. Larson, et al. , arXiv e-prints p. arXiv:2211.10035 (2022).\n- 19. X. Fan, et al. , The Astronomical Journal 132 , 117 (2006).\n- 20. J. Miralda-Escud'e, Astrophysical Journal 501 , 15 (1998).\n- 21. T. Totani, et al. , Publications of the Astronomical Society of Japan 58 , 485 (2006).\n- 22. T. Tepper-Garc'ıa, Monthly Notices of the Royal Astronomical Society 369 , 2025 (2006).  \n- 23. J. Buchner, et al. , Astronomy & Astrophysics 564 , A125 (2014).\n- 24. F. Feroz, M. P. Hobson, M. Bridges, Monthly Notices of the Royal Astronomical Society 398 , 1601 (2009).\n- 25. E. Leonova, et al. , Monthly Notices of the Royal Astronomical Society 515 , 5790 (2022).\n- 26. L. Whitler, et al. , arXiv e-prints p. arXiv:2305.16670 (2023).\n- 27. P. Laursen, J. Sommer-Larsen, B. Milvang-Jensen, J. P. U. Fynbo, A. O. Razoumov, Astronomy & Astrophysics 627 , A84 (2019).\n- 28. S. Hernandez, et al. , Astrophysical Journal 892 , 19 (2020).\n- 29. H. Yang, et al. , Astrophysical Journal 844 , 171 (2017).\n- 30. A. E. Jaskot, T. Dowd, M. S. Oey, C. Scarlata, J. McKinney, Astrophysical Journal 885 , 96 (2019).\n- 31. R. Cen, T. Kimm, Astrophysical Journal Letters 801 , L25 (2015).\n- 32. J. Rosdahl, et al. , Monthly Notices of the Royal Astronomical Society 515 , 2386 (2022).\n- 33. J. Y. C. Yeh, et al. , Monthly Notices of the Royal Astronomical Society 520 , 2757 (2023).\n- 34. Y. Fudamoto, et al. , Astrophysical Journal 934 , 144 (2022).\n- 35. R. C. Kennicutt, N. J. Evans, Annual Review of Astronomy and Astrophysics 50 , 531 (2012).\n- 36. F. Lelli, M. Verheijen, F. Fraternali, R. Sancisi, Astronomy & Astrophysics 537 , A72 (2012).\n- 37. A. Hutter, et al. , Monthly Notices of the Royal Astronomical Society 503 , 3698 (2021).\n- 38. K. Boyett, et al. , arXiv e-prints p. arXiv:2303.00306 (2023).  \n- 39. K. E. Heintz, et al. , Astrophysical Journal Letters 934 , L27 (2022).\n- 40. K. E. Heintz, et al. , Astrophysical Journal Letters 944 , L30 (2023).\n- 41. R. P. Naidu, et al. , Astrophysical Journal 892 , 109 (2020).\n- 42. T.-Y. Lu, et al. , arXiv e-prints p. arXiv:2304.11192 (2023).\n- 43. S. L. Finkelstein, et al. , Astrophysical Journal 879 , 36 (2019).\n- 44. M. McQuinn, A. Lidz, M. Zaldarriaga, L. Hernquist, S. Dutta, Monthly Notices of the Royal Astronomical Society 388 , 1101 (2008).\n- 45. G. Brammer, Grizli: Grism redshift and line analysis software, Astrophysics Source Code Library, record ascl:1905.001 (2019).\n- 46. G. Brammer, msaexp: NIRSpec analyis tools, Zenodo (2022).\n- 47. Planck Collaboration, et al. , Astronomy & Astrophysics 641 , A6 (2020).\n- 48. Astropy Collaboration, et al. , Astronomy & Astrophysics 558 , A33 (2013).\n- 49. K. Horne, Publications of the Astronomical Society of the Pacific 98 , 609 (1986).\n- 50. P. Jakobsen, et al. , Astronomy & Astrophysics 661 , A80 (2022).\n- 51. S. Gonzaga, W. Hack, A. Fruchter, J. Mack, The DrizzlePac Handbook (2012).\n- 52. A. C. Carnall, R. J. McLure, J. S. Dunlop, R. Dav'e, Monthly Notices of the Royal Astronomical Society 480 , 4379 (2018).\n- 53. G. Bruzual, S. Charlot, Monthly Notices of the Royal Astronomical Society 344 , 1000 (2003).  \n- 54. G. J. Ferland, et al. , Revista Mexicana de Astronom'ıa y Astrof'ısica 53 , 385 (2017).\n- 55. S. Salim, M. Boquien, J. C. Lee, Astrophysical Journal 859 , 11 (2018).\n- 56. P. Kroupa, Monthly Notices of the Royal Astronomical Society 322 , 231 (2001).\n- 57. K. G. Iyer, et al. , Astrophysical Journal 879 , 116 (2019).\n- 58. J. Leja, A. C. Carnall, B. D. Johnson, C. Conroy, J. S. Speagle, Astrophysical Journal 876 , 3 (2019).\n- 59. S. Tacchella, et al. , Monthly Notices of the Royal Astronomical Society 522 , 6236 (2023).\n- 60. D. E. Osterbrock, G. J. Ferland, Astrophysics of gaseous nebulae and active galactic nuclei (2006).\n- 61. J. Kennicutt, Robert C., Annual Review of Astronomy and Astrophysics 36 , 189 (1998).\n- 62. M. W. Topping, et al. , Monthly Notices of the Royal Astronomical Society 516 , 975 (2022).\n- 63. Y. I. Izotov, G. Stasi'nska, G. Meynet, N. G. Guseva, T. X. Thuan, Astronomy & Astrophysics 448 , 955 (2006).\n- 64. R. L. Sanders, A. E. Shapley, M. W. Topping, N. A. Reddy, G. B. Brammer, arXiv e-prints p. arXiv:2303.08149 (2023).\n- 65. D. Langeroodi, et al. , arXiv e-prints p. arXiv:2212.02491 (2022).\n- 66. K. Nakajima, et al. , arXiv e-prints p. arXiv:2301.12825 (2023).\n- 67. M. Curti, et al. , arXiv e-prints p. arXiv:2304.08516 (2023).\n- 68. D. Watson, Astronomy & Astrophysics 533 , A16 (2011).  \n- 69. K. D. Gordon, G. C. Clayton, K. A. Misselt, A. U. Landolt, M. J. Wolff, Astrophysical Journal 594 , 279 (2003).\n- 70. P. De Vis, et al. , Astronomy & Astrophysics 623 , A5 (2019).\n- 71. D. Foreman-Mackey, D. W. Hogg, D. Lang, J. Goodman, Publications of the Astronomical Society of the Pacific 125 , 306 (2013).\n- 72. C. L. Steinhardt, V. Kokorev, V. Rusakov, E. Garcia, A. Sneppen, arXiv e-prints p. arXiv:2208.07879 (2022).\n- 73. N. Byler, J. J. Dalcanton, C. Conroy, B. D. Johnson, Astrophysical Journal 840 , 44 (2017).\n- 74. A. Raiter, D. Schaerer, R. A. E. Fosbury, Astronomy & Astrophysics 523 , A64 (2010).\n- 75. N. R. Tanvir, et al. , Monthly Notices of the Royal Astronomical Society 483 , 5380 (2019).\n- 76. O. E. Hartoog, et al. , Astronomy & Astrophysics 580 , A139 (2015).\n- 77. A. Saccardi, et al. , Astronomy & Astrophysics 671 , A84 (2023).\n- 78. R. A. Simcoe, et al. , arXiv e-prints p. arXiv:2011.10582 (2020).\n- 79. A. Hutter, et al. , arXiv e-prints p. arXiv:2209.14592 (2022).\n- 80. Planck Collaboration, et al. , Astronomy & Astrophysics 641 , A6 (2020).\n- 81. D. Nelson, et al. , Monthly Notices of the Royal Astronomical Society 490 , 3234 (2019).\n- 82. I. Langan, D. Ceverino, K. Finlator, Monthly Notices of the Royal Astronomical Society 494 , 1988 (2020).\n- 83. C. C. Lovell, et al. , Monthly Notices of the Royal Astronomical Society 500 , 2127 (2021).  \n- 84. G. Ucci, et al. , Monthly Notices of the Royal Astronomical Society 518 , 3557 (2023).\n- 85. M. Curti, F. Mannucci, G. Cresci, R. Maiolino, Monthly Notices of the Royal Astronomical Society 491 , 944 (2020).", 'Acknowledgements': 'We would like to thank Johan P. U. Fynbo and Mark Dickinson for their constructive and enlightening comments on the interpretations of the results presented in this work. We further acknowledge the significant work and effort of the CEERS collaboration in obtaining parts of the observations presented here and are grateful that their early data is publicly available. K.E.H. acknowledges support from the Carlsberg Foundation Reintegration Fellowship Grant CF21-0103. A.H. acknowledges support by the VILLUM FONDEN under grant 37459. C.M. acknowledges support by the VILLUM FONDEN under grant 37459 and the Carlsberg Foundation under grant CF22-1322. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant No. 140. RPN acknowledges funding from JWST programs GO-1933 and GO-2279. 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. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST.', 'Author contributions': 'K.E.H. and D.W. drafted the manuscript and led the main analysis. G.B. extracted and reduced the observational data. S.V. performed the IGM and DLA modelling. A.H. provided and interpreted the Astraeus simulations output. V.B.S. performed the SED modelling. All authors contributed to the interpretation and analysis of the results presented in this work.', 'Competing interests': 'The authors declare no competing interests.', 'Data and materials availability': 'The JWST imaging and spectroscopic data are publicly available on the JWST MAST archive at https://mast.stsci.edu . The data has been processed using public software codes; grizli v. 1.8.3, ( 45 ) and MsaExp v. 0.6.7, ( 46 ). The reduced data are available from the corresponding author upon reasonable request.', 'Supplementary Materials': "Cosmology. Throughout this paper we assume the concordance, flat Λ CDM cosmological model with H 0 = 67 . 4 kms -1 Mpc -1 , Ω m = 0 . 315 , and Ω Λ = 0 . 685 ( 47 ). Cosmological measurements such as the luminosity distances d L , age of the universe at a given z , and transverse proper distances are computed using the cosmology distance calculator from Astropy ( 48 ).  \nObservations and data reduction. This work is mainly based on spectroscopic observations collected through the dedicated JWST Director's Discretionary time program DD-2750 (PI: P. Arrabal Haro) and the General Observer program GO 1433 (PI: D. Coe). The spectroscopic data obtained of CEERS-43833 and CEERS-16943 (MSA ID's 28 and 1, respectively) have previously been presented in earlier works ( 12 ), and likewise for MACS0647-JD ( 15 ). For consistency and internal homogeneity, we here reprocess the full spectroscopic dataset, using the custom-made pipeline MsaExp v. 0.6.7 ( 46 ). This code utilizes the Stage 2 output from the MASTJWSTarchive and performs standard wavelength, flat-field and photometric calibrations on the individual NIRSpec exposure files. This is done with the relevant JWST reference files associated with the CRDS context jwst1027.pmap . MsaExp further corrects for the 1 /f noise and the bias levels in individual exposures. We generate the full combined 2D spectra from these output files and optimally extract the 1D spectra using an inverse-weighted sum of the 2D spectra in the dispersion directions ( 49 ). We scale the flux densities of the extracted 1D spectra to match the integrated flux within the available JWST/NIRcam passband in each filter, using a wavelength-dependent polynomial function. This is to improve the absolute flux calibration of the spectra and take into account any potential slit-losses. The required correction factors are typically less than ≈ 10% . The NIRSpec prism observations all have wavelength coverage from 0 . 7 µ mto 5 . 3 µ m, with a varying spectral resolution from R≈ 50 in the blue end to R≈ 400 in the red end ( 50 ). The reduced and photometrically-calibrated spectra of the three  \ngalaxies at z = 8 . 8 -11 . 4 analyzed here are shown in Fig. 1.  \nWe further include JWST/NIRCam imaging obtained through the CEERS survey (ERS1345, PI: S. Finkelstein) ( 9-11 ), in which CEERS-43833 and CEERS-16943 were first identified as a high-redshift galaxy candidates (the latter dubbed Maisie's galaxy ( 9 )). MACS0647JD was discovered as part of the Cluster Lensing And Supernova survey with Hubble (CLASH) ( 13 ), and was later observed with JWST/NIRCam (GO 1433, PI: D. Coe). These observations revealed a complex morphology of the multiply-lensed source ( 14 ). For each source, we adopt the photometric measurements from the catalog presented by Brammer et al. (in prep.), currently available online 1 . This repository includes reduced images of the raw data obtained from MAST using the public software package grizli , which provides astrometric calibrations based on the Gaia -DR3 catalog, masks imaging artifacts, and drizzles the images to a common pixel scale of 0 . '' 04 / pixel using astrodrizzle ( 51 ). The final photometric catalog includes measurements that are based on the most recent JWST/NIRCam photometric zeropoints and have been corrected for extinction due to dust in the Galaxy. For the modelling of the galaxy SEDs we use the aperture-matched photometry from this catalog, derived using a circular aperture with diameter 0 . '' 5 .  \nMorphology and galaxy sizes. Weadopt the size measurements of MACS0647-JD from ( 15 ), assuming that the galaxy spans the edges of the A and B components, yielding an effective diameter of ∼ 0 . 5 kpc. For CEERS-43833 we measure full-width-half-maximum (FWHM) angular major and minor axes of 0 . '' 21 × 0 . '' 19 (in F200W, i.e. rest-frame UV), corresponding to 0 . 97 × 0 . 88 kpc at z = 8 . 7622 . For CEERS-16943, we measure FWHM angular major and minor axes of 0 . '' 17 × 0 . '' 14 (F200W), corresponding to 0 . 66 × 0 . 54 kpc at z = 11 . 409 . We estimate the physical sizes enclosing half of the total rest-frame UV luminosity estimated from the Gaussian modelling (half-light semi axes = a maj , min = FWHM maj , min / 2 ), yielding effective  \nareas of the projected ellipses of 0 . 67 kpc -2 (CEERS-43833) and 0 . 28 kpc -2 (CEERS-16943), respectively.  \nSED modelling. We adopt the outputs from the SED modelling of the physical properties of the MACS0647-JD complex as a whole, as presented in ( 15 ), and summarized in Table 1. Their analysis reveals a low-mass, log( M ⋆ /M ⊙ ) = 8 . 1 ± 0 . 3 , galaxy with SFR SED = 8 ± 3 M ⊙ yr -1 , and a blue continuum slope with log U = -1 . 9 ± 0 . 2 and A V , SED = 0 . 07 ± 0 . 06 mag. To infer the physical properties of CEERS-43833 and CEERS-16943, we model their SEDs using the Python code Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation ( Bagpipes ) ( 52 ), which incorporates the stellar population models from Bruzual & Charlot ( 53 ). We perform spectrophotometric fitting of the two galaxy SEDs, fixing the redshifts to z spec inferred from the nebular emission lines. The emission lines are added to the model spectra in Bagpipes using a grid computed from Cloudy models ( 54 ). We allow the ionization parameter, U , to vary between -3 < log U < -1 and assume the attenuation curve from Salim et al. ( 55 ) to account for the typically higher ionization parameters and the more steep attenuation curves of highz galaxies. Throughout, we assume a Kroupa initial mass function (IMF) ( 56 ). The best-fit SED models for CEERS-43833 and CEERS-16943 are shown in Fig. S1.  \nFor simple, constant SFHs we infer young ( ≲ 100 Myr) stellar populations on average in both galaxies. We caution, however, that the inferred ages are heavily dependent on the assumed SFH parameterization. The accurate spectrophotometric fitting allow us to incorporate more flexible, non-parametric SFHs for each galaxy which more accurately capture the stellar ages and mass build-up ( 57 ). This effectively increases the inferred stellar masses (by 0.20.3 dex for these two galaxies) compared to the best-fit models assuming constant SFHs, as also demonstrated in previous works ( 58, 59 ). This parameterization fits for quartiles of total mass formed, and includes a prolonged, less active period of star formation to take into account pos-  \nble hidden mass from an older stellar population due to 'outshining' from a potential younger, burstier star-forming population. Table 1 reports the median and 16th and 84th percentiles from the resulting posterior distributions on the stellar mass, SFR, dust attenuation, and ionization parameter for each galaxy. The results for CEERS-16943 are overall consistent with those obtained by Finkelstein et al. ( 9 ), which were, however, based on the photometric measurements only.  \nEmission line measurements. The emission line detections and spectroscopic redshift for CEERS-43833 and MACS0647-JD have previously been reported based on several nebular emission lines ( 12, 15 ), but for consistency we here rederive the line fluxes and redshifts. We measure z spec = 8 . 7622 ± 0 . 0002 and z spec = 10 . 170 ± 0 . 003 , respectively, as listed in Table 1, see also Fig. S2. The line fluxes are derived by jointly modelling the local continuum and the emission lines represented by Gaussian profiles, tying the redshift and width of the lines across transitions, but allow the individual line fluxes to vary. The measured line fluxes are summarized in Table S1.  \nThe spectroscopic redshifts of CEERS-16943 have also previously been reported ( 12 ), however, only based on the approx. onset of the Ly α break. As demonstrated in this work, this will overestimate the redshift due to broadening of Ly α from the strong DLA feature. We here determine the spectroscopic redshift independently and provide the line fluxes for the new detected emission lines in the spectrum. We measure z = 11 . 409 ± 0 . 001 for CEERS-16943 based on the detection of [O II] λ 3727 , [Ne III] λ 3870 , and He I λ 3889 , see Fig. S1. This is the to-date most distant detection of nebular emission lines and thus also the earliest detection of single metals or elements heavier than hydrogen. We note that in the previous analysis of CEERS-16943 ( 12 ) the feature at the location of [O II] was proposed to be caused by image defects. However, based on our independent reduction and the detection of two other features at the same redshift, we argue that this emission line feature is real.  \nPhysical Properties. In the typical case B recombination scenario, we expect Balmer line ratios H γ/ H β = 0 . 47 and H δ/ H β = 0 . 26 for electron temperatures T e = 10 4 K ( 60 ), which is consistent within 1 σ of our measurements H γ/ H β = 0 . 48 ± 0 . 08 and H δ/ H β = 0 . 21 ± 0 . 07 for CEERS-43833. This implies negligible extinction due to dust, as also supported by the SED modelling. For MACS0647-JD and CEERS-16943 we are not able to measure the Balmer decrement, but their SEDs imply similar low dust attenuation.  \nFor CEERS-43833, we derive the SFR from H β ( 61 ), where SFR = 5 . 5 × 10 -42 L H β (erg / s) × f H α/ H β , assuming a Kroupa IMF ( 56 ). Here, f H α/ H β = 2 . 87 is the predicted ratio from the Case B recombination scenario at T e = 10 4 K. For MACS0647-JD and CEERS-16943, where H β is redshifted outside of the JWST/NIRSpec wavelength coverage, we use the relation between SFR and [O II] λ 3727 line luminosity ( 61 ); SFR = 1 . 0 × 10 -41 L [O II ] (erg / s) , again corrected to a Kroupa IMF. The exact choice of the IMF introduces a 0.3 dex systematic uncertainty on these measurement, which we include in our derivation. The SFRs range from 1 to 15 M ⊙ yr -1 , and are summarized in Table 1. Combined with the stellar mass inferred from the SED fitting, we compute specific SFRs in the range log(sSFR/yr -1 ) = -7 . 55 ± 0 . 30 (CEERS-43833) to log(sSFR/yr -1 ) = -8 . 25 ± 0 . 30 (CEERS-16943). These are consistent with the typical starforming galaxy population in the reionization epoch ( 17,62 ).  \nWe infer the gas-phase metallicities for each galaxy based on strong-line diagnostics, since we do not detect the auroral [O III] λ 4363 emission line enabling direct T e -based measurements ( 63 ), except for in the spectrum of MACS0647-JD (see also ( 15 )). This spectrum, however, does not cover the [O III] λλ 4960 , 5008 transition doublet as well, required to determine the metallicity through the direct T e -method. We use the recent set of strong-line calibrations based on the direct metallicities inferred from a large sample of galaxies at z = 2 -9 spectroscopically observed with JWST ( 64 ). For CEERS-43833, we report the metallicity inferred using the O3 = [O III] λ 5008 /H β line ratio, 12 + log (O/H) = 7 . 46 ± 0 . 09 , as this calibration is  \nless affected by dust extinction and the ionization state of the gas, in addition to any potential issues with the flux calibration as a function of wavelength. Further, it shows the least intrinsic scatter as derived by Sanders et al. ( 64 ). For MACS0647-JD and CEERS-16943, where these features are redshifted out of the NIRSpec wavelength coverage, we instead adopt the Ne3O2 = [Ne III] λ 3870 / [O II] λ 3727 line ratio, yielding 12 + log (O/H) = 7 . 40 ± 0 . 3 and 7 . 70 ± 0 . 30 , respectively. The inferred metallicity of MACS0647-JD is consistent with previous estimates ( 15 ), based on various methods and line diagnostics. Using the Ne3O2 calibration to infer the metallicity of CEERS-43833 yields a consistent result to the O3 calibration within the uncertainties. The inferred metallicities are in the range 12 + log (O/H) = 7.40 - 7.70, i.e. 5-10% solar, as listed in Table 1. In Fig. S3, we show the physical properties of the three young, z = 9 -11 galaxies, compared to the typical population of star-forming galaxies at z = 7 -10 ( 17 ), in terms of their SFRs, stellar masses, and gas-phase metallicities. Their inferred physical properties are consistent with the established mass-metallicity (MZ) relation at z > 7 ( 17, 65, 66 ), and furthers shows a significant offset by 0.3-0.5 dex at 2 -5 σ from the fundamental-metallicity relation (FMR) at z ≈ 0 as also observed for other star-forming galaxies at z > 7 ( 17,66,67 ).  \nA limit on the metallicity of the DLAs can be determined from their lack of significant dust extinction. At a Milky Way dust-to-gas (DTG) ratio ( 68 ), such DLAs would produce approximately 50 magnitudes of extinction at the Ly α wavelength using an SMC-like extinction curve ( 69 ). Using the inferred galaxy attenuation as an estimate of the extinction, we find DTG ratios consistent with the power-law relation of ( 70 ) for our inferred emission-line metallicities. This also allows us partially to address the question of whether such extreme DLAs do exist at lower redshift but are not detected due to dust extinction. Assuming such galaxies could be detected with at most one magnitude of extinction at Ly α and a column density of N H I ≃ 10 22 cm -2 , this corresponds to DTG ratios ≲ 2 . 5 × 10 -4 . Such low DTG ratios are only observed  \nin galaxies with metallicities of 12+log(O/H) ≲ 8 . 15 . This suggests that the absence of these extreme DLAs at lower redshift is real and due to a lack of huge foreground H I reservoirs rather than dust bias.  \nThe damped Ly α absorption feature. We model the broad DLA features in each galaxy at z = 9 -11 by jointly fitting the optical depth from the Gunn-Peterson trough due to absorption from H I in the IGM τ IGM ( z ) and the H I located in the immediate surroundings of the galaxies, τ ISM , as described in the Main Text. This allows us to independently model the neutral hydrogen fraction, x H I, in the IGM and the H I column density in the near-proximity to the galaxies, most likely in the interstellar or circumgalactic medium. We use a set of galaxy template spectra designed to match the blue rest-frame UV colors of galaxies at z > 8 ( 18 ) as the intrinsic continuum galaxy models. We leave the visual extinction A V , the average IGM neutral fraction, x H I, and the H I column density, N H I, at z = z gal as free parameters. We recover extremely strong DLAs in all three z = 9 -11 galaxies, with column densities N H I > 10 22 cm -2 . The modelling is largely insensitive to x H I, as the DLA feature dominates the observed line profile in the three sources (see Fig. 1). From the modelling, we infer lower bounds on x H I from the 68% and 95% highest density intervals of the posterior distributions, x H I = 0 . 35 -0 . 4 and x H I = 0 . 06 -0 . 08 , respectively.  \nTo verify the robustness of the models, we also run a Markov Chain Monte Carlo (MCMC) algorithm ( 71 ) on the spectra. This fitting method finds similar values for the neutral hydrogen fraction, x H I, of the IGM and for the H I column density. We further consider an expanded set of intrinsic galaxy models, including models with bursty SFHs and continous SFHs on longer timescales (10-30 Myr) from ( 18 ), and models with more physically motivated IMFs for galaxies at z = 8 -12 ( 72 ). The variations in these intrinsic galaxy models introduce at most a systematic uncertainty of 0.1 dex on the derived H I column densities, which we propagate with the statistical uncertainty to the total uncertainty listed in Table 1. This is because the strong  \nDLA feature dominates the spectral shape around Ly α over the more subtle variations in the various galaxy continuum emission models.  \nTo exclude any potential false-positive detections of DLAs due to instrumental issues with the flux or wavelength calibration of NIRSpec, we compare CEERS-43833 to a set of other highz galaxies in Fig. S2 with similar spectroscopic observations (from ( 17 )). These comparison galaxies are selected to be at similar redshifts, with stellar masses and metallicities ranging from M ⋆ = 10 8 . 6 -10 10 M ⊙ and 12 + log(O / H) = 7 . 4 -8 . 0 , and have equally high S/N in the spectral regions around Ly α . We convert and scale the spectroscopic redshifts to z = 8 . 7622 and normalize the continuum for each galaxy to that of CEERS-43833 to improve the visual comparison. It is evident that none of these comparison galaxies exhibit similar broadening of the Ly α trough, beyond that expected for the IGM at x H I = 0 . 1 -1 . In Fig. S4 we also overplot the various IGM models convolved with the wavelength-dependent spectral resolution of the JWST/NIRSpec prism observations. While the IGM models are less easily discerned, the broad DLA feature is still readily recovered even in the low spectral-resolution observations. This excludes the broad Ly α feature as being induced by instrumental broadening effects.  \nWe also consider other physical effects that might mimic the broad Ly α troughs observed in the galaxy spectra. For instance, nebular continuum emission models with varying stellar mass ranges have been shown to produce similar damped Ly α profiles ( 73 ) as observed. However, this effect is only recovered for stellar populations dominated by older, low-mass ( < 5 M ⊙ ) stars in these models, which we argue are unlikely to be the case for the three highz galaxies studied here, and is also inconsistent with the inferred blue spectral slopes and young stellar populations. Further, 2 γ continuum emission may potentially imprint a distinct continuum feature near Ly α as well ( 74 ). However, the probability of these transitions only becomes significant in extremely metal-poor galaxies (less than × 10 the metallicities of the z = 9 -11 galaxies), where the physical state of the star-forming gas would depart from Case B, so we  \nargue that this mechanism is also unlikely to explain the observed broad Ly α troughs.  \nThe extended DLAs covering the UV emission of the three young, z = 9 -11 galaxies are the strongest observed to-date in emission-selected star-forming galaxies, considering some of the most metal-poor galaxies at z ≈ 0 ( 7,28 ) or Lyman-break galaxies at z ≈ 3 ( 8 ). Further, only about ≲ 15% of γ -ray burst absorption systems at redshifts z = 2 -6 show H I column densities exceeding N H I > 10 22 cm -2 ( 75 ). Moreover, the most distant DLAs previously detected have been observed out to redshifts z ≈ 6 in γ -ray burst absorbers ( 76,77 ), and indirectly at z = 6 . 84 in a single quasar absorption system ( 78 ). For Lyman-break galaxies, the DLA feature has only been seen out to z ≈ 3 ( 8 ) due to the lacking NIR sensitivity of ground-based telescopes. With these JWST observations, we are now able to push this observational barrier and detect DLAs out to redshifts beyond z > 11 , i.e. within 400 Myr after the Big Bang at almost a factor of two larger redshifts than previously possible. In Fig. S5 we compare the H I column densities of the z = 9 -11 galaxies to representative samples of metal-poor, z ≈ 0 'highz analogs' ( 7, 28 ) and the most recent collection of γ -ray burst absorbers expanded upon the initial sample presented in ( 75 ). We further mark the maximum H I column densities based on the lowest Ly α equivalent width (EW) sources in the z ≈ 3 Lyman-break galaxy sample by ( 8 ), assuming that N H I = 1 . 88 × 10 18 EW 2 Ly α cm -2 in the Ly α damping regime.  \nIGM modelling. To include a more physically motivated transmission of the IGM, we consider a typical sightline from a galaxy at z ≃ 8 . 8 (matching approximate CEERS-43833), based on high-resolution cosmological hydro-simulations, post-processed with radiative transfer of ionizing UV and Ly α photons ( 27 ). The transmission T ( λ ) as a function of wavelength is calculated as the median value (and 16/84 and 5/95 percentiles) of 10 4 sightlines emanating from ∼ 300 galaxies in a zoomed-in region, 27 comoving Mpc across, resolved down to 10 pc scales. Because of the large neutral fraction of the IGM, the scattering of Ly α photons into the line of sight cannot be neglected out to several times the virial radius of a galaxy, and hence the  \nsightlines start at 10 r vir . At this point the radiative transfer is dominated by scattering out of the line of sight. More details on both the hydro-simulations and the radiative transfer can be found in ( 27 ).  \nIn Fig. S6, we show the median and the '1 σ ' (given by the 16th and the 84th percentiles) and the '2 σ ' (given by the 5th and the 95th percentiles) distributions of these transmission curves shifted to the redshift of the galaxy at z = 8 . 7622 . Evidently, none of these modelled galaxy sightlines are able to reproduce the large Ly α damping observed in this particular galaxy spectrum, supporting a large abundance of H I in the interstellar or circumgalactic medium of the galaxy. We caution, however, that the models in ( 27 ) assume a relatively low IGM neutral gas fraction, x H I = 0 . 13 , so the Ly α transmission curves should effectively be treated as lower limits for these particular models.  \nAstraeus simulation. The Astraeus simulations included in this work, couples a semianalytical galaxy evolution model to a semi-numerical reionization scheme and runs on the outputs of a dark-matter-only N -body simulation. It follows the key baryonic processes of gas accretion, gas and stellar mass being brought in by mergers, star formation, supernovae feedback, metal and dust enrichment, and radiative feedback from reionization, as well as the spatial ionization of the IGM (see ( 37, 79 ) for details). In this paper, we analyse the MHDEC simulation (described in ( 79 )), which assumes the escape fraction of ionising photons to decrease with rising halo mass. This simulation reproduces the observed UV luminosity and stellar mass functions at z = 5 -10 as well as the electron optical depth and the IGM's neutral hydrogen fraction at z < 6 inferred from Planck ( 80 ) and quasar absorption line measurements, respectively.  \nTable S1. Line flux measurements reported in units of 10 -19 erg s -1 cm -2 ˚ A -1 . The measurements for MACS0647-JD have not been corrected for the magnification factor µ = 8 ± 1 .  \nFig. S1. Imaging and best-fit SED models. In the top panels for CEERS-43833 and CEERS16943 are shown the reduced, drizzled images of the source in each available JWST/NIRCam filter. In the bottom panels are shown the reduced and photometrically-calibrated 1D spectra (black curves) with the associated error spectra (grey-shaded regions). The best-fit SED models from the spectrophotometric fitting using Bagpipes are overplotted in orange.  \n<!-- image -->  \nFig. S2. Detected line emission features. The reduced and photometrically-calibrated 1D spectra for each galaxy at z = 9 -11 (marked in the top left of each panel) are shown in black with the associated error spectra shown by the grey-shaded regions. The best-fit local continuum emission and Gaussian line profiles are overplotted in red.  \n<!-- image -->  \nFig. S3. Physical properties of the young, z = 9 -11 star-forming galaxies. In the left panel, the inferred stellar masses and gas-phase metallicities derived in this work (red star symbols) are compared to individual observations (black) and the best-fit mass-metallicity relation of main-sequence star-forming galaxies at z = 7 -10 ( 17 ), in addition to a suite of predictions from simulations at similar redshifts ( 81-84 ), as indicated by the labels in the bottom right. In the right panel, is shown the fundamental-metallicity relation again including individual starforming galaxies at z = 7 -10 ( 17 ), and predictions from the Astraeus simulation at similar redshifts ( 84 ). Both panels show the best-fit mass-metallicity (left) and fundamental-metallicity relation (right) of galaxies at z ≈ 0 ( 85 ) for comparison.  \n<!-- image -->  \nFig. S4. Ly α comparison plots of star-forming galaxies at z ≈ 8 . CEERS-43833 (black) at z = 8 . 7622 with a clear DLA feature is compared to a set of four galaxies from the literature ( 17 ) at similar redshifts and with equal or higher spectral quality. The IGM and DLA models are identical to Fig. 1, but here convolved with the spectral resolution of the JWST/NIRSpec prism observations. The models and comparable spectra exclude the broad Ly α feature observed in CEERS-43833, and the two other galaxies at z > 9 , as being induced by instrumental broadening effects.  \n<!-- image -->  \nFig. S5. H I column density distribution as a function of redshift. The symbol notation follow Fig. 2, but here we include additional constraints on the maximum N H I observed for Lyman-break galaxies at z ≈ 3 ( 8 ) shown as the green triangle and the expanded set of γ -ray burst H I absorbers at z ∼ 2 -6 ( 75 ) shown as the orange circles. The DLA demarcation region, N H I > 2 × 10 20 cm -2 , is marked by the light-grey area and the dark-grey region highlights the lower bound on N H I > 10 22 cm -2 , corresponding to an optical depth of τ > 10 5 at the Lyman limit, in the three z = 9 -11 galaxies analyzed in this work.  \n<!-- image -->  \nFig. S6. IGM transmission curves from model galaxy sightlines. The normalized spectrum of CEERS-43833 (black) and associated error spectrum (grey) is compared to the median (blue line) and 1 σ and 2 σ distributions of IGM transimissions curves modelled from typical galaxy sightlines at z ≃ 8 . 8 , out to 10 r vir , based on cosmological hydro-simulations and theoretical radiative transfer of ionizing UV and Ly α photons ( 27 ).  \n<!-- image -->"}