Early Results from GLASS-JWST. XIV. A Spectroscopically Confirmed Protocluster 650 Million Years after the Big Bang

We base our primary analysis on data acquired through NIRSpec micro‐shutter assembly (MSA) observations in two programs, the GLASS-JWST Early Release Science Program (PID 1324, PI Treu; Treu et al. 2022) and a JWST Director's Discretionary Time (DDT) program (PID 2756, PI. W. Chen; Roberts-Borsani et al. 2022b). The GLASS-JWST observations were executed on 2022 November 10 with three spectral resolution configurations with three high-resolution gratings, G140H/F100LP, G235H/F170LP, and G395H/F290LP, which also provide total wavelength coverage of 0.81–5.14 μm, at R ∼ 1000–3000. The on-source exposure time was 4.9 hr in each spectral configuration. The DDT NIRSpec observations were executed on 2022 October 23, with 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 was 1.23 hr.

For the MSA target selection, we started with the same source catalog for both programs. Specifically for the z ∼ 8 protocluster sources, z/Y-dropout galaxies (hereafter ZDs and YDs, respectively) were included (4 in the DDT and 4 in GLASS-JWST), all within the vicinity of the overdensity (Zheng et al. 2014; Ishigaki et al. 2016) including the spectroscopically confirmed galaxies YD4, GLASSZ8-1 (ZD2), and GLASSZ8-2 from Laporte et al. (2017), Roberts-Borsani et al. (2022a), respectively. Considering the overlap between the two programs, a total of seven distinct protocluster targets were observed, but data was corrupted for one target due a nonnominal operation of a microshutter, leaving seven targets suitable for analysis. The choice of protocluster targets in each MSA was based on three primary factors, namely (i) the central pointing of the MSA, (ii) the position of the MSA ensuring no spectral overlap in the detector, and (iii) the preferential selection of brighter objects to maximize the probability of emission line or continuum detections.

The data were reduced using the official STScI JWST pipeline (ver.1.8.2) ³³ for Level 1 data products, and the msaexp ³⁴ code for levels (2) and (3) data products, the latter of which is built on the STScI pipeline routines but also includes custom routines for additional corrections. Briefly, we begin our data reduction on the uncalibrated files with the Detector1Pipeline routine and the latest set of reference files (jwst_1014.pmap) to correct for detector-level artifacts and convert to count-rate images. We then utilize custom preprocessing routines from msaexp to correct for 1/f noise, identify and remove snowballs, and remove bias on an exposure-by-exposure basis, before running a number of STScI routines from the Spec2Pipeline to produce the final 2D cutout images. These include the AssignWcs, Extract2dStep, FlatFieldStep, PathLossStep, and PhotomStep routines to perform World Coordinate System registration, flat-fielding, pathloss corrections, and flux calibration. Background subtraction is performed locally using a three-shutter nod pattern before drizzling the resulting images onto a common grid. From there, we optimally extract the spectra via an inverse-variance-weighted kernel, derived by summing the 2D spectrum along the dispersion axis and fitting the resulting signal along the spatial axis with a Gaussian profile by following the recipe of Horne (1986). We visually inspect all kernels to ensure spurious events are not included (or limited) where possible. The kernel then extracts the 1D spectrum along the dispersion axis.

Deep NIRCam images are available from DDT program (PID 2756; PI. W. Chen) and General Observers program UNCOVER (Bezanson et al. 2022), including F115W, F150W, F200W, F277W, F356W, F410M, and F444W filters. The imaging data are reduced in the same way as presented by Merlin et al. (2022), using the official STScI JWST pipeline, including the most recent version of the photometric zero-points and reference files. Images are PSF-matched to the F444W filter for the flux estimates below.

To supplement our photometric wavelength coverage, we include ancillary HST data taken by several programs (Postman et al. 2012; Treu et al. 2015; Lotz et al. 2017; Steinhardt et al. 2020). The HST data have been uniformly rereduced using Grizli (Brammer et al. 2022). The HST images are PSF-matched to the F160W filter instead of the NIRCam F444W. The choice was made because—despite their similar PSF FWHMs—there are significant differences in the PSF profile of the two different telescopes, which make it challenging to obtain a satisfying convolution kernel. The remaining systematic offset in fluxes caused by this is corrected in the following process.

A photometric catalog is constructed following Morishita & Stiavelli (2022), using borgpipe (Morishita et al. 2021). Briefly, fluxes are estimated in the PSF-matched images with a r = 032 aperture by using SExtractor (Bertin & Arnouts 1996). The flux offsets between NIRCam and HST filters are corrected with a rescaling based on the mean offset between NIRCam F150W and a pseudo F150W filter derived for the same sources using the HST F140W and F160W fluxes, whose broadband filters straddle the NIRCam F150W. The correcting factor is found to be 1.25, which is consistent with Morishita & Stiavelli (2022).

Lastly, the fluxes are scaled to a total flux by applying C = f_auto,F444W/f_aper,F444W, where f_auto,F444W is FLUX_AUTO of SExtractor, measured for individual sources.

We present our spectroscopic analyses of the seven galaxies in our sample in Figure 2, which shows the 2D spectra and 1D extraction. Remarkably, all galaxies show clear [O iii] 5007 lines at ∼4.4 μm, and tentative Hβ, [O ii], and [Ne iii] lines in a few galaxies (see Mascia et al. 2023 and G. Roberts-Borsani et al. 2023, in preparation. for a dedicated analysis of the emission lines). Here we focus on redshift determinations using the [O iii] 4959,5007-doublet and Hβ line.

Zoom In Zoom Out Reset image size

The redshift of each source is determined by fitting a three-component Gaussian to Hβ and the [O iii]-doublet after subtracting the underlying continuum. Before subtraction, we first scale the observed spectrum by matching the continuum level to the best-fit model, for slitloss and any remaining offset in absolute flux calibration. We use the wavelength range at 4–4.8 μm while emission lines (i.e., Hβ and [O iii]) are masked. For continuum subtraction, we use a best-fit spectral energy distribution (SED) template derived with broadband photometry (Section 3.2). We fix the line ratio of the [O iii]-doublet lines to 1:3 and set the width to a single parameter for the two components of the doublet. For Hβ, the amplitude and line width are set as free parameters. Including the redshift, we have five free parameters. The uncertainties and parameter posterior distribution functions are estimated by sampling the parameter space via emcee (Foreman-Mackey et al. 2013).

In order to assess the detection significance of Hβ and [O iii] emission, we estimate the noise level in the spectrum from 3.6 to 4.8 μm. We measure total fluxes at various wavelengths integrated over ±2σ (where σ is the best-fit Gaussian width of each emission line). We then compare the standard deviation of these fluxes (i.e., noise) with the emission lines integrated over ±2σ around the central wavelength. This test indicates secure (>5σ confidence) detections of [O iii]_λ5007 in all seven photometric candidates that were targeted, along with [O iii]_λ4959 in five and Hβ in three. The resulting line fluxes and spectroscopic redshifts are presented in Table 1. The total line fluxes are measured by integrating the best-fit Gaussian component for each line when detected at >5σ. It is noted that the slit loss is corrected by multiplying a median ratio of the best-fit SED (Section 3.2) and the observed spectrum at 4.0–4.8 μm after masking the wavelengths of Hβ and [O iii]-doublet lines. The line fluxes and spectroscopic redshifts are presented in Table 1.

Table 1. Spectroscopically Confirmed Protocluster Member Galaxies behind Abell 2744

ID	R.A.	Decl.	mF150W	mF444W	Redshift	fHβ	f[OIII]4959+5007	EW0(Lyα) a
	degree	degree	mag	mag		10−19 erg s−1 cm−2	10−18 erg s−1 cm−2	Å
YD4	3.6038544	−30.3822365	26.5	25.6	${7.8758}_{-0.0001}^{+0.0001}$	${8.6}_{-1.4}^{+1.6}$	${20.2}_{-0.7}^{+0.7}$	<15.9
YD7	3.6033909	−30.3822289	26.4	25.8	${7.8772}_{-0.0041}^{+0.0041}$	<76.3	${40.4}_{-8.1}^{+8.8}$	<23.9
ZD6	3.6065702	−30.3808918	26.9	26.1	${7.8843}_{-0.0013}^{+0.0013}$	<13.6	${44.0}_{-3.4}^{+3.3}$	<26.2
YD8	3.5960841	−30.3858051	26.7	26.0	${7.8869}_{-0.0004}^{+0.0004}$	${109.0}_{-13.2}^{+14.6}$	${275.1}_{-6.1}^{+6.2}$	<16.9
ZD2	3.6045184	−30.3804321	26.2	25.2	${7.8800}_{-0.0001}^{+0.0001}$	${19.6}_{-0.8}^{+0.8}$	${70.1}_{-0.5}^{+0.5}$	⋯
ZD3	3.6064637	−30.3809624	26.7	26.6	${7.8816}_{-0.0001}^{+0.0001}$	${4.5}_{-0.7}^{+0.8}$	${6.4}_{-0.2}^{+0.2}$	⋯
GLASSz8-2	3.6013467	−30.3791904	26.5	25.4	${7.8831}_{-0.0001}^{+0.0001}$	${8.2}_{-0.6}^{+0.7}$	${16.1}_{-0.4}^{+0.4}$	<27.6

Note.

^a 2σ rest-frame equivalent width of Lyα over Δ ∼ 100 Å (∼2700 km s⁻¹). Lyα equivalent width measurements of ZD2 and ZD3 are not available as the Lyα wavelength falls in the detector gap. Measurements here are not corrected for magnification. For those with Hβ not detected, 2σ flux limit (assuming 100 Å for the line width) are presented.

Download table as:  ASCIITypeset image

For the four galaxies observed with the prism configuration (YD4, YD7, YD8, and ZD6) and one with the high-resolution grating (GLASSz8-2), we have spectroscopic coverage at the wavelength of Lyα, allowing us an independent check of the inferred redshift measurements. Indeed, we confidently detect the Lyα break at the expected wavelength for the redshift derived above for all of the four galaxies (Figure 2). For the other two galaxies (ZD2 and ZD3), while a small part of the wavelength range near Lyα falls in the detector gap, the break is still consistent with the inferred redshift.

We analyze the SED of the individual galaxies by using photometric data that covers 0.4–5 μm. We use the SED fitting code gsf (Morishita et al. 2019), which allows flexible determinations of star formation histories in a nonparametric form, by finding an optimal combination of stellar and interstellar medium (ISM) templates. We generate templates of different ages, [1, 3, 10, 30, 100, 300] Myr, and metallicities $\mathrm{log}{Z}_{* }/{Z}_{\odot }\in [-2:0.4]$ at an increment of 0.1, by using fsps (Conroy et al. 2009; Foreman-Mackey et al. 2014). A nebular component (emission lines and continuum) that is characterized by an ionization parameter $\mathrm{log}U\in [-3:-1]$ is also generated by fsps (see also Byler et al. 2017) and added to the template after multiplication by an amplitude parameter. Dust attenuation and metallicity of the stellar templates are treated as free parameters during the fit, whereas the metallicity of the nebular component is fixed to the same value of the stellar component during the fitting process.

The posterior distribution function of the parameters is sampled by using emcee for 10⁵ iterations with the number of walkers set to 100. The final posterior is collected after excluding the first half of the realizations (known as burn-in). The resulting physical parameters are quoted as the median of the posterior distribution, along with the 16th–84th percentile uncertainty ranges. The star formation rate is calculated by averaging the last 100 Myr of the posterior star formation history. The inferred physical properties are presented in Table 2.

Table 2. Physical Properties of the Final Candidates

ID	μa	MUV	βλ	$\mathrm{log}{M}_{* }$	SFR	$\mathrm{log}{t}_{* }$	$\mathrm{log}{Z}_{* }$	AV	RH II b	ξion(1 − fesc) c
		mag		M⊙	M⊙ yr−1	Gyr	Z⊙	mag	pMpc	1025 Hz erg−1
Spectroscopic Sample
YD4	${2.01}_{-0.04}^{+0.05}$	$-{19.69}_{-0.06}^{+0.08}$	$-{1.28}_{-0.11}^{+0.07}$	${9.17}_{-0.17}^{+0.27}$	${6.46}_{-1.99}^{+1.94}$	$-{1.44}_{-0.53}^{+0.69}$	$-{1.30}_{-0.38}^{+0.42}$	${1.10}_{-0.15}^{+0.16}$	${0.10}_{-0.09}^{+0.31}$	${4.06}_{-1.36}^{+1.49}$
YD7	${2.02}_{-0.04}^{+0.05}$	$-{19.96}_{-0.05}^{+0.05}$	$-{1.86}_{-0.04}^{+0.02}$	${9.39}_{-0.15}^{+0.08}$	${4.04}_{-1.57}^{+1.51}$	$-{0.71}_{-0.26}^{+0.22}$	$-{1.16}_{-0.49}^{+0.37}$	${0.25}_{-0.14}^{+0.16}$	${0.76}_{-0.17}^{+0.11}$	⋯
ZD6	${1.94}_{-0.04}^{+0.05}$	$-{19.39}_{-0.10}^{+0.09}$	$-{1.53}_{-0.07}^{+0.08}$	${8.86}_{-0.19}^{+0.27}$	${3.34}_{-1.08}^{+2.98}$	$-{1.41}_{-0.38}^{+0.40}$	$-{0.86}_{-0.69}^{+0.76}$	${0.71}_{-0.26}^{+0.23}$	${0.26}_{-0.17}^{+0.18}$	⋯
YD8	${2.63}_{-0.06}^{+0.07}$	$-{19.28}_{-0.05}^{+0.04}$	$-{2.25}_{-0.05}^{+0.02}$	${8.36}_{-0.26}^{+0.09}$	${1.34}_{-0.66}^{+0.30}$	$-{1.48}_{-0.27}^{+0.18}$	$-{0.43}_{-0.42}^{+0.13}$	${0.14}_{-0.04}^{+0.06}$	${0.53}_{-0.18}^{+0.07}$	${19.59}_{-2.33}^{+2.59}$
ZD2	${1.98}_{-0.04}^{+0.05}$	$-{20.13}_{-0.02}^{+0.02}$	$-{2.24}_{-0.01}^{+0.01}$	${8.73}_{-0.05}^{+0.03}$	${3.13}_{-0.37}^{+0.24}$	$-{1.49}_{-0.31}^{+0.16}$	$-{0.29}_{-0.05}^{+0.08}$	${0.20}_{-0.01}^{+0.01}$	${0.58}_{-0.02}^{+0.03}$	${10.98}_{-0.61}^{+0.63}$
ZD3	${1.94}_{-0.04}^{+0.05}$	$-{19.67}_{-0.07}^{+0.06}$	$-{2.07}_{-0.03}^{+0.05}$	${8.33}_{-0.15}^{+0.32}$	${1.00}_{-0.28}^{+0.59}$	$-{1.44}_{-0.47}^{+0.61}$	$-{0.64}_{-0.68}^{+0.59}$	${0.22}_{-0.14}^{+0.14}$	${0.07}_{-0.05}^{+0.23}$	${1.62}_{-0.87}^{+0.59}$
GLASSZ8-2	${2.15}_{-0.04}^{+0.06}$	$-{19.75}_{-0.04}^{+0.03}$	$-{2.18}_{-0.02}^{+0.02}$	${8.56}_{-0.04}^{+0.05}$	${2.13}_{-0.18}^{+0.25}$	$-{1.46}_{-0.32}^{+0.17}$	$-{0.21}_{-0.05}^{+0.04}$	${0.24}_{-0.02}^{+0.02}$	${0.38}_{-0.02}^{+0.03}$	${4.59}_{-0.54}^{+0.58}$
Photometric Sample d
YD3	${2.79}_{-0.05}^{+0.08}$	$-{17.51}_{-0.19}^{+0.16}$	$-{2.04}_{-0.11}^{+0.09}$	${7.78}_{-0.28}^{+0.19}$	${0.27}_{-0.13}^{+0.18}$	$-{1.33}_{-0.47}^{+0.43}$	$-{1.38}_{-0.42}^{+0.59}$	${0.39}_{-0.26}^{+0.32}$	${0.09}_{-0.07}^{+0.18}$	⋯
YD6	${2.00}_{-0.04}^{+0.05}$	$-{18.23}_{-0.12}^{+0.11}$	$-{1.61}_{-0.08}^{+0.06}$	${8.61}_{-0.27}^{+0.19}$	${1.58}_{-0.74}^{+1.28}$	$-{1.32}_{-0.32}^{+0.57}$	$-{1.29}_{-0.44}^{+0.83}$	${0.69}_{-0.35}^{+0.24}$	${0.23}_{-0.13}^{+0.16}$	⋯
ZD1	${2.02}_{-0.04}^{+0.05}$	$-{17.08}_{-0.34}^{+0.25}$	$-{1.64}_{-0.13}^{+0.29}$	${8.19}_{-0.28}^{+0.28}$	${0.70}_{-0.37}^{+0.78}$	$-{1.35}_{-0.27}^{+0.40}$	$-{1.14}_{-0.61}^{+0.86}$	${0.55}_{-0.33}^{+0.51}$	${0.17}_{-0.05}^{+0.12}$	⋯
ZD4	${1.96}_{-0.04}^{+0.05}$	$-{17.59}_{-0.19}^{+0.20}$	$-{1.89}_{-0.16}^{+0.09}$	${8.09}_{-0.28}^{+0.31}$	${0.59}_{-0.30}^{+0.60}$	$-{1.30}_{-0.27}^{+0.49}$	$-{1.16}_{-0.51}^{+1.04}$	${0.35}_{-0.23}^{+0.36}$	${0.20}_{-0.07}^{+0.13}$	⋯
ZD5	${3.30}_{-0.13}^{+0.13}$	$-{17.55}_{-0.19}^{+0.15}$	$-{2.13}_{-0.04}^{+0.08}$	${7.84}_{-0.37}^{+0.27}$	${0.36}_{-0.20}^{+0.31}$	$-{1.29}_{-0.33}^{+0.35}$	$-{0.84}_{-0.62}^{+0.54}$	${0.16}_{-0.10}^{+0.21}$	${0.20}_{-0.11}^{+0.12}$	⋯
ZD7	${10.85}_{-0.38}^{+0.39}$	$-{17.44}_{-0.12}^{+0.10}$	$-{2.39}_{-0.04}^{+0.03}$	${7.27}_{-0.30}^{+0.26}$	${0.09}_{-0.05}^{+0.07}$	$-{1.49}_{-0.45}^{+0.52}$	$-{1.59}_{-0.31}^{+0.48}$	${0.10}_{-0.08}^{+0.20}$	${0.07}_{-0.05}^{+0.13}$	⋯
ZD9	${4.28}_{-0.11}^{+0.14}$	$-{17.53}_{-0.12}^{+0.11}$	$-{2.08}_{-0.02}^{+0.03}$	${7.53}_{-0.29}^{+0.42}$	${0.14}_{-0.07}^{+0.13}$	$-{1.42}_{-0.48}^{+0.66}$	$-{0.72}_{-0.90}^{+0.85}$	${0.20}_{-0.12}^{+0.24}$	${0.07}_{-0.04}^{+0.19}$	⋯
ZD10	${6.29}_{-0.34}^{+0.27}$	$-{15.33}_{-0.52}^{+0.39}$	$-{1.54}_{-0.29}^{+0.44}$	${7.50}_{-0.37}^{+0.28}$	${0.09}_{-0.06}^{+0.11}$	$-{1.07}_{-0.59}^{+0.46}$	$-{0.86}_{-0.92}^{+0.85}$	${0.55}_{-0.43}^{+0.61}$	${0.09}_{-0.07}^{+0.16}$	⋯

Notes.

^a Median magnification of the lens model by P. Bergamini 2023, (in preparation), calculated at the position of the source. Measurements are associated with 1σ uncertainties. ^b H ii bubble size, calculated with Equation (8). ^c Ionizing photon production efficiency (Section 3.3). ^d Sources selected by Zheng et al. (2014) that were not observed in our NIRSpec programs. Only those with photometric redshift consistent with z = 7.88 at 2σ are included. Redshift of the photometric sample is fixed to z = 7.88 in the SED analysis. t_*: mass-weighted age calculated over the posterior star formation history. Star formation rate is calculated by averaging the last 100 Myr of the posterior star formation history. All measurements are corrected for magnification.

Download table as:  ASCIITypeset image

To supplement our characterization of the overdensity (Section 4.1), we include the remaining eight photometric candidates presented by Zheng et al. (2014; Figure 1). As revealed by the spectroscopy, the confirmed sample consists both of ZDs and YDs. The ambiguity is likely due to the fact that the redshift of interest falls in the middle of the effective redshift ranges probed by the two color selections, which define 7 < z < 8 and 8 < z < 9 samples, respectively.

We fit the redshifts of the remaining photometric candidates with EAzY (Brammer et al. 2008). To exclude possible outliers, we only include those where the 2σ redshift uncertainty overlaps with z = 7.88. After this selection, we have nine photometric sources. The SEDs of the selected photometric sources are fitted as described above, with the redshift fixed to z = 7.88 (Table 2).

In Figure 3, we show the derived star formation histories of the confirmed member galaxies and the phot-z sample. All galaxies experience the peak of star formation in the last <100 Myr. The exception is YD7, which formed about a half of the total mass at ∼300 Myr prior to the observed redshift, making it a relatively old (with mass-weighted age t_* ∼ 200 Myr) and massive system.

Zoom In Zoom Out Reset image size

Overdense regions are generally known as a place of accelerated evolution (e.g., Dressler 1980; Thomas et al. 2010). It is thus of particular interest to investigate if there exists a luminous galaxy or quasar (M_UV < −22) in and/or near an overdense region at high redshift. The sources presented here are relatively faint in rest-frame UV, with UV absolute magnitude, M_UV, ranging from −15.3 to −20.1 mag. On the other hand, some of the sources have moderate dust attenuation (≳0.5 mag). Among those, YD4 is the most significant case with A_V = 1.1, which is consistent with the ∼4σ detection in dust continuum at its position (Laporte et al. 2017; see also Section 4.4). We compare A_V values for our sample to those of a reference field sample that consists of 13 photometric sources at z ∼ 8 from Leethochawalit et al. (2023). While the distribution of our sample is skewed toward higher A_V , the shift is not statistically significant. We find no significant differences in other properties.

We note that spectroscopic data are not included in the SED fitting process above; despite the predicted line fluxes of the best-fit SED model overall showing good consistency with the values measured in Section 3.1. The only exception is YD8, where the observed flux is ∼5× larger than what is predicted by the best-fit model for both Hβ and [O iii]. While detail investigation on this discrepancy is deferred to future work, we note that the emission line templates used here are optimized for standard stellar populations and do not include extreme components such as, e.g., active galactic nucleus.

Having direct measurements of optical emission lines brings us not only solid redshift confirmation but also direct insight into ionizing properties of the ISM. One of such measurements is the ionizing photon production per unit UV luminosity, or ionizing photon efficiency,

$$\begin{eqnarray}&&{\xi }_{\mathrm{ion}}={{\rm{N}}}_{\mathrm{LyC}}/{L}_{\mathrm{UV},1500}\,[\mathrm{Hz}\,{\mathrm{erg}}^{-1}]\end{eqnarray} \tag{ 1 }$$

where N_Lyc is the total ionizing photons of Lyman continuum and L_UV,1500 intrinsic UV luminosity density measured at rest-frame 1500 Å (e.g., Schaerer et al. 2016; Prieto-Lyon et al. 2022). The production efficiency is an important parameter, as the total ionizing output of galaxies can be simply parameterized (e.g., Madau et al. 1999; Robertson et al. 2013) as the product of ξ_ion and the fraction of ionizing photons that escape the interstellar medium into the intergalactic medium, f_esc. The direct measurement of N_LyC is not available for the redshift range of interest because the luminosity in the optical hydrogen recombination lines is proportional to the number of Lyman continuum photons absorbed in the galaxy. As a proxy, we adopt the following formula:

$$\begin{eqnarray}&&{N}_{\mathrm{LyC}}=2.1\times {10}^{12}{(1-{f}_{\mathrm{esc}})}^{-1}L{({\rm{H}}\beta )}_{\mathrm{unatten}.}\end{eqnarray} \tag{ 2 }$$

The intrinsic UV luminosity density is inferred from the best-fit SED. We correct the measured Hβ flux for attenuation by using the dust attenuation from the same SED modeling and assuming E(B − V)_neb/E(B − V)_stel. = 2.27 (Shivaei et al. 2020). Using the two equations above, we derive the production rate of ionizing photons that did not escape the galaxy, ξ_ion(1 − f_esc). The median value for the five objects available for the measurement is 〈ξ_ion(1 − f_esc)〉 = 4.6 × 10²⁵ erg Hz⁻¹. The measurements for individual galaxies are reported in Table 2.

We first investigate the spatial distribution of the member galaxies in A2744-z7p9OD. We use an updated version of the lens model presented by Bergamini et al. (2023), which includes the recent spectroscopic confirmation of a triply imaged z ∼ 10 Lyman break galaxies (Roberts-Borsani et al. 2022b) in the field (P. Bergamini et al. 2023, in preparation), to correct for the magnification by the foreground cluster. The 2D distribution of our sources in physical units is shown in Figure 4. After correcting for the lens magnification, we find that the confirmed sources are located within a circle of radius R ∼ 60 kpc in the source plane. The distribution on the sky is fairly elongated.

Zoom In Zoom Out Reset image size

Based on the derived spatial distribution in the source plane, we estimate overdensity by including only spectroscopically confirmed members. For the field reference, we use the luminosity function at z ∼ 8 derived in Bouwens et al. (2021) and integrated it down to M_UV = −19. Within a projected area of r = 60 kpc (∼12''), we find $\bar{n}={0.3}_{-0.1}^{+0.1}$, where the associated uncertainties reflect the 16th–84th percentile ranges of the luminosity function adopted. This gives an estimate of the overdensity $\delta ={24}_{-8}^{+12}$. While our estimate is lower than the one derived by Ishigaki et al. (2016, $\delta ={132}_{-51}^{+66}$), they included eight galaxies (including photometric candidates that are not confirmed in our study) in the central smaller region (r = 6''). Therefore, given that only the spectroscopic sample is included, our estimate is likely a conservative lower limit.

Second, we attempt to estimate the mass of the structure. Following previous work (e.g., Laporte et al. 2022), we can estimate the halo mass of the individual components from the halo-mass–galaxy-luminosity relation. Using the relation derived by Mason et al. (2022), we infer that the brightest member of the overdensity (ZD2, M_UV = −20.1) lives in a M_h ≈ (7 ± 2) × 10¹⁰ M_⊙ halo. Summing the halo mass of all the confirmed members, we obtain a lower limit to the total halo mass ≳4 × 10¹¹ M_⊙.

Lastly, we can take advantage of the spectroscopic data to obtain for the first time an estimate of the velocity dispersion of a protocluster at such high redshift. Given the small number of measured redshifts, we adopt a simple Gaussian estimator and bootstrap method to derive the uncertainty (Beers et al. 1990), obtaining 1100 ± 200 km s⁻¹. We caution the reader that the estimate should be treated with a degree of caution since the system is likely not virialized, and that in computing this quantity we are assuming the spread in redshift with respect to the mean is due to motion as opposed to distance along the line of sight. Nevertheless, we report it to assist future theoretical investigations.

With seven members being spectroscopically confirmed, it is of extreme interest to estimate the present-day halo mass of a system like A2744-z7p9OD. We attempt to estimate the total cluster mass at z = 0 by following the widely used formula

$$\begin{eqnarray}&&{M}_{z=0}=\bar{\rho }{V}_{\mathrm{cor}}(1+{\delta }_{m}),\end{eqnarray} \tag{ 3 }$$

(e.g., Steidel et al. 1998; Chiang et al. 2013), where ${V}_{\mathrm{cor}}=\tfrac{{V}_{\mathrm{obs}}}{C}$ is the redshift-space distortion-corrected comoving volume of the system, $\bar{\rho }$ is the mean matter density of the universe, and δ_m is the mass overdensity. The correction coefficients C and δ_m are linked as

$$\begin{eqnarray}&&1+b{\delta }_{m}=C(1+\delta ),\end{eqnarray} \tag{ 4 }$$

where b is the bias parameter, and δ is galaxy overdensity (Section 4.1). By adopting the linear interpolation presented in Ouchi et al. (2018), we adopt the bias parameter b = 6.5, at z = 7.9 of the total halo mass 1 × 10¹¹ M_⊙. The correction coefficient is expressed as

$$\begin{eqnarray}&&C=1+f-f{(1+{\delta }_{m})}^{1/3},\end{eqnarray} \tag{ 5 }$$

where f is a function of the mass density parameter, Ω_M (z), at z

$$\begin{eqnarray}&&f={{\rm{\Omega }}}_{M}^{4/7}(z).\end{eqnarray} \tag{ 6 }$$

By solving Equations (4) and (5), we obtain $C={0.56}_{-0.08}^{+0.08}$, and ${\delta }_{m}={2.00}_{-0.48}^{+0.55}$ for the galaxy overdensity measurement estimated in Section 4.1. For the same area used in Section 4.1 and the redshift interval δ z = 0.011, we estimate the comoving volume of A2744-z7p9OD to be V_obs = 9920 cMpc³. By substituting this in Equation (3), we obtain ${M}_{z=0}\,={2.2}_{-0.6}^{+0.9}\times {10}^{15}\,{M}_{\odot }$. This implies that A2744-z7p9OD would be expected on average to become a Coma-like system at z = 0, whereas 10¹⁴ M_⊙ is typically used as the threshold for a system to be called a "cluster" (e.g., Rosati et al. 2002; Overzier 2016).

Another way of predicting the present-day mass is to find and compare with overdense systems like A2744-z7p9OD in a simulation. Here we look into the EAGLE 100 Mpc "Reference" simulation (Schaye et al. 2015). We first extract all galaxies in the simulation with M_* > 10⁸ M_⊙ at z = 8 (N = 557), and trace them to their z = 0 descendants. To limit our analysis to the likely analogs of A2744-z7p9OD, we trace overdense regions that have six or more galaxies in a spherical region of r = 300 pkpc, at z = 8 (N = 19). At z = 0, the descendants of these 19 galaxies are hosted by halos in the mass range $\mathrm{log}{M}_{200c}/{M}_{\odot }$ = 13.5–14.5, which is consistent with a model overdensity identified by Ishigaki et al. (2016) in their cosmological simulation. While this indicates that A2744-z7p9OD could evolve into a system at the lower bound of typical clusters or even a group, it should be noted that the upper side of the resulting mass distribution above is likely limited by the simulation volume, as there is no object with $\mathrm{log}{M}_{200c}/{M}_{\odot }\gt 14.5$ in the entire simulation box. Estimating a precise mass of the structure and its future evolution would require simulations with sufficient resolution and astrophysical detail to resolve individual galaxy components matched in luminosity or stellar mass, while simultaneously probing sufficient volumes, to include multiple structures of this kind to average out the expected stochasticity (Chiang et al. 2013), or mapping to a typical dark-matter halo rarity then followed across cosmic time (Trenti et al. 2008). A dedicated study is beyond the scope of this paper and left for future work.

Our very first attempt of spectroscopic follow-up on the bright members in the core region already confirmed seven member galaxies at 100% success rate. For further characterization of the confirmed overdensity, a sample of galaxies at larger extent would be required. The progenitors of massive clusters are typically spread over several Mpc, and thus to robustly estimate the mass of the descendant, one would require a survey covering a much larger area of the sky (e.g., Overzier et al. 2009; Contini et al. 2016). In fact, we found in the simulation that the z = 0 mass distribution skews at a higher mass when the search region is defined by a larger radius (∼2 pMpc) that contains more numerous galaxies (N = 35 or more).

The absence of strong Lyα emission provides insight into the intergalactic medium (IGM) properties surrounding the protocluster. None of our spectroscopically confirmed sources show a clear Lyα line (Figure 2). To quantify the nondetections, we estimate the upper limit on rest-frame equivalent widths of the line, EW₀(Lyα), following Hoag et al. (2019), Morishita et al. (2020):

$$\begin{eqnarray}&&{\mathrm{EW}}_{0}(\mathrm{Ly}\alpha )=\displaystyle \frac{{F}_{\mathrm{Ly}\alpha }}{{f}_{\mathrm{cont}.}}\displaystyle \frac{1}{(1+z)}.\end{eqnarray} \tag{ 7 }$$

For f_cont. we use the continuum model derived from our SED fitting analysis. We replace the nondetected Lyα flux with the limiting flux estimated over the instrumental resolution (∼100 Å, or 2700 km s⁻¹) at the wavelength of Lyα. The resulting range of upper limits is ∼16–28 Å (2σ; Table 1). The two galaxies (ZD2, ZD3) observed in the Early Science Release programs are not available for the EW measurement as the wavelength range of interest falls in the detector gap.

The lack of strong Lyα emission is perhaps not surprising given the redshift z = 7.9 of the host galaxies, where inferences on the ionization state of the IGM find neutral fractions in excess of x_HI > 70% (Hoag et al. 2019; Mason et al. 2019). From the measured Lyα EW limits of the four galaxies, and their absolute magnitudes, we estimate the volume-averaged neutral hydrogen fraction of IGM to be >0.45 (at 68% CL) using the same methodology presented by Mason et al. (2018). This is consistent with previous work on the cosmic average, within the uncertainties (Mason et al. 2019). We note that this analysis assumes the observations are independent sightlines. A more realistic analysis including their correlation within the same physical region would likely recover a slightly lower limit, but is beyond the scope of this work. A larger number of spectroscopic measurements in the protocluster are needed to refine this limit and identify potential differences with regard to the cosmic average.

In such a highly neutral environment, large ionized bubbles are expected to be extremely rare (e.g., Mesinger & Furlanetto 2007). Even around regions of comparable overdensity containing sources of similar magnitude ranges, large reionization simulations predict median bubble sizes to be smaller than 1 pMpc where x_HI > 70% (T.-Y. Lu et al. 2023, in preparation, using the evolution of structure reionization simulations; Mesinger et al. 2016). If the bubble size is below ∼1 pMpc, the redshift along the line of sight is not sufficient for Lyα to escape, and in fact Lyα transmission is ≲20% at its line center (Mason & Gronke 2020; Qin et al. 2021). Assuming a Gaussian-shape emission line with velocity offset of Lyα from systemic of 200 km s⁻¹ (which is likely an overestimate in these low luminosity galaxies, Mason et al. 2018) and FWHM equal to the velocity offset, the total fraction of transmitted Lyα flux is expected to be <40% for a 1 pMpc bubble and <30% for a 0.75 Mpc bubble. Assuming the average Lyα EW = 30 Å for M_UV ∼ −19.75 galaxies at z ∼ 6 from De Barros et al. (2017) as the "emitted" EW, we would thus expect to observe Lyα, after transmission through the IGM, with EW < 12 Å, below our detection threshold (Table 1).

We can verify our theoretical expectation by estimating the radius of an H ii region, R_HII , ionized by UV photons of a single galaxy using the equation in Haiman & Loeb (1997; also Endsley & Stark 2022):

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{dR}}_{{\rm{H}}{\rm\small{II}}}}{{dt}} & = & \displaystyle \frac{\langle {\xi }_{\mathrm{ion}}\rangle \langle {f}_{\mathrm{esc}}\rangle {L}_{\mathrm{UV}}}{4\pi {R}_{{}_{{\rm{H}}{\rm\small{II}}}}^{2}\ {\bar{n}}_{{\rm{H}}{\rm\small{I}}}(z)}+{R}_{{}_{{\rm{H}}{\rm\small{II}}}}\ H(z)\\ & & -\,{R}_{{}_{{\rm{H}}{\rm\small{II}}}}\ {\alpha }_{B}\ {\bar{n}}_{{\rm{H}}{\rm\small{I}}}(z)\ \displaystyle \frac{{C}_{{\rm{H}}{\rm\small{I}}}}{3}.\end{array}\end{eqnarray} \tag{ 8 }$$

For simplicity, this equation assumes that the ionizing bubble is spherically symmetric and created by a single source at its center. We adopt a case (B) recombination coefficient, α_B = 2.59 × 10⁻¹³ cm³s^–1 (Osterbrock 1989), and ionizing photon escape fraction of 〈f_esc〉 = 0.2, and an IGM H i clumping factor C_HI = 3.0 (Shull et al. 2012; Robertson et al. 2013). For the ionizing photon production efficiency, we use those derived in Section 3.3 for those with Hβ flux measurements available. We adopt the median value for the remaining sample.

By using the derived star formation history and luminosity presented in Section 3.2, we estimate the bubble size for the spectroscopically confirmed individual sources and for the photometrically selected sample identified in Section 3.2. The estimated bubble sizes are ≪1 Mpc for most of the sample, due to the relatively low UV luminosity of the galaxies (Table 2). Given that their separation in the source plane is of order 60 pkpc, we also estimate the bubble size by considering the cumulative effect of all the confirmed sources as if they were colocated. Even in this case, we find it to be R ∼ 0.78 Mpc, i.e., insufficient to allow Lyα to escape. Even the inclusion of all photometric candidates is not significantly changing the estimate, as those are mostly fainter than the confirmed members (also see their individual estimates in Table 2).

In conclusion, our bubble size estimate is consistent with the nondetection of Lyα in our confirmed sources. By comparison, at slightly lower redshifts, Endsley & Stark (2022) estimated bubble sizes of 0.7–1.1 Mpc for UV-bright (−M_UV ∼ 20–22 mag) galaxies at z = 6.6–6.9. Lyα was detected in nine out of ten galaxies, showing that both an overdense environment and sufficient ionizing photon flux are required to produce an ionized bubble large enough to allow significant transmission of Lyα photons (see the comparison of Lyα detections in UV-bright and fainter galaxies in Roberts-Borsani et al. 2022c).

Our estimated lower limit of the mass of A2744-z7p9OD is comparable to those of previously known protoclusters and protoscluster candidates using similar methods, including the recently reported candidate behind the SMACS0723 cluster (Laporte et al. 2022), where two of the member candidates are spectroscopically confirmed to be z = 7.66. The real breakthrough of our observations, however, is the sheer number of spectroscopically confirmed redshift measurements, which allow us to establish secure membership to the protocluste and get a first estimate of the its velocity dispersion. This clearly provides a glimpse of the power of JWST to add unprecedented detail to studies of the progenitors of today's large-scale structures.

The spectroscopic redshift of A2744-z7p9OD is in agreement with previous photometric redshift estimates, but in apparent tension with the redshift reported for YD4 (z = 8.38) based on Lyα and ALMA [C ii] 158 μm and [O iii] 88 μm emission (Laporte et al. 2017, 2019; Carniani et al. 2020). A likely explanation of the apparent tension is a line-of-sight superimposition of sources at similar redshifts. This is a common occurrence in the photometric identification of overdensities and should be kept in mind when considering protocluster candidates without spectroscopy. As shown in Figure 5, YD4 is close on the sky (separation ∼05) to a secondary Y-dropout source (YD6; Zheng et al. 2014), which falls outside our NIRSpec spectroscopic apertures (in the GLASS and DDT observations) but lies within the VLT/X-Shooter long-slit and indistinguishable from YD4 at ALMA resolution (in the case of [C ii]). We hypothesize that the source detected at z = 8.38 is actually in the background of the protocluster and likely associated with YD6 (estimated to be at z_phot = 8.3 ± 0.2; Zheng et al. 2014). In fact, the [O iii] 88 μm flux appears better centered on the fainter counterpart while the detection of dust appears associated with YD4 (see Figure 2 in Laporte et al. 2019) is consistent with both the large dust quantity (A_V = 1.1) and red UV slope (β_UV = −1.3) estimated for YD4 in Section 3.2 and the discrepant spectroscopic redshifts. More extensive spectroscopic coverage is needed to confirm the hypothesis.

Zoom In Zoom Out Reset image size