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

The escape fraction of Lyman-continuum (LyC) photons ( f esc ) is a key parameter for determining the sources of cosmic reionization at z ≥ 6. At these redshifts, owing to the opacity of the intergalactic medium, the LyC emission cannot be measured directly. However, LyC leakers during the epoch of reionization could be identiﬁed using indirect indicators that have been extensively tested at low and intermediate redshifts. These include a high [O iii ] / [O ii ] ﬂux ratio, high star-formation surface density, and compact sizes. In this work, we present observations of 29 4 . 5 ≤ z ≤ 8 gravitationally lensed galaxies in the Abell2744 cluster ﬁeld. From a combined analysis of JWST-NIRSpec and NIRCam data, we accurately derived their physical and spectroscopic properties: our galaxies have low masses (log( M (cid:63) ) ∼ 8 . 5), blue UV spectral slopes ( β ∼ − 2 . 1), compact sizes ( r e ∼ 0 . 3 − 0 . 5kpc), and high [O iii ] / [O ii ] ﬂux ratios. We conﬁrm that these properties are similar to those characterizing low-redshift LyC leakers. Indirectly inferring the fraction of escaping ionizing photons, we ﬁnd that more than 80% of our galaxies have predicted f esc values larger than 0.05, indicating that they would be considered leakers. The average predicted f esc value of our sample is 0.12, suggesting that similar galaxies at z ≥ 6 have provided a substantial contribution to cosmic reionization.


Introduction
The Lyman-continuum (LyC, λ < 912 Å) photons escaping from star-forming galaxies into the neutral intergalactic medium (IGM) can account for the photon budget required to complete reionization only if a substantial fraction of them escape from the galaxies' interstellar and circumgalactic media (ISM and CGM).Given the number density of star-forming galaxies in the Epoch of Reionization (EoR), an average LyC escape fraction ( f esc ) of ∼ 10% across all galaxies would be required (e.g., Finkelstein et al. 2019;Robertson et al. 2015) to both reionize the Universe by z = 6 and match the Thomson optical depth of electron scattering observed in the cosmic microwave background (CMB, Planck Collaboration et al. 2020).
However, at z ≥ 4.5 it is impossible to detect the LyC photons escaping from galaxies, since they are absorbed and scattered by the IGM along the line of sight (Inoue et al. 2014).Therefore, efforts at low redshift, where LyC can be detected, have been focused on identifying other observable properties that trace physical conditions facilitating the escape of LyC photons.These indirect indicators could then be used in the EoR to identify the cosmic ionizers.Several indirect diagnostics have been proposed (e.g., Yamanaka et al. 2020;Izotov et al. 2018b;Marchi et al. 2018;Verhamme et al. 2017), but they are all characterized by a large scatter.One of the best indicators is the presence of a strong Lyα emission (e.g., Pahl et al. 2021;Gazagnes et al. 2020), often characterized by two emission peaks with a small velocity separation.However, at z > 6.5, this line is attenuated due to its resonant nature and by the increasingly neutral IGM as we approach the EoR (Pentericci et al. 2018;Mason et al. 2019;Jung et al. 2020;Ouchi et al. 2020;Bolan et al. 2021).Emission from the [Mgii]λλ2796, 2803 doublet has been proposed as a promising LyC proxy (Chisholm et al. 2020), as the escape of this line is controlled by resonant scattering in the same low column-density gas as the Lyα (see also Xu et al. 2022;Izotov et al. 2022).More recently, the nebular C iv emission line, requiring comparably high ionization energies to the He ii line (E > 47.9 eV and > 54.4 eV respectively), has attracted attention since its presence might be strongly linked to the escape of Lyman continuum photons from galaxies (e.g., Schaerer et al. 2022;Senchyna et al. 2022;Saxena et al. 2022;Mascia et al. 2023), although this line is in general much fainter than Lyα.Another very popular indicator is the emission line ratio [O iii]λλ4959, 5007/[Oii]λ3727 (hereafter, O32), as Nakajima & Ouchi (2014) first found evidence for its correlation with f esc : high values of O32 would reflect partially incomplete Hii regions, where some LyC photons could escape from (Marchi et al. 2018).Later on, it was found that this correlation is characterized by a substantial scatter, as highlighted by low-redshift studies (e.g., Izotov et al. 2018b;Nakajima et al. 2020;Flury et al. 2022).As a result, a high O32 flux ratio is still a necessary condition for a significant measurement of f esc , although it is not sufficient by itself to define a LyC leaker, as viewing angles might play a role as well as variation in metallicity and ionization parameter (e.g., Bassett et al. 2019;Katz et al. 2020).Further additional properties, such as low values of Balmer lines' rest-frame equivalent widths (EW 0 ), are thus required.As a matter of fact, measuring low values of EW in these lines could indicate lower optical depth in the Hii region (e.g., Bergvall et al. 2013).However, many local LyC emitters exhibiting high values of Balmer emission line EW 0 s are present in the literature (Izotov et al. 2016a(Izotov et al. ,b, 2018a,b),b).A possible explanation for this discrepancy lies in the fact that Balmerline EW 0 s are sensitive not only to Hii region size, but also to starburst age and star formation history (Zackrisson et al. 2017;Binggeli et al. 2018;Alavi et al. 2020).It is thus imperative to pair the EW diagnostic with another indicator in order to prevent degeneracy.Finally, another diagnostic for a high f esc is the SFR surface density (Σ S FR ): feedback from star formation can blow bubbles or chimneys into the host galaxy's ISM, linking a high Σ S FR value to a high f esc .This connection is supported by the detection of some compact LyC emitters (LCEs, e.g., Izotov et al. 2018b;Marchi et al. 2018), even though compactness may not be a defining characteristic of all of them (e.g., Marchi et al. 2018).
Recently, Flury et al. (2022) presented the first statistical test of many of the diagnostics just described on a large sample of low-redshift LCEs.Their conclusion is that indicators based on Lyα emission properties (such as peak velocity separation and equivalent width) perform well, and that diagnostics such as the [Oiii]/[Oii] flux ratio and the SFR surface density, Σ S FR , could also be used to determine if a galaxy is an LCE or not.
With the advent of JWST, we now have the opportunity to constrain LyC diagnostics during the epoch of reionization with the first NIRSpec observations of high-redshift galaxies.In this paper, we make use of NIRSpec spectra obtained in the Abell 2744 cluster region for a sample of galaxies at 4.5 ≤ z ≤ 8 observed as part of the JWST Early Release Science program GLASS (Treu et al. 2022) and the JWST Director Discretionary Time program (JWST-GO-2756; P.I.Chen; Roberts-Borsani et al. 2022a).With the wavelength coverage of 0.9 − 5.3 µm offered by the NIRSpec observations, we can confirm not only the redshifts of tens of photometrically selected candidates using intense emission lines, but we can also measure optical line ratios and rest-frame EWs with high precision.Our unique data set also includes deep JWST/NIRCam images, enabling us to better characterize those cosmic reionizer candidates based on their physical properties.
This paper is organized as follows.We present the data set in Sect. 2. We characterize the selected sample and compare the physical and spectroscopic proprieties with cosmological models and other samples at lower redshifts in Sect.3. In Sect.4, we indirectly infer f esc for the high-redshift sample and in Sect.5, we summarize our key conclusions.Throughout this work, we assume a flat ΛCDM cosmology with H 0 = 67.7 km s −1 Mpc −1 and Ω m = 0.307 (Planck Collaboration et al. 2020) and the Chabrier (2003) initial mass function.All magnitudes are expressed in the AB system (Oke & Gunn 1983).

Observations and method
The final target list of the GLASS-JWST program and the way in which targets have been prioritized will be presented in a future paper.However, Morishita et al. (2022) have already described the observations and data reduction strategy.Here, we present a brief summary, highlighting the points that are most relevant for this work and describing the methods used to study the properties of the galaxies in our sample.

JWST/NIRSpec MSA observations and data reduction
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 µm, at R ∼ 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 µm at R ∼ 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)1 for Level 1 data products and the msaexp2 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, flatfielding, 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 threeshutter 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).

Imaging data
Deep NIRCam images were acquired from the GO program UN-COVER (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.

Emission line and redshift identifications
The focus of the study is on all sources at z ≥ 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α 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 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β, [O iii]λλ4959, 5007, and (when present) Hα lines.In 29 cases, the [Oii]λ3727, [O iii]λλ4959, 5007 and Hβ were detected and their line fluxes were measured.We also measured the Hα 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 specutils3 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α emission (flux on the order of 1 − 3 × 10 −19 erg/s/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 F150W 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 F150W magnitudes derived from the imaging data when present.

Dust correction and flux measurements
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β ≤ 2, we set 2σ as an upper limit.Prior to carry out a quantitative analysis, it is necessary to consider corrections for dust reddening.For 22 galaxies, Hα and Hβ 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 Hα/Hβ = 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α is not in the observed range, we applied the average correction E(B − V) neb ∼ 0.1 derived from all other sources.
With the dust corrected values, we calculated the R23 = ([Oiii] + [Oii])/Hβ and O32 line fluxes ratios and the Hβ restframe 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 highz star-forming galaxies (e.g., Curti et al. 2023;Nakajima et al. 2023), however the Hα/Hβ 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.

Measurements of physical parameters
Following Santini et al. (2022), we measured the stellar masses M ,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)∝ (t 2 /τ) • exp(−t/τ) -with τ 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σ uncertainties on the physical parameters by retaining, for each object, the minimum and maximum fitted masses among all the solutions with a probability P(χ 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 µ 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 µ values range from 1.6 to 12, with a median value of ∼ 2. The results on our sources' stellar masses, M ,true , their M 1500 , and their lensing magnification factors µ are reported in Table 2.

Sizes and SFR surface density estimates
The SFRs were derived from the Hα emission line using the standard conversion by Kennicutt (1998), that is, S FR(Hα)[M yr −1 ] = 7.9 10 42 L Hα , and then corrected for the Chabrier IMF.For those few sources at z ≥ 6.6, where the Hα line falls outside the observable window we used dust-corrected Hβ 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α line in our sample.Thanks to the high resolution of the GLASS-JWST spectra, no correction for [Nii] contamination is required for Hα measurements.For galaxies observed by the DDT programs, Hα 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α 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 Σ S FR = S FR/2πr 2 e (e.g., Naidu et al. 2020;Flury et al. 2022).The SFR surface densities Σ S FR are listed in Table 2.

UV-β slopes
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α line and Lybreak, 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 powerlaw of the form f(λ) ∝ λ β .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 β 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 β for each galaxy.
The results on β with associated errors are reported in Table 2.We find a median β of −2.1, with a 1σ 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.

Possible AGN contamination
Our sample was selected solely based on known spectroscopic redshift (via faint and narrow Lyα 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]λ5007/Hβ 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 dia-gram that compares the [Oiii]λ5007/Hβ to the [Nii]/Hα 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α.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]λ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. 2016aIzotov et al. ,b, 2018aIzotov et al. ,b, 2021;;Wang et al. 2019).

LyC indirect diagnostics
As discussed above, Lyα is possibly the best indirect diagnostic of LyC escape, since the conditions that favour the escape of Lyα 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 ≥ 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α 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, β 1200 , r 50 , and Σ S FR as well as EW 0 (Hβ) 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 (Hβ) 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. 2016aIzotov et al. ,b, 2018aIzotov et al. ,b, 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 ≥ 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 β with the Article number, page 7 of 12 A&A proofs: manuscript no.main measurement of the EW of a Balmer line, such as Hβ, and that the predictions are almost independent of the model assumed (i.e., radiation or density bounded nebula).We know that both We show the properties of our sample in Fig. 6.Specifically, at a given value of EW 0 (Hβ), 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 β slopes for our galaxies are perfectly consistent with the average values found for LBGs at z∼ 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 β 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 ≤ f esc ≤ 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 ≤ 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.65kpc), 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 lowz 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 6 analyzed by Shibuya et al. (2015).Specifically, for M UV = −18(−20), the average UV rest-frame sizes are r e 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 (Σ S FR ) for our galaxies span a very wide range, with an average Σ S FR 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 Σ S FR values, as discussed extensively by Naidu et al. (2020), who actually proposed a physically motivated model in which f esc would scale as Σ 0.4 S FR .Flury et al. (2022) identify a threshold value of Σ S FR = 10 M yr −1 kpc −2 above which their LCE fraction changes from 10% to 60%, and where indeed we find half of our high-redshift galaxies.

O32-R23 diagnostics
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. 2016aIzotov et al. ,b, 2018a,b) ,b) 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) 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.

Predicting escape fractions of EoR galaxies
In the sections above, we have shown that our high-redshift galaxies have properties that are largely overlapping with those of low-z LCEs.The next step is to try and give an indirect estimate of the f esc values for our galaxies, so as to understand if typical low-mass galaxies at z 5 − 7 could be really the drivers of reionization.We assume that the mechanisms that drive the escape of LyC photons are the same at all redshifts and depend only on the physical properties of the sources.
We proceeded as follows: we used the spectroscopic and physical properties of the 66 galaxies that are part of the Flury et al. (2022) sample, with the additional 22 LCEs from previous studies (Izotov et al. 2016a(Izotov et al. ,b, 2018a(Izotov et al. ,b, 2021;;Wang et al. 2019) to calibrate an empirical relation with the measured f esc values.Specifically, we focused on the following properties: O32, EW 0 (Hβ), β, r e (in kpc), Σ S FR , and M 1500 .For all 88 galaxies, these parameters are accurately measured and, of course, accurate measurements of f esc are available.To this end, we specifically used the f esc derived by the COS UV spectral fits (see definition in Flury et al. 2022).Recently Chisholm et al. (2022) followed a somewhat similar approach, also using the same set of low-redshift observations, but limited to the UV-β slope and an indirect proxy.They provided a scaling relation between β and f esc , although the relation has appreciable scatter that scales with f esc .
We first ran the Spearman rank correlation and found that O32, r e , EW 0 (Hβ), and β (in that order) are the properties that are best correlated with f esc .We then identified one useful fitting linear relation to obtain an estimate of log 10 ( f esc ) on the basis of the other four measured physical properties.A fully data-driven regression analysis was carried out by performing a regularized minimization of the root-mean-square error (MSE), computed between the values provided by the equation and the dataset, for several possible combinations of the above properties.We tested both a scheme in which the error on the escape fraction, f esc , is not considered and a scheme in which values are weighted by the error.Since O32 and EW 0 (Hβ) exhibit a very tight correlation (Spearman correlation between them > 0.9, see also Fig. 4), the information they provide is redundant and therefore we decided to use only O32.We checked that the remaining three parameters are reasonably independent (Spearman correlation < 0.5) and therefore provide complementary information.We finally found a best-fit relation of the form: log 10 ( f esc ) = A + B log 10 (O32) + Cr e + Dβ. (1) After identifying the best type of equation, we repeated the minimization process 100 times, following a bootstrap approach every time a random number of sources (between 1 and 25, to avoid being left with too few sources) is randomly removed from the sample.In this way, we could constrain the confidence interval for the above best fit parameters.We Testing the relation on the residuals, we find that it tends to slightly overestimates the f esc at very low values (much lower than 0.01), while it tends to underestimate the f esc at values that are higher than ∼0.1-0.2.This is due to the fact that in the sample used to fit the relation, there are very few galaxies with high f esc values and, in general, their measured f esc values have higher errors.We finally applied the above relation to our sample of highredshift galaxies: out of 29 sources, 3 galaxies do not have the r e measurement since they are outside the UNCOVER footprint; for another 2, we do not have an estimate for the O32 parameter.We therefore could apply the best-fitting relation only to 24 sources.In Fig. 8, we show the predicted f esc for our sources as a function of r e , as well as the low-redshift comparison sample (for which the f esc are measured values).Most of our galaxies have predicted f esc values larger than 0.05, meaning that they would be considered leakers.The average f esc is 0.12 with the bluest and most compact sources having f esc as large as 0.2−0.4,which could be lower limits for the reasons discussed above.
Clearly, the main limitations of the above analysis are the fact that the low-z sample used to calibrate the relation is small and, most importantly, it is not evenly populated as it contains mostly objects with rather low f esc .Also, it would be important to test these predictions at intermediate redshift (z ∼ 3 − 4), that is, much closer to the EoR, where it is possible to both directly detect LyC emission and determine all other physical and spectroscopic properties.However, at the moment, measurements of f esc at intermediate redshift are still sparse and most samples lack near-infrared (NIR) spectroscopic follow-up results.In spite of these limitations, a consistent picture emerges from our results, suggesting low-mass galaxies at z ∼ 5 − 7 have mostly properties which indicate moderate values of f esc .Interestingly, the average f esc value inferred for our sample is equal to the one predicted by Naidu et al. (2020) for galaxies in the mass range log 10 M = 8 − 9 at z ∼ 6, according to their simple model, where f esc scales with Σ S FR , which was constrained to reproduce

Summary and conclusions
Thanks to the magnification power of the Abell 2744 cluster, in this paper, we present the first JWST/NIRSpec observations of 29 gravitationally lensed galaxies with photometric or spectroscopic redshift in the range 4.5 ≤ z ≤ 8. From a combined NIRSpec and NIRCam analysis, we were able to derive accurate physical properties of the galaxies, including M , SFR, r e , UV-β slopes, and measurements of the most prominent optical emission lines ([Oii], [Oiii], Hβ, Hα).We summarize our findings as follows: -Our sample is composed of purely star-forming galaxies, as inferred from the MEx diagram that excludes the presence of any AGN.-The galaxies in our sample have blue UV slopes (median β = -2.08)and are mostly very compact with r e 0.2 − 0.5 kpc.These properties are consistent with those of the general population of LBGs at similar redshift and with similarly faint M UV (Bouwens et al. 2014;Shibuya et al. 2015).-Compared to the low-z sample by Flury et al. (2022) as well as to the low-z LCEs from previous studies, our galaxies present properties (in terms of O32, EW 0 (Hβ), UV-β slopes, r e , and Σ S FR ) that are entirely consistent with those of low-z galaxies with measured f esc larger than 0.05 and which are considered to be LyC leakers.-Using a linear analysis that employs the minimization of MSE, we found a best fitting relation between f esc and the three most correlated and independent parameters (O32, UVβ, and r e ) for the low-redshift sample of 88 galaxies.Applying this relation to the 24 galaxies from our sample where these parameters are all measured, we find that 20/24 have predicted escape fractions larger than 0.05, that is, they would be considered leakers and the average f esc of our sample is 0.12.
In conclusion, our results show that indeed the average lowmass galaxies around the epoch of reionization have physical and spectroscopic properties consistent with moderate values of escaping ionizing photons ( f esc = 0.1 − 0.2).With upcoming JWST observations from GTO and GO programs, this analysis could be easily extended to larger samples of galaxies with higher masses and/or at higher redshift to determine if the bulk of the reionizing photons have been provided by even more massive and brighter galaxies (as predicted, e.g., by Sharma et al. 2016;Naidu et al. 2020) or if we need a more substantial contribution by the fainter and more numerous galaxies with M UV > −18, as predicted by most other popular models (e.g., Trebitsch et al. 2022;Finkelstein et al. 2019;Atek et al. 2015).

Fig. 1 :
Fig. 1: Spatial location of the 29 selected sources, color coded by their spectroscopic redshift.They are superimposed on the RGB image of the UNCOVER program, made with NIRCam filters (blue = F115W + F150W, green = F200W + F277W, and red = F356W + F410M + F444W).The MUSE footprint is shown in white, the NIRSpec GLASS-JWST pointing is shown in cyan, and the NIRSpec DDT pointing is shown in red.
Fig.2: Example 2D and 1D spectra of two representative galaxies in our sample.Top: 2D (top) and 1D (bottom) NIRSpec GLASS-JWST spectrum of 110000 at z = 5.765 (green for the G235H/F170LP and purple for the G395H/F290LP configuration) and 1σ uncertainty (gray).Bottom: 2D (top) and 1D (bottom) NIRSpec DDT spectrum of 150015 at z = 5.041 (black) and 1σ uncertainty (gray).Vertical dashed lines mark the position of well-detected rest-optical emission lines.The continuum emission is also detected as seen in the 2D spectrum.On the X-axis in the bottom panel, the observed wavelength (Å) is reported.

Fig. 4 :
Fig. 3: MEx diagram for the sample of galaxies analyzed in this paper:black dots and green squares are for the GLASS-JWST and DDT samples, respectively.For reference, we plot also the galaxies at z = 0.3 − 0.4 from Flury et al. (2022) (diamonds) and the LCE candidates from previous studies (stars, Izotov et al. 2016a,b, 2018a,b, 2021; Wang et al. 2019).The two orange demarcation lines from Juneau et al. (2014) show the boundaries of the AGN-and-star-forming transition region.All objects above the upper line are AGN-dominated; all galaxies below or rightward of the lower line are presumptively dominated by star formation.We also show the separation from Coil et al. (2015) (dotted lines), which is the adaptation of the Juneau et al. model for galaxies and AGNs at z ∼ 2.3 from the MOSDEF survey.
showed that z ∼ 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∼ 3 LyC leakers presented in Nakajima et al. (

Fig. 6 :
Fig. 6: β vs. rest-frame EW 0 (Hβ).Models are from Zackrisson et al. (2013) and simulate the expected trend for galaxies with an exponential declining SFR (Z = 0.02, solid lines) and various values of escape fractions.Symbols are as in Fig. 4.
Fig. 7: O32 vs. diagnostic diagram for our sample.Black dots and green squares show the GLASS-JWST and DDT samples, respectively.For reference, we plot also the SDSS galaxies from Thomas et al. (2013) as dark grey points.The LCEs from Nakajima et al. (2020) are shown red dots.The Ion2 at z = 3.2 comes from Vanzella et al. (2016); de Barros et al. (2016) as a red x sign, the low-z LyC leakers from Izotov et al. (2016a,b, 2018a,b); Wang et al. (2019) as stars, and the LCEs at z = 0.3 − 0.4 from the Lowredshift Lyman Continuum Survey (as pink diamonds, Flury et al. 2022).We show also the locus of high-redshift LCEs predicted from cosmological simulations by Katz et al. (2020) as an orange dashed line and the threshold of O32 > 5 determined in Flury et al. (2022) as a pink dashed line.

Fig. 8 :
Fig. 8: Predicted f esc vs. size r e for our sample (blue and black symbols) and measured f esc vs. size r e from the literature.Symbols are as in Fig. 4.

Table 1 :
Sample and spectroscopic measurements.
a from MUSE catalog, magnitude F140W HST.Missing O32 and R23 values are due to[Oii]undetection because of its position in un-acquired part of the spectrum.Hβ and [Oiii] are always detected.

Table 2 :
Physical parameters derived from multi-band photometry.All measurements are corrected for magnification where needed.
* Median magnification of the lens model byBergamini (in prep.), calculated at the position of the source.Measurements are associated with 1 σ uncertainties.a measured from HST photometry