Closing in on the sources of cosmic reionization: First results from the GLASS-JWST program

Our spectra were acquired through NIRSpec MSA observations in two programs: the GLASS-JWST Early Release Science Program (PID 1324, PI Treu; Treu et al., 2022) and a JWST DDT program (PID 2756, PI. W. Chen; Roberts-Borsani et al., 2022a), which obtained NIRSpec observations for a subset of targets residing over the central regions of the Frontier Field galaxy cluster Abell 2744.

The GLASS-JWST observations were carried out on November 10, 2022, with three spectral configurations (G140H/F100LP, G235H/F170LP, and G395H/F290LP). These configurations cover wavelengths between $1-5.14\ \mu$m, at $R\sim 2000-3000$. We exposed each of the three high-resolution gratings for a total of 4.9 hours. Specifically, in this work, we use the G235H/F170LP and G395H/F290LP observations, which contain the bright emission lines we will analyze.

DDT NIRSpec observations were carried out on October 23 2022, using the CLEAR filter+PRISM configuration, which provides continuous wavelength coverage of $0.6-5.3\ \mu$m at $R\sim 30-300$ spectral resolution. The on-source exposure time is 1.23 hours.

Data were reduced using the official STScI JWST pipeline (ver.1.8.2)¹¹1https://github.com/spacetelescope/jwst for Level 1 data products and the msaexp²²2https://github.com/gbrammer/msaexp code for Level 2 and 3 data products, which is based on the STScI pipeline but also includes additional correction routines. In summary, we initially reduced the uncalibrated data using the Detector1Pipeline routine and the latest set of reference files (jwst_1023.pmap) to correct for detector-level artifacts and convert them to count-rate images. Then, we applied custom preprocessing routines from msaexp to remove residual $1/f$ noise that is not corrected by the IRS2 readout, to identify and remove ”snowballs”, and to remove bias exposure by exposure before running STScI routines from Spec2Pipeline for the final 2D cutout images. To perform WCS registration, flat-fielding, path-loss corrections, and flux calibration, these routines include AssignWcs, Extract2dStep, FlatFieldStep, PathLossStep, and PhotomStep. Of note, our chosen reference files include an in-flight flux calibration, accounting for NIRSpec’s better-than-expected throughput at blue wavelengths. Local background subtraction was performed using a three-shutter nod pattern before the resulting images are drizzled onto a common grid. We optimally extracted the spectra using an inverse-variance weighted kernel, which is derived by summing the 2D spectrum along the dispersion axis and fitting the signal along the spatial axis to a Gaussian profile. We visually inspected all kernels to make sure spurious events are not included. As a result, the kernel extracts the 1D spectrum along the dispersion axis. The final step was to verify the default wavelength calibration for the gratings, which is accurate within 1 Å (Williams et al. in prep).

Deep NIRCam images were acquired from the GO program UNCOVER (GO 2561; PI I. Labbe) and included observations in the F115W, F150W, F200W, F277W, F356W, F410M, and F444W filters. The imaging data were reduced using the STScI JWST pipeline and the latest versions of photometric zero points and reference files. A detailed description of all the reduction and calibration steps is presented in Merlin et al. (2022). Of the 29 sources analyzed in this work (see next section), 27 were observed within the UNCOVER pointings. Their positions and the UNCOVER footprints are presented in Fig. 1.

The focus of the study is on all sources at $z\geq 4.5$. Specifically, we analyzed the spectra of: 13 galaxies with a spectroscopic redshift larger than 4.5 previously confirmed by the MUSE observations (Mahler et al., 2018; Richard et al., 2021), all showing Ly$\alpha$ in emission in their optical spectra (the footprint of the MUSE observations is also shown in Fig. 1); 29 galaxies with a photometric redshift in the range 4.5-8, of which 5 were selected as part of the $z\simeq 7.9$ candidate protocluster and whose confirmation has been recently presented by Morishita et al. (2022). Finally, 23 of the 42 galaxies were observed as part of GLASS-JWST, while 19 were part of the DDT program.

Spectra were visually examined for detectable optical lines using the spectroscopic or photometric redshift information. For photometric sources we determined the spectroscopic redshift when possible using the H$\beta$, [O iii]$\lambda\lambda 4959,5007$, and (when present) H$\alpha$ lines. In 29 cases, the [Oii]$\lambda 3727$, [O iii]$\lambda\lambda 4959,5007$ and H$\beta$ were detected and their line fluxes were measured. We also measured the H$\alpha$ line in 17 out of 29 cases, since it falls within the observed spectrum (examples are shown in Fig. 2). For this part of our analysis, we used the latest version of the specutils³³3https://specutils.readthedocs.io/en/stable/index.html packages in python.

From our initial target list of 42 galaxies, there were eight sources with a spectroscopic redshift confirmed from previous MUSE observations, which were placed in the MSA, but for which we are unable to confirm any emission line from JWST data. Of these, five were observed as part of GLASS-JWST and three during the DDT. Most of these sources are extremely faint (fainter than 28 F150W), thus, their redshift confirmation was solely based on the presence of faint Ly$\alpha$ emission (flux on the order of $1-3\times 10^{-19}\ \text{erg}/\text{s}/\text{cm}^{2}$) in the MUSE cubes. We were also unable to detect any emission lines for five galaxies with a photometric redshift between 4.5 and 7: these galaxies are also faint ($m_{\text{F150W}}\simeq 27.2-28$). However, in these cases it is possible that the photometric redshift was incorrect and this is why we are unable to confirm any emission line in their spectra.

In Table 1, we report the GLASS-JWST IDs, the coordinates and the spectroscopic redshifts of all sources together with their $m_{\text{F150W}}$ magnitudes derived from the imaging data when present.

We measured the total flux of all detected lines (Balmer lines, [Oii], and of [Oiii]) with a Gaussian fit of each line component. From the flux measurement we subtracted a constant continuum emission, which is estimated as the signal averaged over regions in which there are no lines near the one being measured. When the signal-to-noise ratio (S/N) of [Oii] or H$\beta$  $\leq 2$, we set 2$\sigma$ as an upper limit. Prior to carry out a quantitative analysis, it is necessary to consider corrections for dust reddening. For 22 galaxies, H$\alpha$ and H$\beta$ are both available: for these, we calculated the correction for dust extinction on the basis of the Balmer decrement, assuming a Calzetti et al. (2000) extinction law and an intrinsic ratio ${\rm H}\alpha/{\rm H}\beta=2.86$ (see e.g., Domínguez et al., 2013; Kashino et al., 2013; Price et al., 2014), which is valid for an electric temperature of 10000 K. For the other 7 sources where H$\alpha$ is not in the observed range, we applied the average correction $E(B-V)_{neb}\sim 0.1$ derived from all other sources.

With the dust corrected values, we calculated the $R23=([\textrm{O}\textsc{iii}]+[\textrm{O}\textsc{ii}])/{\rm H}\beta$ and $O32$ line fluxes ratios and the H$\beta$ rest-frame $EW_{0}$s. We list all these values in in Table 1. It is important to highlight that the $O32$ values slightly depend on the dust correction. The assumed temperature may be too low for high-z star-forming galaxies (e.g., Curti et al., 2023; Nakajima et al., 2023), however the H$\alpha$/H$\beta$ ratio varies by less than 5% for a 20000 K temperature (Dopita & Sutherland, 2003) and therefore the $O32$ variation are well below the current uncertainties.

Following Santini et al. (2022), we measured the stellar masses $M_{\star,obs}$, the observed absolute UV magnitudes at 1500Å ($M_{1500,obs}$), the star formation rates (SFRs), and dust reddening $E(B-V)$ by fitting synthetic stellar templates to the seven-band NIRCam photometry and the released HST photometry (Castellano et al., 2016) with zphot (Fontana et al., 2000). We adopted Bruzual & Charlot (2003) models and assumed delayed exponentially declining star formation histories – SFH($t$)$\propto(t^{2}/\tau)\cdot\exp(-t/\tau)$ – with $\tau$ ranging from 0.1 to 7 Gyr. The age ranges from 10 Myr to the age of the Universe at each galaxy redshift, while metallicity can assume values of 0.02, 0.2 or 1 times Solar metallicity. For the dust extinction, we used the Calzetti et al. (2000) law with $E(B-V)$ ranging from 0 to 1.1. We computed $1\sigma$ uncertainties on the physical parameters by retaining, for each object, the minimum and maximum fitted masses among all the solutions with a probability $P(\chi^{2})>32\%$ of being correct, fixing the redshift to the best-fit value.

As a final step, we needed to adjust the stellar masses and $M_{1500}$ and the SFR values for the lensing amplification factor. The magnification factor $\mu$ was derived by combining the LM-model (Bergamini et al., 2022), the model used by Roberts-Borsani et al. (2022b), with a new spatially more extensive model that fully covers the JWST field of view (Bergamini et al. in prep). The $\mu$ values range from 1.6 to 12, with a median value of $\sim 2$. The results on our sources’ stellar masses, $M_{\star,true}$, their $M_{1500}$, and their lensing magnification factors $\mu$ are reported in Table 2.

The SFRs were derived from the H$\alpha$ emission line using the standard conversion by Kennicutt (1998), that is, $SFR({\rm H}\alpha)[M_{\odot}\ \text{yr}^{-1}]=7.9\ 10^{42}\ L_{{\rm H}\alpha}$, and then corrected for the Chabrier IMF. For those few sources at $z\geq 6.6,$ where the H$\alpha$ line falls outside the observable window we used dust-corrected H$\beta$ fluxes instead. To correct for slit losses and possible residual uncertainties on flux calibration, we normalized the spectra to the F444W filter by integrating the spectrum under the F444W bandpass, the closest to the H$\alpha$ line in our sample. Thanks to the high resolution of the GLASS-JWST spectra, no correction for [Nii] contamination is required for H$\alpha$ measurements. For galaxies observed by the DDT programs, H$\alpha$ is blended with the [Nii] doublet. However our sources are all very-low-mass galaxies and, based on low-redshift results, we expect contamination to be less than 10% (e.g., Faisst et al., 2018). Therefore, we did not attempt to make any correction on the fluxes.

The results obtained on the SFR were also compared to the values determined by the SED fitting procedure described in Sect. 2.5, finding a reasonable agreement between the two data sets, indicating that there are no systematic issues. Obviously, the SED fitting SFR is heavily dependent on the choice of SFH, and also on the dust correction, while the value derived from the H$\alpha$ luminosity is sensitive to short timescales. That is to say that it only probes the presence of short-lived, massive stars, so their ratio is actually an indication of the burstiness of the galaxy.

We measured the size of each galaxy $r_{e}$ in the F115W band, corresponding to the UV rest-frame of the galaxies. Adopting forward-modeling technique, we assumed that galaxies are well represented by a Sersic profile (Sersic, 1968) as Yang et al. (2022). Then we modeled the appearance of the source in the image plane considering lensing and PSF effects. The source reconstruction was performed via python software Lenstruction (see details in Yang et al., 2020), which is built on Lenstronomy (Birrer et al., 2021). In this way we obtained the intrinsic properties of galaxies in the source plane, hence, the intrinsic size. Finally, we calculated the SFR surface density using the relation $\Sigma_{SFR}=SFR/2\pi r_{e}^{2}$ (e.g., Naidu et al., 2020; Flury et al., 2022). The SFR surface densities $\Sigma_{SFR}$ are listed in Table 2.

We measured the UV slope of our galaxies from the NIRCam photometry and/or the previously available HST photometry (Castellano et al., 2016), with the approach detailed in Calabrò et al. (2021). In brief, we considered all the photometric bands whose entire bandwidths are between a $1216$ and $3000$ Å rest frame. The former limit is set to exclude the Ly$\alpha$ line and Ly-break, while the latter limit is slightly larger than that adopted in Calabrò et al. (2021) to ensure a larger range.

Then, we fitted the selected photometry with a single power-law of the form f($\lambda$) $\propto\lambda^{\beta}$. In practice, we fitted two or three photometric bands amongst HST F814W and JWST-NIRCam F115W, F150W or F200W depending on the exact redshift of the sources. This choice allows us to uniformly probe the spectral range between $1500$ and $3000$ Å for most of the galaxies. We measured the $\beta$ and associated uncertainty for each source using a bootstrap method. By using $n=500$ Monte Carlo simulations, the fluxes in each band were varied according to their error. The results provided a resultant slope distribution from which we calculated the mean and standard deviation of $\beta$ for each galaxy. The results on $\beta$ with associated errors are reported in Table 2. We find a median $\beta$ of $-2.1$, with a $1\sigma$ dispersion of $0.4$. The UV magnitudes $M_{1600}$ that can be derived simultaneously with this approach are consistent with the values obtained from the SED fitting and described in the previous section.

Our sample was selected solely based on known spectroscopic redshift (via faint and narrow Ly$\alpha$ emission) or through photometric redshifts. It could therefore be affected by AGN contamination. Since we are interested in searching for candidate LyC emitters amongst the star-forming population, we first consider whether the primary source of ionization might not be star formation, but AGN activity instead.

We employed the mass-excitation (MEx) diagram, first proposed by Juneau et al. (2011) to combine the measurements of the [Oiii]$\lambda 5007$/H$\beta$ emission line ratio with the stellar mass to discriminate between AGNs and star-forming galaxies. This diagram was proposed as an alternative to the classical BPT diagram that compares the [Oiii]$\lambda 5007$/H$\beta$ to the [Nii]/H$\alpha$ emission line ratios (Baldwin et al., 1981) when these last two lines fall out of the visibility window and it is not possible to use them to characterize the ionization mechanism in the galaxies – as is the case for our low-resolution DDT spectra where the [Nii] doublet is blended with H$\alpha$. Due to the high resolution of the GLASS-JWST spectra, for the galaxies with these observations, we used the dust-corrected flux measurement of the 5007 Å component of the [Oiii] doublet. Instead, for the galaxies observed as part of the DDT program for which the doublet cannot be resolved, we determined the [Oiii]$\lambda 5007$ flux, assuming the expected line ratio of 1:3 for the two components fixed by atomic physics. As shown in Fig. 3, the position of our sources in the MEx diagram indicates that our sample contains essentially star forming galaxies, lying below or around the division line identified by Coil et al. (2015) for $z=2.3$ galaxies and AGN from the MOSDEF survey. For reference, we also plot the galaxies at $z=0.3-0.4$ from the Low-redshift Lyman Continuum Survey from Flury et al. (2022) and the compilation of low-z LCE also used by the same authors (from Izotov et al., 2016a, b, 2018a, 2018b, 2021; Wang et al., 2019).

As discussed above, Ly$\alpha$ is possibly the best indirect diagnostic of LyC escape, since the conditions that favour the escape of Ly$\alpha$ photons are often the same that allow for the escape of LyC photons. However, we cannot use this parameter for our sample, since the NIRSpec data do not cover the 1216Å region for most of our galaxies and only a small subset of our sources have been covered by previous MUSE observations (see also Fig. 1). In addition, for the sources at $z\geq 7$, the IGM becomes significantly neutral, thus absorbing the emission even in sources where the line would be intrinsically bright (e.g., Stark et al., 2010; Pentericci et al., 2011; Mason et al., 2018). Indeed, Morishita et al. (2022) recently showed that none of the galaxies in the protocluster candidate at $z=7.89$ presents bright Ly$\alpha$ emission and estimate an average neutral hydrogen fraction of the IGM in the region to be $>0.45$. Therefore in this work we concentrate on the other most promising diagnostics tested at low redshift by Flury et al. (2022) and available for our galaxies: these authors showed that $O32$, $\beta_{1200}$, $r_{50}$, and $\Sigma_{SFR}$ as well as $EW_{0}({\rm H}\beta)$ and $M_{1500}$ exhibit some of the strongest and most significant correlations with $f_{esc}$. This would indicate that certain characteristics, such as concentrated star formation, young stellar populations, and high ionization states, play an essential role in the escape of LyC photons.

We began by analyzing the $O32$ and the rest-frame $EW_{0}({\rm H}\beta)$ relation, as shown in Fig. 4. We plot the values for the JWST high-redshift sample and compare them to the local galaxies with measured $f_{esc}$ from previous works (Flury et al., 2022; Izotov et al., 2016a, b, 2018a, 2018b, 2021; Malkan & Malkan, 2021; Wang et al., 2019). These values are color-coded as a function of their $f_{esc}$ (COS UV) measured values. We can see that the large majority of the high-redshift sources have values of $O32$ larger than 5, which has been indicated as a threshold for LyC leakers (LCEs from here onwards) by Flury et al. (2022). We note that Flury et al. (2022) define leakers as galaxies with an $f_{esc}>0.05$ measured with $S/N\geq 5$. Our sample indeed mostly lies in the region of the plot populated by low-redshift leakers. In Fig. 5, we show $O32$ as a function of total stellar mass, again comparing our sample to the low-z galaxies: our sample perfectly overlaps with the low-z galaxies at the low-mass end, where we find most of the LCEs.

Zackrisson et al. (2013) suggested that galaxies ought to be identified with high $f_{esc}$ by combining the UV slope $\beta$ with the measurement of the EW of a Balmer line, such as H$\beta$, and that the predictions are almost independent of the model assumed (i.e., radiation or density bounded nebula). We know that both $\beta$ and $EW_{0}({\rm H}\beta)$ are also good indirect indicators according to Flury et al. (2022), although there does not seem to be a direct correlation between the two values: Flury et al. (2022) showed that their LCEs do not seem to follow the predictions provided by the models from Zackrisson et al. (2013).

We show the properties of our sample in Fig. 6. Specifically, at a given value of $EW_{0}({\rm H}\beta)$, the high-redshift galaxies are, on average, bluer than the majority of the low-z galaxies (consistent with the lower stellar mass probed by our sources); but their slopes are similar to the subset of the low-z sources with moderate-to-high values for $f_{esc}$. The $\beta$ slopes for our galaxies are perfectly consistent with the average values found for LBGs at z$\sim 6$ for galaxies with $M_{UV}$ in the same range Bouwens et al. (2014). In the figure, we also show the models from Zackrisson et al. (2013) for various values of escape fraction, after correcting the intrinsic $\beta$ slopes of the models for dust attenuation and assuming the average reddening of our sample derived from the SED fitting, $E(B-V)=0.097$ (solid line), and also the maximum and minimum $E(B-V)$ in our sample (shaded area). We note that our galaxies show more consistency with model predictions that assume moderate-to-low escape fractions $0.0\leq f_{esc}\leq 0.5$.

Finally, we analyzed the sizes and star formation rate density for our galaxies. Several authors have postulated that concentrated star formation provides the feedback necessary to clear paths in the ISM that allow for the escape of LyC photons. In addition, Marchi et al. (2018) recently found that stacks of galaxies that are UV-compact ($r_{UV}\leq 0.30$ kpc) have much higher LyC flux than the average population (see also Izotov et al., 2018b). We find that, with the exception of one galaxy (ID 160281, which has $r_{e}=1.65\text{kpc}$), our targets are all extremely compact with typical $r_{e}$ in the rest-frame UV around $0.2-0.6$ kpc. These values are similar or even lower than those found for the low-z galaxies and LCEs (Flury et al., 2022); thus, once again the high-redshift sources have equal properties to the low-z LCEs. Since we know that sizes evolve with redshift and galaxies become progressively more compact, we also compared our sample to the general population of star-forming galaxies (LBGs selected) at $z\simeq 6$ analyzed by Shibuya et al. (2015). Specifically, for $M_{UV}=-18(-20)$, the average UV rest-frame sizes are $r_{e}\simeq 0.38(0.56)$ kpc, respectively. Thus, indeed our galaxies are, from this point of view, as compact as the general UV faint population at the same redshift.

As for the star formation rate densities, the values of $\log_{10}(\Sigma_{SFR})$ for our galaxies span a very wide range, with an average $\Sigma_{SFR}$ that is slightly higher than the average expected at their redshift, as measured by Shibuya et al. (2015). We note that their values are inferred from the UV and then dust corrected. LyC leakers at low and intermediate redshift, tend to show high $\Sigma_{SFR}$ values, as discussed extensively by Naidu et al. (2020), who actually proposed a physically motivated model in which $f_{esc}$ would scale as $\Sigma_{SFR}^{0.4}$. Flury et al. (2022) identify a threshold value of $\Sigma_{SFR}=10\ M_{\odot}\ \text{yr}^{-1}\ \text{kpc}^{-2}$ above which their LCE fraction changes from 10% to 60%, and where indeed we find half of our high-redshift galaxies.

The $O32$ vs $R23$ index diagram (Fig. 7) is widely used to examine the gas-phase metallicity and ionization state both in the local universe (e.g., Thomas et al., 2013; Izotov et al., 2016a, b, 2018a, 2018b) and at high redshift (e.g., Flury et al., 2022; Nakajima et al., 2020; Reddy et al., 2022; Vanzella et al., 2019), as both these indices are sensitive to combination of these quantities. Nakajima et al. (2020) showed that $z\sim 3$ LyC leakers tend to populate the upper right part of this diagram. Recently Katz et al. (2020) used high-resolution cosmological radiation hydrodynamics simulations to examine the properties of LyC leakers deep into the epoch of reionization and found that simulated high-redshift galaxies populate the same regions of the $R23-O32$ plane, as the z$\sim$ 3 LyC leakers presented in Nakajima et al. (2020) that tend to have low metallicity. Although they conclude that this plane is not the most useful to differentiate between leakers and non-leakers we note that the $z=3$ leakers by Nakajima et al. (2020), with measured values and Ion2 occupy the same region as the $z=0.3$ low-redshift leakers. Our galaxies occupy a much broader region, which could reflect either a wider range of metallicity and ionization states or simply the fact that we have large measured uncertainties on the diagnostics.