The JWST Advanced Deep Extragalactic Survey: Discovery of an Extreme Galaxy Overdensity at z = 5.4 with JWST/NIRCam in GOODS-S

Our imaging data consist of: (1) deep optical imaging taken with HST/ACS in five photometric bands (F435W, F606W, F775W, F814W, and F850LP); and (2) deep infrared imaging taken with JWST/NIRCam in 14 photometric bands (F090W, F115W, F150W, F182M, F200W, F210M, F277W, F335M, F356W, F410M, F430M, F444W, F460M, and F480M).

The HST/ACS mosaics used here were produced as part of the Hubble Legacy Fields project v2.0 and include observations covering a $25^{\prime} \times 25^{\prime}$ area over the GOODS-S field (Illingworth et al. 2016; Whitaker et al. 2019). The JWST/NIRCam data were obtained by the JWST Advanced Deep Extragalactic Survey (JADES; Eisenstein et al. 2023; PID: 1180 and 1210) and the JWST Extragalactic Medium-band Survey (JEMS; Williams et al. 2023; PID: 1963) in 2022 September and October. The JADES observations consist of a deep mosaic covering a $4\buildrel{\,\prime}\over{.} 4\times 6\buildrel{\,\prime}\over{.} 2$ area with nine filters (F090W, F115W, F150W, F200W, F277W, F335M, F356W, F410M, and F444W) and a medium region covering an additional $6\buildrel{\,\prime}\over{.} 1\times 6\buildrel{\,\prime}\over{.} 5$ area with eight filters (F090W, F115W, F150W, F200W, F277W, F356W, F410M, and F444W). The JEMS observations consist of two $2\buildrel{\,\prime}\over{.} 2\times 2\buildrel{\,\prime}\over{.} 2$ regions with five filters (F182M, F210M, F430M, F460M, and F480M), all lying in JADES coverage. For all of the subsequent analysis, we do not require any of our objects to have JEMS observations, but we use these data when available.

A detailed description of the JWST/NIRCam imaging data reduction and mosaicking will be presented in a forthcoming paper from the JADES Collaboration (S. Tacchella et al. 2023, in preparation). We briefly summarize here the main steps of the reduction and mosaicking process. The data are initially processed with the standard JWST calibration pipeline. ¹⁷ Customized steps are included to aid in the removal of "1/f" noise, "wisp" artifacts, "snowball" artifacts, and persistence from previous observations (see also Rigby et al. 2023). The JWST Calibration Reference Data System context map jwst_1008.pmap is used, including the flux calibration for JWST/NIRCam from Cycle 1. The background from the sky is modeled and removed using the BackGround2D class from photutils (Bradley et al. 2022). Finally, the image mosaics for each of the 14 JWST/NIRCam filters are registered to the Gaia Data Release 3 frame (Gaia Collaboration et al. 2023) and resampled onto the same world coordinate system (WCS) with a 30 mas pixel⁻¹ grid. Assuming a circular aperture with a diameter of 03, the 5σ point-source detection limit in the F200W filter is m ≈ 30.0 AB mag and m ≈ 29.0 AB mag for the deep and medium regions, respectively.

A detailed description of the JWST/NIRCam source detection was outlined in Robertson et al. (2023) and will be presented in detail in another forthcoming paper from the JADES Collaboration (B. Robertson et al. 2023, in preparation). We briefly summarize here the main steps of the source detection process. Six image mosaics (F200W, F277W, F335M, F356W, F410M, and F444W) are initially stacked using the corresponding error images and inverse-variance weighting to produce a single detection image. These filters were chosen in order to avoid biasing our catalog against short-wavelength (SW) dropouts (e.g., dropouts in F090W, F115W, or F150W). In this detection image, we construct a source catalog by selecting contiguous regions of greater than five pixels with signal-to-noise ratios S/N > 3 and applying a standard Source Extractor (SExtractor; Bertin & Arnouts 1996) deblending algorithm with parameters nlevels = 32 and contrast = 0.001 using photutils (Bradley et al. 2022). Finally, we perform forced convolved photometry at the source centroids in all HST/ACS and JWST/NIRCam photometric bands, assuming elliptical Kron apertures with parameter = 1.2 (i.e., Kron small) and parameter = 2.5 (i.e., Kron large). To correct for potential missing light, we rescale the Kron small photometry by the flux ratio of Kron large to Kron small in the F444W filter. Using model point-spread functions from the TinyTim (Krist et al. 2011) package for HST/ACS and the WebbPSF (Perrin et al. 2014) package for JWST/NIRCam, we apply aperture corrections assuming point-source morphologies. Uncertainties are estimated by placing random apertures across regions of the image mosaics to compute a flux variance (e.g., Labbé et al. 2005; Quadri et al. 2007; Whitaker et al. 2011), which are summed in quadrature with the associated Poisson uncertainty for each detected source.

Using the previously described photometry, we measure photometric redshifts with the template-fitting code EAZY (Brammer et al. 2008). A more detailed description of this procedure will be discussed in a forthcoming paper from the JADES Collaboration (K. N. Hainline et al. 2023, in preparation). We briefly summarize here the main steps of the photometric redshift process. EAZY uses a chi-square (χ²) minimization technique to model the broadband spectral energy distributions (SEDs) for galaxies using linear combinations of galaxy templates. It is designed to be both fast and flexible, and has been used extensively in the literature to model the photometric redshifts of galaxies (e.g., Newman et al. 2013; Skelton et al. 2014; Bouwens et al. 2015).

We fit all of the available photometry for each object with EAZY, assuming the rescaled Kron small photometry described in Section 2.1. For objects in the ultra-deep field (UDF), this includes photometry in 19 filters. For objects in the JADES deep region but not in the UDF, this includes photometry in 14 filters. For objects in the JADES medium region, this includes photometry in 13 filters. We utilize 16 templates in total to perform the fitting, including the nine EAZY "v1.3" templates, two additional templates for simple stellar populations with ages of 5 and 25 Myr, and five more templates with strong nebular continuum emission that were created using the Flexible Stellar Population Synthesis code (FSPS; Conroy et al. 2009; Conroy & Gunn 2010). These templates span a large range of stellar population properties and include contributions from both nebular continuum and line emission, as well as obscuration from dust.

While the photometric calibration of JWST/NIRCam has improved significantly over the course of the last few months, there is still some uncertainty with these calibrations that needs to be taken into account. To this end, we iteratively calculate the photometric offset from the EAZY templates compared to the true JWST/NIRCam photometry, using a sample of galaxies with S/N between 5 and 20 in F200W. These photometric offsets are relatively small (on the order of a few percent for both HST/ACS and JWST/NIRCam) and are subsequently applied to the entire photometric catalog. We choose not to adopt any apparent magnitude priors, but we do make use of the template error file "template_error.v2.0.zfourge."

The primary measurements used here are the EAZY "z_a" and "z_peak" redshifts. The former corresponds to the fit where the likelihood is maximized (χ² is minimized), while the latter corresponds to the weighted average of the fits where the likelihood is maximized (the probability is equal to the integral of the likelihood). We allow EAZY to fit across the redshift range of z = 0.2–22 with a redshift step size of ${\rm{\Delta }}z=0.01\left(1+z\right)$. To test the accuracy of these photometric redshifts, we compare these predictions with existing spectroscopic redshifts in the GOODS-S field from the Multi Unit Spectroscopic Explorer (MUSE; Inami et al. 2017; Urrutia et al. 2019). While the MUSE spectroscopic redshifts are biased toward the brightest objects detected by JWST at z < 7, we found catastrophic outlier fractions of only 5% (Rieke et al. 2023b) when comparing to the highest-quality spectroscopic redshifts available with MUSE. The catastrophic outlier fraction is defined to be the fraction of objects that satisfy Equation (1):

$$\begin{eqnarray}&&\displaystyle \frac{| {z}_{\,\mathrm{spec}}-{z}_{\,\mathrm{phot}}| }{1+{z}_{\,\mathrm{spec}}}\gt 0.15.\end{eqnarray} \tag{ 1 }$$

To perform an accurate and efficient targeted emission-line search within the available spectroscopic data, we require a sample of relatively bright objects, since these are the only objects for which we expect to detect an emission line (see Section 2.3). We also require these objects to have tight photometric redshift constraints, which allow for spectroscopic redshift confirmation using only a single line detection. Most objects with tight photometric redshift constraints have emission lines that fall in one of the medium-band filters (e.g., F410M), which allows for tight constraints when paired with the broadband filter coverage of JADES (e.g., F444W). Our selection criteria for the final photometric catalog consist of the following: m < 28.5 AB mag in F444W assuming elliptical Kron apertures with parameter = 2.5, 4.5 < z_a < 9.5, 4.5 < z_peak < 9.5, Δz ₁ < 1, and Δz ₂ < 2. The first EAZY confidence interval (Δz ₁) is defined to be the difference between the 16th and 84th percentiles of the photometric redshift posterior distribution and is roughly twice the standard deviation. The second EAZY confidence interval (Δz ₂) is defined to be the difference between the 5th and 95th percentiles of the photometric redshift posterior distribution and is roughly four times the standard deviation.

Our spectroscopic data consist of WFSS observations taken with JWST/NIRCam in the F444W filter (λ = 3.9–5.0 μm). These data were obtained by the First Reionization Epoch Spectroscopic COmplete Survey (FRESCO; PI: Oesch; PID: 1895) in 2022 November. The FRESCO observations cover an $8\buildrel{\,\prime}\over{.} 2\times 8\buildrel{\,\prime}\over{.} 6$ area using the row-direction grisms on both modules of JWST/NIRCam (Grism R; R ≈ 1600). The total overlapping area between the JADES and FRESCO footprints is ≈41 square arcminutes. The total spectroscopic observing time for FRESCO in GOODS-S is ≈16 hr with a typical on-source time of ≈2 hr. The 3σ unresolved emission-line detection limit around 4.2 μm in the F444W filter is ∼1.2 × 10⁻¹⁸ erg s⁻¹ cm⁻², which corresponds to a star formation rate (SFR) detection limit of ∼2.1 M_⊙ yr⁻¹ at z = 5.4 using the conversion factor from Kennicutt & Evans (2012).

A detailed description of the JWST/NIRCam grism data reduction can be found in Sun et al. (2023). We briefly summarize here the main steps of the reduction process. The data are initially processed with the standard JWST calibration pipeline. ¹⁸ We assign WCS to the rate files, perform flat-fielding, and subtract out the sigma-clipped median sky background from each individual exposure after the "ramp-to-slope" fitting in the calibration pipeline. Because we are interested in conducting a targeted emission-line search, and we do not expect any of our sources to have a strong continuum due to their general faintness (m = 27–28 AB mag), we utilize a median-filtering technique to subtract out any remaining continuum or background on a row-by-row basis, following the methodology of Kashino et al. (2023). This produces emission-line maps for each grism exposure that are void of any continuum. Although this median-filtering technique is able to properly remove continuum contamination, it sometimes oversubtracts signal in the spectral regions immediately surrounding the brightest emission lines (e.g., emission from [N ii] on either side of Hα). The full widths at half maximum (FWHMs) of these emission lines are relatively small and typically on the order of a few pixels, which means that the Hα line flux is preserved by the 9 pixel central gap of the median filter and is not affected by the aforementioned oversubtraction. We further remove the "1/f" noise using the tshirt/roeba algorithm ¹⁹ in both the row and column directions.

We extract two-dimensional (2D) grism spectra using the reduced emission-line maps for all of the objects that are part of the final photometric catalog discussed in Section 2.2. SW parallel observations were conducted in two photometric bands (F182M and F210M) and are used for both astrometric and wavelength calibration of the long-wavelength spectroscopic data. We use the spectral tracing and grism dispersion models (Sun et al. 2023) that were produced using the JWST/NIRCam commissioning data of the Large Magellanic Cloud (PID: 1076), which are also outlined in Wang et al. (2023). We additionally use the flux calibration models that were produced using JWST/NIRCam Cycle 1 absolute flux calibration observations (PID: 1536/1537/1538).

Using the already extracted 2D spectra, we further extract one-dimensional (1D) grism spectra using a boxcar aperture, assuming a height of 5 pixels (031). We subsequently identify >3σ peaks automatically in the 1D spectra, assuming various bin sizes (integer units of nanometers from 1 to 8) and fit these detected peaks with Gaussian profiles. For each line that is detected with S/N > 3, we tentatively assign a line identification of either Hα or [O III]λ5008, whichever one minimizes the difference between the best-fit photometric redshift and the tentative spectroscopic redshift. For example, if a line were detected at λ = 4.2 μm and the best-fit photometric redshift is z_phot = 5.8, then the initial line identification would be Hα, since the predicted wavelength of this line would be at λ = 4.5 μm, which is closer to the observed wavelength than the predicted wavelength of [O III]λ5008 (λ = 3.4 μm). Visual inspection is performed on each of these tentative spectroscopic redshift solutions to remove spurious detections caused by either noise or contamination. For sources that pass our visual inspection and have secure line detections, we optimally re-extract the 1D spectra using the F444W surface brightness profile (Horne 1986) and once again fit these detected peaks with Gaussian profiles. According to the grism wavelength calibration uncertainty, the typical absolute uncertainties of our spectroscopic redshifts are Δz_spec = 0.001.

Our final spectroscopic sample includes N = 81 objects at z = 5.2–5.5 with >3σ detections of Hα from the FRESCO spectra. This redshift range was chosen to ensure that Hα would fall in the F410M filter, which is the only medium-band filter for which we have uniform coverage, providing a sanity check for the derived emission-line fluxes through a comparison with the F410M excess relative to F444W. Our final spectroscopic sample represents a subset of a larger spectroscopic sample of galaxies from both GOODS fields across a much broader redshift range (F. Sun et al. 2023, in preparation). For the majority of galaxies in our final spectroscopic sample, neither of the [N ii] lines were detected, partially as a result of the aforementioned median-filtering technique (Kashino et al. 2023), but primarily because the line ratio [N ii]/Hα is typically low at these redshifts (e.g., Cameron et al. 2023). The NIRCam cutout images alongside the 2D and 1D extracted spectra for these objects are shown in Appendix A. Figure 1 shows the distribution of spectroscopic redshifts for these N = 81 objects, while Figure 2 shows the on-sky distribution in angular units. These distributions enabled us to visually identify an overdensity of galaxies around z = 5.4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Following the technique described in Calvi et al. (2021), we use a Friends-of-Friends (FoF) algorithm to identify the overdensity after looking for three-dimensional (3D; two spatial, one spectral) structural groupings (see also Huchra & Geller 1982; Eke et al. 2004; Berlind et al. 2006). This algorithm iteratively selects groups, which consist of one or more galaxies that have projected separations and line-of-sight (LOS) velocity dispersions below the adopted linking parameters (projected separation d_link = 500 kpc, chosen to be roughly the virial radius of a typical galaxy cluster, and LOS velocity dispersion σ_link = 500 km s⁻¹, chosen to be roughly the velocity dispersion of such a cluster). These groupings do not depend strongly on the adopted linking parameters, producing similar results when varying either the projected separation or the LOS velocity dispersion by a factor of a few.

We identify one large-scale structure consisting of N = 39 galaxies out of the N = 81 galaxies that are part of our final spectroscopic sample. Throughout the rest of this work, we refer to the remaining N = 42 galaxies at z = 5.2–5.5 as field galaxies, which consist of: (1) isolated galaxies; and (2) those in smaller groups, as determined by the FoF algorithm. ²⁰ The average spectroscopic redshift of the large-scale structure is z = 5.386, spanning a relatively narrow redshift range of 5.374 < z < 5.398. The maximum on-sky separation of the clustered galaxies is roughly 36, corresponding to a physical separation of 8.6 cMpc.

Following the methodology of Chiang et al. (2013), we calculate 3D galaxy overdensities (δ_gal = n_gal/〈n_gal〉 − 1), where n_gal is the number density of galaxies and 〈n_gal〉 is the ensemble average number density of galaxies, for each of the N = 81 galaxies that are part of our final spectroscopic sample. To calculate these values, we assume a tophat-weighted spherical window that has a comoving volume equal to ${\left(15\,\mathrm{cMpc}\right)}^{3}$. The average and standard deviation of the constituent overdensity values for the large-scale structure identified here are 〈δ_gal〉 = 9.2 ± 1.0. This structure is an extreme overdensity at z = 5.4, ∼10 times more dense in 3D than the ensemble average at z = 5.2–5.5. For comparison, Chiang et al. (2013) found that an average overdensity value of 〈δ_gal〉 = 3.04 identifies structures within cosmological simulations as protocluster candidates with 80% confidence. This value is for the z = 5 SFR-limited sample (SFR > 1 M_⊙ yr⁻¹), which is the sample that is most similar to our own in terms of selection. Throughout this work, we define a protocluster as a structure that will eventually collapse into a galaxy cluster at z = 0.

Throughout the rest of Section 3, the distinction between the N = 39 "confirmed members of overdensity" and the N = 42 "confirmed members of field" will be used. The confirmed members of the overdensity are shown in Figure 1 (Figure 2) by the turquoise histograms (turquoise pluses), while the confirmed members of the field are shown by the gray histograms (gray points). The median redshift or position of the overdensity is given by the solid magenta line or the magenta point. In both of these figures, the overdensity members appear much more clustered when compared to the field members. The confirmed members of the overdensity are additionally shown in the left panel of Figure 3, color-coded by their spectroscopic redshift. We identify a spatial and kinematic bimodality within the overdensity at z = 5.4, which we return to in Section 3.3. To assign objects between the two components of the bimodality, we adopt an iterative process that minimizes the 3D separation within each of these components, finding N = 14 galaxies that are part of the first component (${\rm{\Delta }}v=-{340}_{-50}^{+140}\,\mathrm{km}\,{{\rm{s}}}^{-1}$) and N = 25 galaxies that are part of the second component (${\rm{\Delta }}v=+{190}_{-150}^{+120}\,\mathrm{km}\,{{\rm{s}}}^{-1}$). The median positions of these two groups are given by the black plus and cross in the left panel of Figure 3. Velocity offsets are calculated relative to the median spectroscopic redshift of the overdensity.

Zoom In Zoom Out Reset image size

Following the methodology outlined in Tacchella et al. (2022), we utilize the SED fitting code Prospector (v1.1.0; Johnson et al. 2021) to infer the stellar populations for the N = 81 objects that are part of our final spectroscopic sample. Fits are performed on the rescaled Kron small photometry described in Section 2.1, while the redshift is fixed at the spectroscopic redshift determined in Section 2.2. Prospector uses a Bayesian inference framework and we choose to sample posterior distributions with the dynamic nested sampling code dynesty (v1.2.3; Speagle 2020).

We use the Flexible Stellar Population Synthesis (FSPS; Conroy et al. 2009; Conroy & Gunn 2010) code via python-FSPS (Foreman-Mackey et al. 2014) with the Modules for Experiments in Stellar Astrophysics Isochrones and Stellar Tracks (MIST; Choi et al. 2016; Dotter 2016), which make use of the Modules for Experiments in Stellar Astrophysics (MESA) stellar evolution package (Paxton et al. 2011, 2013, 2015, 2018). We assume the MILES stellar spectral library (Falcón-Barroso et al. 2011; Vazdekis et al. 2015) and adopt a Chabrier (2003) initial mass function. Absorption by the intergalactic medium (IGM) is modeled after Madau (1995), where the overall scaling of the IGM attenuation curve is set to be a free parameter. Dust attenuation is modeled using a two-component dust attenuation model (Charlot & Fall 2000) with a flexible attenuation curve, where the slope is tied to the strength of the ultraviolet (UV) bump (Kriek & Conroy 2013). Nebular emission (both from emission lines and continuum) is self-consistently modeled with the spectral synthesis code Cloudy (Byler et al. 2017).

To test the robustness of the inferred stellar populations and star formation histories (SFHs), we assume three different models for the SFH to see how our results depend on the assumed prior. Two of these models are nonparametric (one with the "continuity" prior, the other with the "bursty continuity" prior), while the third is parametric (with the shape of a delayed-tau function). For each of the nonparametric models, we assume that the SFH can be described by N_SFR distinct time bins of constant star formation. The time bins are specified in units of lookback time and the number of distinct bins is fixed at N_SFR = 6. The first two bins are fixed at 0–30 Myr and 30–100 Myr, while the last bin is fixed between 0.15 t_univ and t_univ, where t_univ is the age of the Universe at the galaxy's spectroscopic redshift, measured with respect to the formation redshift z_form = 20. The rest of the bins are spaced equally in logarithmic time between 100 Myr and 0.85 t_univ. Both the total stellar mass and the ratios of adjacent time bins are set to be free parameters. A summary of the parameters and priors associated with this Prospector model is presented in Table B1. We adopt the results from the nonparametric SFH with the "continuity" prior as fiducial, since the Bayesian evidence does not strongly favor one model over another for the vast majority of the galaxies considered here.

Figure 4 shows the star-forming main sequence for the N = 81 objects at z = 5.2–5.5 that are part of the final spectroscopic sample identified in Section 2.3. The confirmed members of the field are given by the gray points, while the confirmed members of the overdensity are given by the turquoise points. The reported stellar masses and SFRs are derived from the Prospector fits using the nonparametric SFH with the "continuity" prior, where the SFRs are averaged over the last 100 Myr of lookback time. These stellar masses and SFRs are reported in Table C1, which gives a summary of the physical properties for our final spectroscopic sample. Based on these stellar masses and the mass–size relation reported in Shibuya et al. (2015), we find that the minimum separation between the sources in our sample is always larger than twice the effective radius of the sources in their sample. Therefore, each object that is part of our final spectroscopic sample is likely an individual star-forming galaxy rather than an individual star-forming clump within a much larger galaxy, although some are clearly merging systems in the final phase of coalescence. To be used as a point of comparison, the empirical star-forming main sequence at z = 5.4 derived by Popesso et al. (2023) is given by the solid black line. Additionally, the maximum allowed SFR assuming all of the stellar mass was formed in the last 100 Myr of lookback time is given by the solid magenta line.

Zoom In Zoom Out Reset image size

We find that nearly all of the N = 81 objects at z = 5.2–5.5 that are part of our final spectroscopic sample agree with the empirically derived star-forming main sequence at z = 5.4 given by Popesso et al. (2023) within 1σ, despite our sample being biased as a result of our requirement to detect the Hα emission line at greater than 3σ. If we were to instead use Hα-based SFRs derived with the conversion factor from Kennicutt & Evans (2012), the objects within our final spectroscopic sample would be shifted upward by Δ ∼ 0.5 dex. However, this is likely because the canonical hydrogen ionizing photon production efficiency (ξ_ion ∝ L_Hα /L_UV) used in Kennicutt & Evans (2012) is only ξ_ion ∼ 10^25.1 erg⁻¹Hz, lower than measurements at z = 5–6 by Δξ_ion ∼ 0.5 dex (e.g., Bouwens et al. 2016; Ning et al. 2023; Sun et al. 2023). Furthermore, if we were to instead use UV-based SFRs derived with the conversion factor from Kennicutt & Evans (2012) alongside measurements of the rest-UV magnitudes derived as ν L_ν at λ_rest = 1500 Å, it would not change the distribution of our final spectroscopic sample in Figure 4. This is expected, since the Prospector-based SFRs largely rely on the available rest-UV photometry.

The most notable examples of galaxies disagreeing with the empirical relation include one well above the main sequence at high stellar mass (JADES−GS+53.13859−27.79025) and two somewhat below the main sequence at intermediate stellar mass (JADES−GS+53.06799−27.80816 and JADES−GS+53.18328−27.77894), despite large uncertainties in the inferred SFRs. We note that the galaxy well above the main sequence at high stellar mass is an active galactic nucleus (AGN), originally identified as a broadline Hα emitter in Matthee et al. (2023) with a line width of 2200 ± 500 km s⁻¹. Combined with our relatively low SFR detection limit (∼2.1 M_⊙ yr⁻¹; see Section 2.3), the fact that our sample agrees with the empirical relation suggests that we are sampling the bulk of the star-forming population at these redshifts, despite our selection criteria. However, we should mention that we are likely missing some amount of dusty star-forming galaxies (DSFGs), which historically have been used as tracers to identify protocluster candidates at these redshifts (for a review of environmental galaxy evolution, see Alberts & Noble 2022).

To test the impact of the assumed SFH, we compare the Prospector-derived stellar masses and SFRs for the three different SFH models. We remind the reader that for Figure 4 and Table C1, the assumed SFH is nonparametric with the "continuity" prior. Compared to the nonparametric SFH with the "bursty continuity" prior, the derived stellar masses are a bit larger (Δ ≈ 0.1 dex) with large scatter (σ ≈ 0.7 dex), while the derived SFRs are broadly consistent (Δ ≈ 0.0 dex) with large scatter (σ ≈ 0.8 dex). Compared to the parametric SFH, both the derived stellar masses and SFRs are broadly consistent (Δ ≈ 0.0 dex) with large scatter (σ ≈ 0.8 dex and σ ≈ 1.0 dex, respectively). We find that the inferred Prospector parameters considered here do not depend strongly on the assumed SFH.

Figure 5 shows the comparison of SFRs in the two most recent time bins (corresponding to lookback times of 0–30 Myr and 30–100 Myr, respectively) for the N = 81 objects at z = 5.2–5.5. Once again, the confirmed members of the field are given by the gray points, while the confirmed members of the overdensity are given by the turquoise points. A constant SFH is given by the solid black line, which assumes that the average SFRs in the two most recent time bins are equal. Values above this line represent a falling SFH over the last 100 Myr, while values below this line represent a rising SFH. The twin axes to the top and right are similar to those shown in Figure 4, where the median values for the confirmed members of the field (overdensity) are given by the gray (turquoise) dashed lines.

Zoom In Zoom Out Reset image size

By comparing members of the overdensity with members of the field in Figures 4 and 5, we can begin to explore the impact of environment on galaxy formation and evolution at z = 5.2–5.5. Table 1 gives a summary of percentiles for some of the Prospector-inferred physical parameters for both the members of the field and members of the overdensity. In the star-forming main sequence and the summary of percentiles, overdensity members appear to have larger inferred stellar masses and SFRs when compared to the field members. We further compare the distributions of these inferred parameters by performing Kolmogorov–Smirnov (KS) and Anderson–Darling (AD) tests, which are two-sided tests for the null hypothesis that two independent samples are drawn from the same continuous distribution. The results of these tests indicate that the stellar masses and SFRs of the field and overdensity samples exhibit statistically significant differences at roughly the 2σ–3σ level (corresponding to p-values of 0.003 < p < 0.05). However, to fairly compare these two samples, the stellar mass distributions must be similar, such that any observed differences can be attributed to external processes (e.g., environment), rather than internal processes (e.g., feedback). To create a mass-matched sample of field galaxies, we select the field galaxy that is closest in stellar mass to each overdensity galaxy, which produces a sample of galaxies with stellar masses that are consistent with being drawn from the same parent distribution as the overdensity members based on the KS and AD tests. The percentiles for the mass-matched members of the field also appear in Table 1.

Table 1. A Summary of Percentiles for Some of the Prospector-inferred Physical Parameters

Parameter	Field			Overdensity			Mass-matched Field
	16th	50th	84th	16th	50th	84th	16th	50th	84th
AV /mag	0.06	0.16	0.53	0.08	0.23	0.48	0.07	0.20	0.74
${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })$	7.85	8.44	9.13	8.28	8.87	9.35	8.29	8.84	9.33
${\mathrm{log}}_{10}({\mathrm{SFR}}_{0\mbox{--}100\,{Myr}}/\left[{M}_{\odot }\,{\mathrm{yr}}^{-1}\right])$	−0.24	+0.15	+1.05	+0.18	+0.55	+1.12	+0.05	+0.57	+1.05
${\mathrm{log}}_{10}\left({\mathrm{SFR}}_{30\mbox{--}100\,{Myr}}/{\mathrm{SFR}}_{0\mbox{--}30\,{Myr}}\right)$	−1.44	−0.44	+0.11	−0.49	−0.09	+0.52	−1.23	−0.40	+0.11

Download table as:  ASCIITypeset image

When compared to the mass-matched field members, the overdensity members have similar inferred SFRs, consistent with being drawn from the same parent distribution based on the KS and AD tests. The ratio of SFRs in the two most recent time bins (corresponding to lookback times of 0–30 Myr and 30–100 Myr, respectively) also shows an interesting trend, where overdensity members have SFHs consistent with constant, while mass-matched field members have SFHs consistent with rising. Furthermore, the KS and AD tests suggest that the SFHs of the mass-matched field and overdensity sample exhibit statistically significant differences at roughly the 2σ–3σ level (corresponding to p-values of 0.003 < p < 0.05). These results suggest that the physical processes associated with this extreme galaxy overdensity at z = 5.4 have induced earlier star formation and earlier stellar mass assembly relative to the field, although there are large uncertainties associated with all of these parameters, and the distributions of these parameters appear consistent within the associated uncertainties.

To further understand the dynamical state of the overdensity identified in Section 3.1, the right panel of Figure 3 shows the phase-space diagram for the N = 39 confirmed members of the overdensity. Velocity offsets are calculated from the median redshift of the overdensity (see Figure 1), while spatial offsets are calculated from the median on-sky positions of the two groups that make up the bimodality (see the left panel of Figure 3). The median velocity offset of the overdensity is given by the gray dashed line, while the median velocity offsets of the first and second groups are given by the green and blue dotted lines, respectively. For each of the groups, the 1σ confidence intervals of the velocity offsets are given by the shaded regions (for group one, ${\rm{\Delta }}v=-{340}_{-50}^{+140}\,\mathrm{km}\,{{\rm{s}}}^{-1};$ for group two, ${\rm{\Delta }}v=+{190}_{-150}^{+120}\,\mathrm{km}\,{{\rm{s}}}^{-1}$).

Following the methodology described in Long et al. (2020), we derive two different estimates of the dark matter halo mass for the overdensity at z = 5.4. For both of these methods, we use the stellar-to-halo abundance matching relation described in Behroozi et al. (2013) to convert stellar masses into halo masses. Uncertainties on our estimates are calculated by adopting the mean stellar-to-halo abundance matching relation from Behroozi et al. (2013) and propagating the uncertainties on the stellar masses for each individual galaxy. This is an underestimate of the true uncertainties, since there is scatter in the stellar-to-halo abundance matching relation that we do not account for.

The first halo mass estimate is calculated by summing the halo masses for each individual galaxy that is part of the overdensity. This estimate assumes that: (1) each galaxy is formed in a separate dark matter halo; and (2) the individual galaxies are closer to virialization than the large-scale structure associated with the overdensity. We find that the first method gives a total halo mass of ${\mathrm{log}}_{10}\left({M}_{\mathrm{halo}}/{M}_{\odot }\right)=12.8\pm 0.1$. The second halo mass estimate is determined by summing the stellar masses for each individual galaxy that is part of the overdensity and converting to a halo mass. This estimate instead assumes that each galaxy is formed in the same dark matter halo. We find that the second method gives a total halo mass of ${\mathrm{log}}_{10}\left({M}_{\mathrm{halo}}/{M}_{\odot }\right)=12.6\pm 0.3$.

Weighing this kind of large-scale structure in the early Universe is challenging and requires a variety of assumptions. Typical methods for weighing galaxy clusters include: (1) gravitational lensing; (2) the Sunyaev–Zeldovich effect; and (3) using X-ray observations of the hot gas in the intracluster medium (ICM). However, current X-ray and submillimeter observations are not sensitive enough to properly weigh this extreme overdensity at z = 5.4. For this reason, we use the above dark matter halo mass estimates to infer a probable halo mass range of $12.6\lesssim {\mathrm{log}}_{10}\left({M}_{\mathrm{halo}}/{M}_{\odot }\right)\lesssim 12.8$ for the overdensity identified in Section 3.1.

We note that none of these methods account for additional members of the overdensity that were not identified and included in the final spectroscopic sample, including objects that fall outside either the JADES or the FRESCO footprints. This is a non-negligible effect, since the first component of the overdensity (see Figure 3) falls right at the edge of the JADES footprint ($\theta \ll 1^{\prime}$). Additionally, since our final spectroscopic sample only includes galaxies with narrow photometric redshift constraints and secure Hα-line detections, we are likely missing some additional subset of objects with relatively unconstrained photometric redshifts and/or low levels of star formation (e.g., DSFGs and/or obscured AGNs). There are zero known DSFGs at z = 5.4 in GOODS-S in the literature (e.g., Franco et al. 2018; Hatsukade et al. 2018; González-López et al. 2020; Gómez-Guijarro et al. 2022). However, wide and shallow surveys like GOODS-ALMA would miss DSFGs with infrared luminosities of L_IR ≲ 3 × 10¹¹ L_⊙ at z = 5.4. Given that clusters induce earlier quenching for their constituent members, we cannot rule out the existence of a significant population of these kinds of objects. For these reasons, the halo mass range quoted above is likely an underestimate of the true halo mass associated with this extreme galaxy overdensity at z = 5.4.

Figure 6 shows the dark matter halo mass distribution for some representative groups and protoclusters at z < 8. The two halo mass estimates derived here in Section 3 are given by the turquoise triangles, where the triangles indicate that these are likely underestimates of the true halo mass. Groups and protoclusters at z < 6 from Li et al. (2022) are given in grayscale, while protocluster candidates in the COSMOS field at z > 6 from Brinch et al. (2023) are given by the gray points for comparison, both selected based on photometric redshifts, with dark matter halo mass estimates that assume the same stellar-to-halo abundance matching relation used here. The magenta shaded region shows the expected halo mass evolution of a Coma-like cluster (Chiang et al. 2013) assuming a smooth evolution at z > 6. The black dashed line represents the typical threshold for shock stability assuming a spherical infall, below which the flows are predominantly cold and above which a shock-heated ICM is expected (Dekel & Birnboim 2006). The black diagonal dashed line represents the typical threshold for penetrating cold gas flows. The overdensity identified in Section 3.1 is expected to eventually evolve into a Coma-like cluster with ${\mathrm{log}}_{10}({M}_{\mathrm{halo}}/{M}_{\odot })\gt 15$ by z = 0, rivaling some of the most massive galaxy clusters found in the local Universe (Ruel et al. 2014; Buddendiek et al. 2015). For a Coma-like cluster at z = 0, the effective radius at z = 5.4 would be R_e ≈ 10 cMpc (Chiang et al. 2013). This value is much larger than the physical size of this extreme galaxy overdensity (see Figure 3). Thus, we conclude that the classification of this large-scale structure as a protocluster is justified.

Zoom In Zoom Out Reset image size

A number of previous works have found overdensities at similar redshifts to the one reported here. The earliest one identified at z > 5 is from Ouchi et al. (2005), who used narrowband imaging with the Subaru Suprime-Cam to identify galaxies with strong Lyα emission lines at z = 5.7 (Δz = 0.05). The distribution of sources was described as "clumpy," with one prominent overdensity significant at the 4.8σ level and a second one at the 2.2σ level. Follow-up spectroscopy confirmed the narrow redshift range for these clumps (Δz = 0.03), each being about 2' in diameter with a separation of about 9'. The narrow redshift range combined with the small on-sky separations suggests that these two overdensities may be clumps within a large-scale structure analogous to the one reported here.

Similar structures have been reported more recently (e.g., Jiang et al. 2018; Chanchaiworawit et al. 2019; Harikane et al. 2019; Wang et al. 2021), with similar techniques in the first three cases, and based on overdensities in the submillimeter found with the South Pole Telescope in the fourth. Jiang et al. (2018) found two overdensities near that identified by Ouchi et al. (2005), one at z = 5.68 with a diameter of $\sim 15^{\prime}$ and a smaller one at z = 5.75. These both are presumably related to the overdensities reported by Ouchi et al. (2005), given their proximity in redshift and on the sky. Chanchaiworawit et al. (2019) and Harikane et al. (2019) identified two similar structures at z = 6.5. Since the search method for these works starts with narrowband imaging, the samples are from a narrow redshift range and the actual density of such overdensities on the sky must be significantly higher than implied by the existing detections. Thus, massive clusters of galaxies must be well on the way to formation by z = 6, which is something that we expect based on simulations (Chiang et al. 2013, 2017) and more recent observations (Laporte et al. 2022; Morishita et al. 2023).

One of the key breakthroughs illustrated here with respect to these previous works is the efficient spectroscopic confirmation of such a large sample of galaxies at z = 5.2–5.5, made possible by the powerful combination of deep imaging and WFSS provided by JWST/NIRCam. As also presented in Kashino et al. (2023), JWST/NIRCam WFSS has enabled the discovery of three galaxy overdensities along the sightline of quasar J0100 + 2802 at z = 6.19, z_quasar = 6.33, and z = 6.78 through the blind detections of [O iii]-emitting galaxies. It is also worth mentioning that a similar z = 5.2 galaxy overdensity was identified in the GOODS-N field with a probable halo mass range of $12.3\,\lt {\mathrm{log}}_{10}\left({M}_{\mathrm{halo}}/{M}_{\odot }\right)\lt 12.9$ through Lyα and submillimeter spectroscopy (Walter et al. 2012; Calvi et al. 2021), which was also partially observed by FRESCO in 2023 February (see the recent papers by Herard-Demanche et al. 2023; Sun et al. 2024). Future observations with JWST will certainly find more structures similar to those highlighted here, allowing a more complete look at the progenitors of the most massive gravitationally bound structures in the local Universe: galaxy clusters.