Metallicity Scatter Originating from Sub-kiloparsec Starbursting Clumps in the Core of a Protocluster at z=7.88

Our NIRSpec/IFU observations (GTO4553, PI Stiavelli) were executed within a single visit on Oct 27, 2024, centered at (R.A., Decl.) = ($3.60680,-30.38082$). The visit was configured at the V3 position angle PA_V3=41.1 degree with the high-resolution grating G395H/F290LP. The exposure consists of 12 dithers, with the total science time 19,432 sec (7.44 hrs including overhead was charged).

The data are reduced using the official JWST pipeline (ver 1.16.0) with additional custom steps. Briefly, we retrieve the level1 products (_uncal.fits) from MAST. We run the stage1 step and then the stage2 step, with additional stripe elimination using NSClean (Rigby et al., 2023)¹¹1https://github.com/JWST-Templates/NSClean and cosmis ray rejection using lacosmic (van Dokkum, 2001; Bradley, 2023). We also run the source detection tool SExtractor (Bertin & Arnouts, 1996) in each _rate image and flag significant point-like source, i.e., cosmic ray (CR) that are not detected in the previous steps. We run the stage3 step to construct a cube and subtract background in each wavelength frame using off-source regions. We run lacosmic again in each wavelength frame to mask any residual CRs.

We utilize the JWST NIRCam imaging data available in the field. Those include; GLASS-JWST (ERS1324; Treu et al., 2022), UNCOVER (GO2561; Bezanson et al., 2022), DD2756 (Chen et al., 2022), GO2883 (PI F. Sun), All the Little Things (ALT; GO3516; Naidu et al., 2024), GO3538 (PI E. Iani), BEACON (GO3990; Morishita et al., 2024b), and Medium Bands, Mega Science (GO4111; Suess et al., 2024). We also utilize the Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) and Wide Field Camera 3 (WFC3)-IR mosaic images publicly available by the GLASS-JWST team (Merlin et al., 2022), which combined data from multiple HST programs (Postman et al., 2012; Lotz et al., 2017; Kelly et al., 2018; Steinhardt et al., 2020).

All data were compiled and used for photometry as described in Morishita et al. (2024c). Briefly, we identified sources in the IR detection image (F277W+F356W+F444W stack) using SExtractor (Bertin & Arnouts, 1996), and measured fluxes on the mosaic images (matched to the F444W point-spread function, PSF), with a fixed aperture of radius $r=0.\!^{\prime\prime}16$ . The measured fluxes are calibrated using the method introduced in Morishita & Stiavelli (2023). Interested readers are referred to the work for more details. The filters included in our analysis and the limiting magnitudes ($5\,\sigma$, for a point source) are: HST-F435W (28.7), F606W (28.1), F775W (28.3), F814W (27.7), F105W (27.5), F125W (27.5), F140W (28.2), F160W (27.1), NIRCam-F070W (28.5), F090W (28.9), F115W (28.7), F140M (27.9), F150W (28.8), F162M (27.9), F182M (28.1), F200W (28.9), F210M (28.0), F250M (28.0), F277W (29.1), F300M (28.1), F335M (28.1), F356W (29.1), F360M (28.2), F410M (28.6), F430M (27.7), F444W (28.9), F460M (26.8), F480M (27.2), NIRISS-F115W (28.8), F150W (28.7), F200W (28.6), F356W (28.8), F430M (27.6), F444W (28.5), F480M (27.0).

In the following analysis, we adopt the latest magnification model by Bergamini et al. (2023a, b).

Emission lines are extracted using the reduced flux cube. Given the limited Signal-to-Noise ratio (SNR), we mainly base our analysis on the strongest emission line, $[\textrm{O}~{}\textsc{iii}]_{\lambda 5007}$. We extract the $[\textrm{O}~{}\textsc{iii}]_{\lambda 5007}$ emission line map by coadding cube frames (each weighted by the corresponding rms frame) over the velocity range $\pm 300$ km s^-1 centered at $z=7.88$, i.e., the systemic redshift of ZD3 and ZD6.

In Fig. 1, we show the extracted $[\textrm{O}~{}\textsc{iii}]_{\lambda 5007}$ emission line map along with the NIRCam images. We identify five separate regions detected with SNR $>3$. Two are ZD3 and ZD6, previously spectroscopically confirmed galaxies (Morishita et al., 2023). In addition, one emitting region is detected east of ZD6, hereafter referred to as ZD6-E. Two strong [O iii]-emitting regions are identified within the eastern system ZD12²²2Zheng et al. (2014) presented 11 $z$-band dropout sources. ZD12 was not identified in the original photometric study, possibly due to confusion with the foreground galaxy., which we refer to as ZD12-W and ZD12-E, respectively. In what follows, we base our analysis on those five regions. We refer to them as [O iii]-emitting regions, to distinguish them from UV clumps that are defined separately later in this paper.

Spectra are extracted within apertures of $\sim 0.\!^{\prime\prime}15$–$0.\!^{\prime\prime}3$ centered in the [O iii]-emitting region, where the aperture size is determined by a growth curve method. Furthermore, we define larger segments for the entire ZD3 + ZD6, ZD12, and ZD3 + ZD6 + ZD12 systems (regions refined by yellow lines in Fig. 1) and extract the spectra for further analysis. The extracted spectrum of the total system is shown in Fig. 2.

We note a foreground galaxy (ID 51545, 00:14:25.7, $-$30:22:50.6), to the South-East of ZD12 (Fig. 1). This galaxy is clearly detected in the F090W image and thus is not a $z$-band dropout source; instead, the galaxy has a photometric redshift solution of $z\sim 1.6$ and is best fitted with a star-forming galaxy template ($\log M_{*}/M_{\odot}\sim 8.5$), while our IFU cube does not reveal strong emission lines.

We fit each emission line of interest with a Gaussian function. Our basic strategy is to define a wavelength window for each line, model the underlying continuum spectrum (either source origin or residual background) by a polynomial fit, and subtract it from the one-dimensional spectrum before the Gaussian fit is carried out. The model includes amplitude, line width, and redshift parameters. The [O iii]_{λλ4959+5007} doublet lines are modeled with a fixed line ratio (1:3), with common line width and redshift parameters. The [O ii]_{λλ3726+3729} doublet lines are modeled with common line width, but without their line ratio fixed (see Sec. 3.2). For emission lines other than those two doublet pairs, we fix the redshift parameter to the value determined by the [O iii]-doublets and leave the amplitude and width as free parameters.

The flux of each line is measured by integrating the modeled Gaussian component, and the flux error is estimated by summing the error weighted by the amplitude of the Gaussian model in quadrature. The uncertainty associated with the modeled continuum is also added. In the following analysis, we adopt flux measurements when the line is detected (signal-to-noise ratio SNR $\geq 1.5$); for those not detected, we proceed our analysis with the 1.5-$\sigma$ upper limit. We securely detect the [O iii] doublets in all regions. ${\rm H\beta}$ in four and [O ii] in three. In addition, the auroral $[\textrm{O}~{}\textsc{iii}]_{\lambda 4363}$ is detected in ZD12-W, allowing us to measure the electron temperature (Sec. 3.2).

We find that all five regions are within a small distance along the radial direction ($\Delta v=1500$ km s^-1). The separation between ZD3 and ZD6 is even smaller, $\delta v=350$ km s^-1, confirming that they are interacting, rather than by chance alignment in the sky plane.

The measured line fluxes are corrected for dust attenuation. Although it is ideal to determine the correction on the basis of the ratio of Balmer lines, no more than one line is confidently detected in our cube, except for ZD12-W. We thus apply the correction derived from our resolved SED analysis (Sec. 3.3), after rescaling it by a factor of 2.27 (Calzetti et al., 2000, also Morishita et al. 2024d, which found the assumption is valid using $z>3$ galaxies) that accounts for continuum and emission line dust effects. Dust-corrected line fluxes are reported in Table 1.

We infer spectral energy distributions (SEDs) within the detected region of our sample galaxies using the JWST+HST photometric data. Firstly, the spatially resolved segments are defined by the voronoi binning method (Cappellari & Copin, 2003), by requiring a minimum SNR of 5 in the F444W image. Fluxes and flux errors in each region are estimated using the PSF-matched images. We use the SED fitting code gsf (ver1.9; Morishita et al., 2019), which allows flexible determinations of the SED by adopting binned star formation histories (often referred to as non-parametric), and adopt the SMC dust curve (Gordon et al., 2003). We follow the same procedure as described in Morishita et al. (2024d).

Physical quantities of each component, e.g., stellar mass and UV luminosity, are calculated by integrating the derived quantity within the defined region. The UV luminosity is corrected for dust attenuation using the $\beta_{\rm UV}$ slope, which is measured by using the posterior SED, as in Smit et al. (2016):

The attenuation corrected UV luminosity is then converted to SFR via the relation in Kennicutt (1998):

Since the relation is essentially for the Salpeter IMF, the estimated SFR is corrected for the Chabrier IMF, by multiplying a factor of 0.63, the factor estimated from fsps.

By using the ${\rm H\beta}$ flux in Sec. 3.1 and the UV luminosity derived here, we can measure the ionization photon production efficiency in each [O iii]-emitting region (e.g., Schaerer et al., 2016; Prieto-Lyon et al., 2023):

where $N_{\rm LyC}$ is the total ionizing photons of Lyman continuum and $L_{\rm UV}$ intrinsic UV luminosity density measured at rest-frame 1500 Å  both corrected for attenuation. Since the luminosity in the optical hydrogen recombination lines is proportional to the number of Lyman-continuum photons absorbed in the galaxy, the direct measurement of $N_{\rm LyC}$ is proportional to the ionizing photon escape fraction, $f_{\rm esc}$:

where we here adopt $f_{\rm esc}=0$.

Similarly, using the best-fit SED, we estimate rest-frame equivalent widths (EW) of [O iii]-doublet emissions. We note that continuum flux is not detected in the IFU cube in any of the [O iii]-emitting regions. This means that our EW measurement is rather indirect and susceptible to systematics such as absolute flux calibration in the cube. However, we have validated the absolute flux accuracy by comparing the measured line fluxes with the flux excess seen in two broadband photometry, F410M and F444W (where F444W captures the ${\rm H\beta}$+[O iii] and F410M does only continuum). The uncertainty in the continuum flux level, estimated in the SED modeling, is also added to the EW measurement uncertainties. The measurements are reported in Table 2.

We have seen in Sec. 3.2 that the measured metallicity of the [O iii]-emitting regions distributes a wide range, from $\log{\rm O/H}+12<7$ to $\sim 8$ (Fig. 4). Although the sample number here is small, the large scatter observed is qualitatively consistent with previous studies using non-IFU spectra, reporting $\sim 0.3$ dex (e.g., Morishita et al., 2024d). Morishita et al. (2024d) attributed it partially to the elevated star-formation activity in early galaxies (see below; also Sec. 1) but also to the aperture effect by the MSA slit. Our IFU observations confirm that large variation of metallicity is present even within in single galaxy systems in these early times. This qualitatively supports the previous interpretation that observed scatters around the MZ relation are partially due to the limited aperture size of MSA.

Notably, the total metallicities measured in the integrated spectra are more represented by enriched components. The integrated ZD3+ZD6 system is at $\log{\rm(O/H)}$ $+12\sim 7.9$, whereas ZD3 is characterized with a much lower metallicity, $\log{\rm(O/H)}$ $+12\sim 7.3$. Similarly, the integrated ZD12 system is characterized by a higher metallicity, $\log{\rm(O/H)}$ $+12\sim 7.3$, than ZD12-E alone, $\log{\rm(O/H)}$ $+12\,\hbox to0.0pt{\lower 3.5pt\hbox{$\mathchar 0\relax\sim$}\hss}\raise 1.0pt\hbox{$<$}\,6.9$. This is a natural consequence of the MZ relation that more massive regions (i.e., more luminous) outshine less enriched (fainter) regions, particularly in young galaxy systems. Interestingly, the region with the least enrichment, ZD12-E with an upper limit $\log{\rm(O/H)}$ $+12\,\hbox to0.0pt{\lower 3.5pt\hbox{$\mathchar 0\relax\sim$}\hss}\raise 1.0pt\hbox{$<$}\,7$, is offset for $\sim 1$ kpc from the brightest region of ZD12 in NIRCam imaging. Since MSA masks are configured to source light-weighted positions by default, this presents a challenge in current MSA mask configurations in finding pristine regions of galaxies.

Although this type of observational bias is known for low-$z$ observations in the presence of mild metallicity gradients (e.g., Tremonti et al., 2004; Belfiore et al., 2017), the impact is likely more significant at high redshift because of the time available for chemical mixing within ISM. Morishita et al. (2024d) reported the redshift evolution of scatter amplitude, from $\sim 0.1$ at $z\sim 0$ (Curti et al., 2020), to $\sim 0.2$ at $\sim 0.2$ dex at $3<z<5$ and $\sim 0.3$ dex at $z>5$ (also Curti et al., 2023). The impact of star formation feedback, especially in low-mass systems of shallow potential, is therefore expected to be more significant at higher redshifts. Pallottini et al. (2024) show that high stochasticity in high-$z$ star formation (also Shen et al., 2023; Mirocha & Furlanetto, 2023; Muñoz et al., 2024) enables the system to have large variations of metallicity, making the mass–metallicity relationship less constrained (also see Sec. 4.2).

On the other hand, the metallicity of ZD6 is observed considerably high, $\log{\rm(O/H)}$ $+12\sim 8$, or $10\,\%$ solar. This, given the age of the universe $\sim 660$ Myrs, suggests that metal enrichment within some sub-galactic regions could be much more rapid than other regions in the same system. In particular, ZD6 is located near the center of an ongoing merger, where an accelerated gas cycle is expected. Our SED analysis of resolved pixels has revealed an increase in $A_{V}$, which is another evidence of fast enrichment.

As redshift increases, so does the merger rate (Fakhouri & Ma, 2008; Rodriguez-Gomez et al., 2015), and thus galaxy assembly is progressively driven by mergers rather than smooth accretion (see also Hsiao et al., 2023; Boyett et al., 2024; Dalmasso et al., 2024; Sugahara et al., 2024, for recent observations). As a result, the chemical enrichment of the system is significantly affected by individual merger events (and also the properties of independent building blocks), rather than by smooth accretion of material and self-recycling, increasing the metallicity variation. As such, in addition to the general increase in star formation intensity in the early universe, our study showcases that an increased fraction of merger-driven evolution would also contribute to the observed scatter in the mass–metallicity relation. A large data set of IFU observations may shed light on the relative contribution of each evolutionary channel.

NIRSpec IFU is limited to the spatial resolution of $\sim 0.\!^{\prime\prime}15$ ($\sim 370$ pc with $\mu\sim 2$ magnification), leaving the nature of the [O iii]-emitting regions rather undetermined. In the NIRCam images (Fig. 1), which offer $\sim 5\times$ higher resolution, we see multiple UV-bright clumps in both integrated ZD12 and ZD3+ZD6 systems, at a smaller scale than the [O iii] regions.

To study in detail these UV-bright clumps, we fit the two-dimensional light distribution in the F150W image, which corresponds to $\sim 1700$ Å rest-frame, by using galfit (Peng et al., 2002, 2010). Each light profile of 10 UV clumps (and one foreground galaxy) is well fitted by a single-Sérsic component (Fig. 5). Four (# 0, 1, 4, 10) have effective radii smaller than the PSF size ($\,\hbox to0.0pt{\lower 3.5pt\hbox{$\mathchar 0\relax\sim$}\hss}\raise 1.0pt\hbox{$<$}\,2$ pixels, or $R_{e}\,\hbox to0.0pt{\lower 3.5pt\hbox{$\mathchar 0\relax\sim$}\hss}\raise 1.0pt\hbox{$<$}\,0.\!^{\prime\prime}06$). We consider those clumps unresolved and quote the PSF’s FHWM as their upper limit sizes.

In Fig. 6, we show the size distribution of the UV clumps as a function of the stellar mass. The stellar mass of each clump is estimated by using the mass-to-light ratio derived in Sec. 3.3 and scaling it with the total flux modeled with galfit. The four unresolved UV clumps are found in the regime where relatively massive clusters in the local universe are located (Norris et al., 2014). However, with the spatial resolution available here, the detailed classification of those unresolved clumps is undetermined. The other six UV clumps (# 2, 3, 5, 6, 8, 9) have a size comparable to that of individual galaxies (Morishita et al., 2024a).

Those unresolved clumps are characterized with a high star formation surface density $\Sigma_{\rm SFR}$ ($>10$–$100\,M_{\odot}/{\rm yr/kpc^{2}}$) and a surface mass surface density $\Sigma_{\rm*}$ ($>500$–$1600\,M_{\odot}/{\rm pc^{2}}$; Fig. 7). The observed densities are on the upper end of individual galaxy populations at $z>5$ (e.g., Morishita et al., 2024a), even reaching the realm of star clusters.

Recent observations have seen some evidence of star clusters forming in the early environment through gravitational lensing (e.g., Vanzella et al., 2022, 2023a; Meštrić et al., 2022; Adamo et al., 2024; Messa et al., 2024a; Fujimoto et al., 2024). For example, Messa et al. (2024b) recently presented NIRSpec/IFU observations and revealed star clusters of $R_{e}\sim 3$–10 pc within a strongly lensed galaxy system at $z=6.14$. These clusters are strongly line-emitting and are characterized by blue UV slopes, $\beta_{\rm UV}\,\hbox to0.0pt{\lower 3.5pt\hbox{$\mathchar 0\relax\sim$}\hss}\raise 1.0pt\hbox{$<$}\,-2.5$ and $\,\hbox to0.0pt{\lower 3.5pt\hbox{$\mathchar 0\relax\sim$}\hss}\raise 1.0pt\hbox{$<$}\,0.1\,Z_{\odot}$.

In this context, the nature of star formation in our unresolved clumps is of particular interest. The UV clump #0 and 1 are both located within the [O iii]-emitting region (ZD12-W), one of the metal-poor regions ($\log{\rm(O/H)}$ $+12=7.4$ from the direct-T method) with high [O iii]-equivalent width, EW${}_{0}=670\pm 60$ Å. Remarkably, the measured ionizing photon production efficiency (Sec. 3.3) for ZD12-W is considerably high with $\log\xi_{\rm ion}=25.7$ Hz/erg for its $M_{\rm UV}$ when compared to those of the literature (e.g., Vanzella et al., 2024). Using the linear relation between $\xi_{\rm ion}$ and $M_{\rm UV}$ derived by Prieto-Lyon et al. (2023), a rather moderate efficiency, $\sim 25.4$, would be expected for ZD12-W’s $M_{\rm UV}$. This suggests a higher escape fraction in ZD12-W than assumed here, i.e., $f_{\rm esc}\sim 0.5$ would be required to meet the expected $\xi_{\rm ion}$. Alternatively, if the observed ${\rm H\beta}$ originates mainly in one of the two UV clumps (which have $M_{\rm UV}\sim-17.7$ and $-16.9$ mag), it would make the deviation smaller, but still not completely gone.

Star formation in such dense environments plays a crucial role in chemical enrichment and may be responsible for the diverse properties observed among our samples. Recent JWST studies have also illuminated the origin of star clusters by identifying Nitrogen-rich ISM in the early universe (e.g., Topping et al., 2024; Stiavelli et al., 2024). A detailed study of the chemical compositions within those UV clumps would thus give us further insight into the nature, e.g., whether they are giant star-forming clumps (similar to those at lower redshifts) or even smaller-scale populations such as star clusters.

Theoretical studies investigated the nature of star-forming clumps seen at high redshifts in relation with metallicity. Sugimura et al. (2024) show that high gas density can be reached within an environment of metal-poor star formation because of the intense far-UV radiation field therein. This is qualitatively consistent with the nature of ZD12-W — low metallicity, high electron density, and high surface densities. Nakazato et al. (2024) show that short-lived, UV-bright clumps can form through merger, which may also be relevant to the case of the ZD3+ZD6+ZD12 system. Similarly, Marconcini et al. (2024a) investigated the moderate anti-correlation between $\Sigma_{\rm SFR}$ and $\log{\rm(O/H)}$ seen in a galaxy system CR7. These authors attributed it to the inflow of pristine gas, which would lower the metallicity and trigger intense star formation, though the system’s global metallicity is much more enriched ($\log{\rm(O/H)}$ $+12\sim 8$) than our case here. Different mechanisms responsible for high $\Sigma_{\rm SFR}$ may be in action at different metallicity regimes.

Lastly, we revisit the characterization of A2744-z7p9 as a protocluster system. A2744-z7p9 was reported in Morishita et al. (2023) for the spectroscopic confirmation of seven member galaxies at $z=7.88$ (also Roberts-Borsani et al., 2022). The following studies have identified two (YD1, s1; Hashimoto et al., 2023) and one (ZD4; Chen et al., 2024) additional members. With the new addition of ZD12 in this study, now 11 member galaxies are confirmed (Fig. 1), distributed within a small projected region ($r\sim 60$ pkpc, in the source plane) and a redshift slice ($\delta_{z}=0.019$, or $\sim 25$ pMpc). Following Morishita et al. (2023), we here aim to estimate the overdensity factor, the total halo mass, and the system velocity dispersion.

The overdensity factor, defined as $\delta=(n-\bar{n})/\bar{n}$, represents the excess of the surface number density from the field average. We use $R=60$ pkpc, which defines the extent of the member galaxies in the source plane, and the redshift slice of $\delta_{z}=0.019$. One of the member galaxies, s1, is excluded from the calculation, as this source is rather considered as a sub-region and does not qualify for “galaxy” count. For field reference, we use the luminosity function at $z\sim 8$ presented in Bouwens et al. (2021) and integrate it down to $M_{\rm UV}\sim-17.4$ mag, the lowest magnitude among the confirmed member galaxies (except for s1). For the area within $R=60$ pkpc, we expect $\bar{n}=0.2_{-0.4}^{+0.5}$ galaxies in average fields. With the reference number for the average field, we find $\delta=44_{-31}^{+89}$ for A2744-z7p9, where the uncertainties quoted represent Poisson noise.

To estimate the total halo mass of A2744-z7p9, we estimate the halo mass of the individual components from the halo-mass galaxy-luminosity relation (Mason et al., 2023). Using the relation and adding the halo mass estimates of all confirmed members, we obtain the total halo mass estimate $5.8_{-0.3}^{+0.2}\times 10^{11}\,M_{\odot}$, but this is likely a lower limit given that there are candidate selected photometrically that are not spectroscopically followed.

In addition, we can take advantage of the spectroscopic data to obtain an estimate of the velocity dispersion for the whole system. We adopt a simple Gaussian estimator and bootstrap method to derive the uncertainty (Beers et al., 1990), obtaining $1100\pm 500\,{\rm km~{}s^{-1}}$, consistent with the original estimate in Morishita et al. (2023). We caution the reader that 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.

We note that the only ${\rm Ly\alpha}$-emitting member galaxy, ZD4 (Chen et al., 2024), is located in an off-core region of the protocluster (Fig. 1).