The JWST-PRIMAL Legacy Survey

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¹¹1https://dawn-cph.github.io/dja. 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 different programs. With DJA-Spec we reduce and extract all of the spectra from this diverse array of observational programs with a single standardized pipeline.

All data processing is performed with the publicly available Grizli (Brammer, 2023a)²²2https://github.com/gbrammer/grizli and MSAExp (Brammer, 2023b)³³3https://github.com/gbrammer/msaexp software 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.⁴⁴4jwst pipeline version 1.12.5; CRDS_CONTEXT = jwst_1180.pmap. The pipeline is executed with the default parameters, but with snowblind⁵⁵5https://github.com/mpi-astronomy/snowblind 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 $\approx 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 flux- and 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 multi-extension FITS files, and it is here that the subsequent processing with MSAExp departs significantly from the standard jwst pipeline.

The sky background of each source spectrum is removed by taking differences of the 2D cutout spectra taken at different “nod” positions of the telescope. Typically there are three nod offsets that shift the source by $0\aas@@fstack{\prime\prime}5\approx 5\leavevmode\nobreak\ \mathrm{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 offset 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 wavelength-dependent 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 mask. 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 different offset 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 effort to spectroscopically characterize galaxies at $z>5$ observed with NIRSpec Prism.

Currently (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 high-resolution gratings. A full overview can be found on the dedicated webpage⁶⁶6https://s3.amazonaws.com/msaexp-nirspec/extractions/nirspec_graded_v2.html, 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.

In this work, we focus on the JWST/NIRSpec Prism spectroscopic data from DJA due to the higher continuum sensitivity which is essential for quantifying the evidence for Ly$\alpha$ emission or absorption damping wings. JWST’s NIRSpec Prism configuration covers the entire near-infrared passband from 0.6 $\mu$m to 5.3 $\mu$m, with a varying spectral resolution from $\mathcal{R}\approx 50$ in the blue end to $\mathcal{R}\approx 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.

The main objective of this paper is to accurately model the Ly$\alpha$ damping wings and/or measure strong Ly$\alpha$ 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:

1. 1.

   A robust, spectroscopic redshift measurement at $z>5.5$ from a minimum of one emission line detected at $>3\sigma$ (in addition to the Ly$\alpha$ break).

2. 2.

   An integrated signal-to-noise over the redshifted Ly$\alpha$ and rest-frame UV ($1550\,\AA$) regions of ${\rm S/N}>3$.

3. 3.

   Flux density measurements over the entire NIRSpec Prism wavelength coverage from $0.6-5.3\,\mu$m.

These criteria ensure that the redshift of the foreground H i gas absorber, from interstellar to intergalactic scales, and the shape of the Ly$\alpha$ line profile can be accurately measured. Further, the requirement 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.

In total, 494 sources from DJA-Spec meet these criteria. Fig. 1 show the absolute UV magnitude, $M_{\rm 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.

During the reionization epoch at $z\gtrsim 6$, several factors add to the line shape of Ly$\alpha$ on integrated galaxy spectra: They can show strong Ly$\alpha$ 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 et al., 2016, 2018; Fujimoto et al., 2023b; Umeda et al., 2023; Hayes & Scarlata, 2023), strong damped Ly$\alpha$ 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 effects are degenerate, so we have to disentangle them statistically. For this, we define a new simple diagnostic, which we denote the Ly$\alpha$ damping parameter:

which is the equivalent width (EW) of the transmitted flux density, $F_{\lambda}$ over the wavelength region covered by the instrumentally-broadened Ly$\alpha$ transition. The continuum flux, $F_{\rm cont}$, over the same region is estimated by extrapolating the rest-frame UV slope, $\beta_{\rm UV}$, which is derived directly in each spectra from rest-frame $1400-2600\,\AA$ (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 $\lambda_{\rm Ly\alpha,low}=1180\,\AA$ to $\lambda_{\rm Ly\alpha,up}=1350\,\AA$ (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.

A set of physical models, shown as imprints on the Ly$\alpha$ transmission curves, are visualized in Fig. 2 including: varying ionized bubble sizes from $R_{\rm b}=1-50$ cMpc, neutral hydrogen fractions of the IGM from $x_{\rm HI}=0.1-1.0$, and DLA H i column densities $N_{\rm HI}=10^{21}-10^{22}$ cm^-2, all affecting the shape of the Ly$\alpha$ 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_{\rm IGM,upper}<z_{\rm 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$\alpha$ damping parameter $D_{\rm Ly\alpha}$ is listed for each model. The UV bubble and DLA models all assume a benchmark IGM fraction of $x_{\rm H\textsc{i}}=0.5$, and may therefore vary with the same scatter in $D_{\rm Ly\alpha}$ as observed for the IGM-only models.

We derive damping parameters for the full sample in the range $D_{\rm Ly\alpha}=-400$ to $150\,\AA$. With this, we define four parameter spaces probing different physical regimes: (i) cases with excess continuum flux suggesting contributions from ionized bubbles or weak LAEs have $D_{\rm Ly\alpha}=0$–35 $\AA$; (ii) the effect on the Ly$\alpha$ damping wings of an increasing neutral IGM results in $D_{\rm Ly\alpha}=35$–50 $\AA$; (iii) excess damped Ly$\alpha$ absorption from local H i (with $N_{\rm HI}>10^{21}$ cm^-2) have $D_{\rm Ly\alpha}>50\,\AA$; and (iv) strong LAEs are typically observed with $D_{\rm Ly\alpha}<0\,\AA$. Fig. 3 shows three example spectra from the JWST-PRIMAL sample. One shows an example of a broad Ly$\alpha$ damping wing consistent with local H i absorption with $D_{\rm Ly\alpha}=61.9\pm 3.6\,\AA$, consistent with $N_{\rm HI}\approx 10^{22}\,$cm^-2. Another shows strong Ly$\alpha$ emission with a measured $D_{\rm Ly\alpha}=-60.1\pm 8.0\,\AA$.

This 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$\alpha$ 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_{\rm Ly\alpha}$. Further, as demonstrated in the bottom panel of Fig. 3, rare cases ($\lesssim 1\%$ of the sample) that show both a substantial damped Ly$\alpha$ broadening and simultaneously a strong Ly$\alpha$ emission line will not be accurately identified, here with $D_{\rm Ly\alpha}=19.3\pm 0.7\,\AA$, misleadingly suggesting mild excess Ly$\alpha$ 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$\alpha$ emission, thereby underestimating the total Ly$\alpha$ photon output of the target source (Jiang et al., 2023). The $D_{\rm Ly\alpha}$ parameter is thus most powerful in identifying the most extreme cases, such as strong Ly$\alpha$ emitters or galaxy DLAs with $N_{\rm HI}>10^{21}$ cm^-2, which uniquely have $D_{\rm Ly\alpha}<0\,\AA$ and $D_{\rm Ly\alpha}>50\,\AA$, respectively.

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] $\lambda\lambda 3726,3729$ and [O iii] $\lambda\lambda 4960,5008$ doublets, [Ne iii] $\lambda 3870$, He i $\lambda 3889$, the [S ii] $\lambda\lambda 6718,6732$ doublet, and the Balmer lines H$\alpha$, H$\beta$, H$\gamma$, and $H\delta$. The [O ii] $\lambda\lambda 3726,3729$ and [S ii] $\lambda\lambda 6718,6732$ doublets are unresolved at all wavelengths in the NIRSpec Prism spectra, whereas the [O iii] $\lambda\lambda 4960,5008$ doublet is marginally resolved at $z\approx 6$ and fully resolved at $z\gtrsim 7.5$ due to the increasing spectral resolution with wavelength in this particular configuration. H$\alpha$ and [N ii] $\lambda 6585$ are unresolved at all redshifts. [Ne iii] $\lambda 3870$ and He i $\lambda 3889$ are also generally unresolved, except for the highest redshifts at $z\gtrsim 10$. Due to the broad wavelength coverage ($0.6-5.3\mu$m), H$\alpha$ and [S ii] $\lambda\lambda 6718,6732$ can be detected up to $z\approx 7$, H$\beta$ and the [O iii] $\lambda\lambda 4960,5008$ doublet up to $z\approx 9.5$, and Ly$\alpha$ can be detected at all redshifts from $z\gtrsim 5.5$. The auroral [O iii] $\lambda 4363$ line is also in the majority of cases blended with H$\gamma$ but becomes resolved at $z\gtrsim 9.5$, enabling robust measurements of this important transition (Heintz et al., 2023d; Hsiao et al., 2023). Unfortunately, the [O iii] $\lambda\lambda 4960,5008$ doublet is redshifted out of the NIRSpec Prism coverage at approximately the same redshifts, hindering direct $T_{e}$-based metallicity measurements with the JWST/NIRSpec Prism configuration, except for a narrow redshift range ($z\approx 9.3-9.5$).

To 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 set of redshifted Gaussian line profiles at the rest-frame wavelength of each transition and fit for the redshift $z_{\rm spec}$, line equivalent widths (EWs) and fluxes. The line widths (full-width-at-half-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_{\rm spec}$ and the FWHM corrected by the spectral resolution, effectively 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\sigma$ upper limits, are summarized in Table 2 for the JWST-PRIMAL sample and provided in full online⁸⁸8https://github.com/keheintz/jwst-primal, which also includes an expanded line list.

A simple diagnostic of the specific star-formation rate (sSFR) of galaxies is the emission line [O iii]+H$\beta$ 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$\beta$ 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$\beta$ and [O iii] $\lambda\lambda 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.

In Fig. 6, we show the [O iii]+H$\beta$ EW distribution of the full JWST-PRIMAL sample and divided into redshift bins: $z=5.5$–6.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\,\AA$, with 16th and 84th distribution percentiles from $356\,\AA$ to $2543\,\AA$. We observe no apparent evolution 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$\beta$ EWs are much more accurate than photometric estimates, however, the latter potentially being underestimated by up to $\approx 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.

The UV spectral slope, $\beta_{\rm UV}$, where $F_{\lambda}\propto\lambda^{\beta_{\rm UV}}$, of the stellar continuum at wavelengths $\lambda_{\rm rest}=1250$–2600 $\AA$ 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 $\beta_{\rm UV}\approx-2.5$ (Cullen et al., 2017; Reddy et al., 2018a). Steeper spectral power-law indices of $\beta_{\rm UV}\approx-3.0$ would require a very recent onset of star formation ($\lesssim 2\,$Myr) or near-pristine, neutral gas and negligible nebular processing. JWST has enabled measurements of $\beta_{\rm 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$\alpha$ or metal emission (or absorption) line strengths or even low-luminosity AGN.

Here, we measure $\beta_{\rm 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$\alpha$ 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 $\beta_{\rm UV}$ and the normalization constant derived at UV $1550\,\AA$ rest-frame, $m_{\rm UV}$ of the best-fit model, $F_{\lambda}=F_{\rm UV,1550}\lambda^{\beta}$. An example of one of the model fits is shown in Fig. 7. For our full sample, we recover UV continuum slopes ranging from $\beta_{\rm UV}=-1$ down to a few cases with $\beta_{\rm UV}\lesssim-3.0$. This suggests that this set of rare galaxies are particularly dust- and metal-poor ($Z_{\rm gas}/Z_{\odot}<1\%$), have mean stellar populations with ages $\lesssim 2$ Myr (e.g., Cullen et al., 2023a), and likely have ionizing photon escape fractions of $\gtrsim 10\%$ (Chisholm et al., 2022). In Fig. 8, we show the derived $\beta_{\rm UV}$ spectral slopes as a function of redshift. We observe a marginal tendency for galaxies to have steeper UV slopes at 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 iiregions 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).

To determine the UV luminosity, $L_{\rm UV}$, and absolute magnitude, $M_{\rm UV}$, for each source we derive the flux density and corresponding UV magnitude, $m_{\rm UV}$ at rest-frame $1550\,\AA$ from the spectral power-law model. We convert the apparent $m_{\rm UV}$ to the intrinsic, absolute brightness, $M_{\rm UV}$, via

where $\mu(z)$ is the distance modulus at a given redshift derived from the astropy.cosmology module. The derived $M_{\rm 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_{\rm UV}=-22$ to $-16$ mag, respectively (see Fig. 1).

The Balmer recombination lines originating from star-forming H ii regions, H$\alpha$ in particular, are robust tracers of star formation on short ($\lesssim 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$\alpha$ 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$\beta$ for the majority of sources at $z<10$ as:

following Kennicutt (1998), but here assuming the more top-heavy 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$\beta$ line luminosity via $L_{\rm H\beta}=F_{\rm H\beta}\times 4\pi D^{2}_{L}(z)$ with $D_{L}(z)$ the luminosity distance at the given redshift and assume $f_{\rm H\alpha/H\beta}=2.86$ from the Case B recombination scenario at $T_{e}=10^{4}$ K (Osterbrock & Ferland, 2006). We propagate the uncertainty of $\approx 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 $\lesssim 5\%$.

For the sources at $z>10$, where H$\beta$ is redshifted out of the NIRSpec Prism spectral coverage, we instead rely on the [O ii] $\lambda 3727$ line strength

again following Kennicutt (1998), and assuming the Kroupa (2001) IMF. For the full sample, we derive SFRs in the range 0.1–$100\,M_{\odot}\,{\rm 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