A proto-cluster of massive quiescent galaxies at z=4

We report on discovery of a concentration of massive quiescent galaxies located at z=4. The concentration is first identified using high-quality photometric redshifts based on deep, mutli-band data in Subaru/XMM-Newton Deep Field. Follow-up near-infrared spectroscopic observations with MOSFIRE on Keck confirm a massive (~10^{11} Msun) quiescent galaxy at z=3.99. Our spectral energy distribution (SED) analyses reveal that the galaxy experienced an episode of starburst about 500 Myr prior to the observed epoch, followed by rapid quenching. As its spectrum is sufficiently good to measure the stellar velocity dispersion, we infer its dynamical mass and find that it is consistent with its stellar mass. The galaxy is surrounded by 4 massive (>10^{10} Msun) quiescent galaxies on a ~1 physical Mpc scale, all of which are consistent with being located at the same redshift based on high-accuracy spectro-photometric redshifts. This is likely a (proto-)cluster dominated by quiescent galaxies, the first of the kind reported at such a high redshift as z=4. Interestingly, it is in a large-scale structure revealed by spectroscopic redshifts from VANDELS. Furthermore, it exhibits the red sequence, adding further support to the physical concentration of the galaxies. We find no such concentration in the Illustris-TNG300 simulation; it may be that the cluster is such a rare system that the simulation box is not sufficiently large to reproduce it. The total halo mass of the quiescent galaxies is ~10^{13} Msun, suggesting that they form a group-sized halo once they collapse together. We discuss implications of our findings for the quenching physics and conclude with future prospects.

1. INTRODUCTION Galaxies in the local Universe exhibit a large diversity in their morphologies as characterized by the Hubble Sequence.Most massive galaxies tend to be elliptical and they show little sign of on-going star formation activities.They form a tight sequence on a colormagnitude diagram (red sequence; Bower et al. 1992) and the tightness indicates that they are a very homogeneous population.Their spectra indeed show that they are all dominated by old stars formed long time ago.Furthermore, detailed absorption line analyses suggest that they formed in an intense burst of star formation occurred in a short timescale, followed by passive evolution (i.e., simple aging without formation of new stars; Thomas et al. 2010).This is an intriguing observational fact because these galaxies are massive, and we expect massive galaxies to appear only at low redshifts in the current paradigm of hierarchical galaxy formation.These observations pose three major questions: • What physical process(es) drives the intense starburst in the early Universe?
• What causes the rapid shutdown of the starburst activity?
• What physical process(es) prevents subsequent formation of new stars for >10 Gyr?
These questions are indeed among the most outstanding questions in the field of galaxy formation and evolution for decades.This paper aims to address the first two questions.
An interesting observational approach to these questions is to go back in time and study progenitors of nearby massive elliptical galaxies in the early Universe.A lot of efforts have been put in the search for distant quiescent galaxies (i.e., galaxies with no star formation) and quiescent galaxies have been identified up to z ∼ 4 and beyond (Schreiber et al. 2018a;Tanaka et al. 2019;Valentino et al. 2020;Kakimoto et al. 2023) thanks to the advent of sensitive near-infrared spectrographs on the ground.These galaxies are post-starburst galaxies; their star formation histories (SFHs) inferred from extensive fitting of their SEDs with stellar population synthesis models indicate that they all underwent an intense starburst in a recent past and then experienced rapid quenching.We are getting very close to the formation epoch of the nearby elliptical galaxies.However, as we show below, z ∼ 4 is the spectroscopic sensitivity limit from the ground and it is extremely challenging to go further.
JWST has already started to identify quiescent galaxies at z 3 from its very first observations (Carnall et al. 2023a,b;Marchesini et al. 2023;Valentino et al. 2023), and it will surely revolutionize our understanding of the formation of massive quiescent galaxies.Due to its small field coverage, however, one key aspect will be missed in blank field JWST observations: environment.In the local Universe, we know galaxy properties are strongly dependent on their surrounding environment (Dressler 1980).Quiescent earlytype galaxies dominate the dense cluster core, while starforming spiral galaxies are the main population in the low-density field.Also, we expect that galaxies form earlier in higher density peaks of the density fluctuation in the Universe; galaxy formation is accelerated in cluster regions.However, massive clusters are rare objects in the Universe especially at high redshifts.Thus, a wide-field observation is required to understand the role of environment in the formation of massive quiescent galaxies in the distant Universe.
There have been a lot of efforts to search for (proto-)clusters using star-forming galaxies as a tracer (e.g., Miller et al. 2018;Oteo et al. 2018;Toshikawa et al. 2018).
However, much less attention has been paid to quiescent galaxies due primarily to difficulties in their identifications at high redshifts; they are faint and rare objects (∼ 10 −6 Mpc −3 ; Schreiber et al. 2018a;Valentino et al. 2020;Marsan et al. 2022;Valentino et al. 2023), requiring deep imaging data over a wide area.Deep multiwavelength data in fields like the Cosmic Evolution Survey (COSMOS) field make it possible to search for them in the early Universe.For example, Strazzullo et al. (2015) searched for concentrations of quiescent galaxies in the COSMOS field out to z ∼ 2.5.Recently, Ito et al. (2023a) successfully confirmed a concentration of massive quiescent galaxies with deep near-IR spectroscopy at z = 2.7.Kubo et al. (2015Kubo et al. ( , 2021) ) identified a massive quiescent galaxy in a proto-cluster at z = 3.1.The currently known high-z proto-clusters with quiescent galaxies are still only a few and are limited to z 3.In this paper, we focus on quiescent galaxies at z = 4 based on deep spectroscopic observations and discuss an early structure populated by these quiescent galaxies to address the role of environment.
The paper is structured as follows.We describe our spectroscopic observations in Section 2, followed by spectroscopic analyses of a confirmed quiescent galaxy at z = 3.99 in Section 3. We then focus on a concentration of several quiescent galaxies at the same redshift using joint photometric and spectroscopic fits to their observed spectral energy distributions (SEDs) in Section 4. Section 5 summarizes physical properties of the concentration such as the Butcher-Oemler effect (Butcher & Oemler 1984) at z = 4. Finally, we discuss our results in the context of massive galaxy formation and conclude the paper in Section 6.We adopt H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7.Magnitudes are in the AB system (Oke & Gunn 1983).

OBSERVATION
We focus on the Subaru/XMM-Newton Deep Survey (SDXS), where very deep, multi-wavelength imaging data are available over a wide area.We have constructed a multi-band photometric catalog as summarized in Kubo et al. (2018).This is the same catalog that Tanaka et al. (2019) and Valentino et al. (2020) are based on.We apply the photometric redshift code from Tanaka (2015) to the catalog to infer photometric redshifts as well as physical properties of galaxies such as stellar mass and star formation rate (SFR).In brief, the code uses models from Bruzual & Charlot (2003) over a wide parameter range with template error functions to account for systematic uncertainties in the models.We adopt the Chabrier (2003) initial mass function.
The code is calibrated to relatively low-redshift objects with spectroscopy available at that time.Since then, VANDELS (Pentericci et al. 2018;Garilli et al. 2021) has carried out an extensive spectroscopic campaign focused on high-redshift galaxies in the field and released their catalog to the public.We compare our photometric redshifts with secure spectroscopic redshifts from VANDELS DR4 with zflag= 3 or 4 in Fig. 1.We note that the VANDELS targets are selected on the basis of photometric redshifts.There are a few different types of targets, but we are primarily using quiescent galaxies at 1 < z < 2 and star-forming galaxies at 2.5 < z < 5 here.
We emphasize that Fig. 1 is entirely a blind comparison; VANDELS spec-z's were not available when we computed our photometric redshifts.As can be seen, the correspondence is excellent: an outlier rate defined as the fraction of objects with |δz|/(1 + z) > 0.15 (Tanaka et al. 2018) is 4.8% and the dispersion defined as σ((δz)/(1 + z)) is 0.036.The outliers scatter down to low redshifts, but most of them are consistent with spec-z's within the uncertainties.Photometric redshifts plotted against the secure spectroscopic redshifts from VANDELS DR4.The error bars indicate the 68% confidence intervals.The red circles are spectroscopically confirmed quiescent galaxies from Schreiber et al. (2018a), Tanaka et al. (2019), and this paper.We differentiate spectrophotometric galaxies (see Section 4) from the secure spectroscopic redshifts with the open triangles.The dashed line shows the z phot = zspec relation.
As VANDELS does not include quiescent galaxies at high redshifts, we supplement it with high-z quiescent galaxies from the literature (Schreiber et al. 2018a;Tanaka et al. 2019) as well as from this paper and show them as the red circles in Fig. 1.While the redshift distribution of the quiescent galaxies is very sparse, the photo-z accuracy for the quiescent galaxies at z ∼ 4 is similar to that of star-forming galaxies from VANDELS.The figure demonstrates that our photometric redshifts are accurate even for quiescent galaxies at high redshifts.
Motivated by the good accuracy, we search for an overdensity of massive quiescent galaxies at high redshifts.We search for over-density regions by selecting quiescent galaxies in a redshift slice of δz ± 0.3 and evaluate their over-density significance as a function of position.This redshift slice roughly corresponds to ±1.5σ(z phot ).Quiescent galaxies are defined as those with 1σ upper limit of specific SFR (sSFR) being smaller than 10 −9.5 Gyr −1 , which corresponds to ∼ 1 dex below the star formation main sequence, as we have done in our previous papers (Kubo et al. 2018;Valentino et al. 2020;Ito et al. 2022Ito et al. , 2023a)).We identify a significant over-density of quiescent galaxies at z ∼ 4 in SXDS as shown in Fig. 2; based on the kernel density estimate with Gaussian width of 2.4 arcmin (1 physical Mpc at z = 4), the over-density  The open red/blue circles are photo-z selected quiescent/star-forming galaxies at 3.7 < z phot < 4.3.The color contours show the density of quiescent galaxies smoothed over 1 Mpc (physical).The filled red points are galaxies in the over-density region and are the focus of the paper.There is a spectroscopically confirmed z = 4.01 galaxy in the north of the over-density region (Tanaka et al. 2019).The green square is the field field of view of our MOSFIRE observation.The bar on the bottom-left indicates 1 Mpc (physical) at z = 4.
significance is ∼ 20σ.Tanaka et al. (2019) reported the discovery of the z = 4.01 galaxy in the same field and at the same redshift as the over-density region.This is not a coincidence as we discuss below.The over-density region is approximately 5 arcmin (∼ 2 physical Mpc) on a side.

Observation and Data Reduction
We carried out a follow-up spectroscopic observation of the region with MOSFIRE (McLean et al. 2010) on the Keck I telescope on Nov 29 -Dec 2, 2020.As the over-density region is compact, we could observe all the quiescent targets in a single MOSFIRE field of view.We observed in the K-band and applied the ABBA dither along the 0.7 arcsec width slit with an individual exposure time of 3 min and a nod amplitude of ±1.25 arcsec.The conditions were clear with slightly variable seeing between 0.5 and 1.2 arcsec with a median of 0.74 arcsec.The data were reduced with the MOSFIRE DRP 2018 release with custom scripts (Kakimoto et al. 2023).We detail our data processing in Appendix A. The total integration was ∼ 16 hours.
We observed 5 quiescent galaxy candidates as summarized in Table 1.Only one of them is sufficiently bright to measure its redshift from absorption features and its spectrum is shown in Fig. 3.The spectrum is rather noisy (the median S/N per pixel is 1.8), but we can still identify multiple Balmer absorption lines, which yield z spec = 3.987 ± 0.001.The continuum is relatively flat overall, but there is a hint of the 4000 Å break on the blue end, which indicates an evolved stellar population.We will discuss this object in detail in the next Section.We could not measure redshifts for the reminder of the quiescent galaxy candidates observed in the same mask as they are fainter (see Table 1).We attempt to constrain their redshifts using both the spectra and photometry and discuss their properties in Sections 4 and 5.

DETAILED PROPERTIES OF THE CONFIRMED Z = 4 GALAXIES
In this Section, we examine properties of the z = 3.99 galaxy in detail and compare them with those of the z = 4.01 galaxy presented in Tanaka et al. (2019) and Valentino et al. (2020).We closely follow the analysis of Tanaka et al. (2019) in this Section.We measure its physical size (Section 3.1), and discuss the dynamical properties (Section 3.2).We then infer its star formation histories (Section 3.3).As we show below, the two galaxies are both massive galaxies with suppressed SFRs, but there are interesting differences.The z = 3.99 galaxy is in an over-density region (Section 5), while the z = 4.01 galaxy is not.We compare their properties with their environmental differences in mind.

Size-Mass Relation
We first focus on the size vs. stellar mass relation.Galaxies are known to exhibit a tight size-mass relation in the sense that more massive galaxies are physically larger (Shen et al. 2003).Interestingly, quiescent galaxies tend to be more compact than star-forming galaxies of similar stellar mass.This relation evolves with redshift; galaxies at fixed stellar mass are smaller at higher redshifts (e.g., van der Wel et al. 2016).The evolution has been tracked up to z ∼ 4 and beyond (Kubo et al. 2018;Ito et al. 2023a) and massive quiescent galaxies there are extremely compact and dense.We briefly address this size-mass relation at z ∼ 4.
We do not have high-resolution imaging data of the z = 3.99 galaxy, but as shown in Tanaka et al. (2019), the ground-based high-quality images deliver sizes consistent with those from Hubble Space Telescope.Furthermore, the AO-assisted K s -band imaging of the z =  Figure 3.The top and bottom panels show the 2d and 1d spectrum, respectively.In the bottom panel, the black and blue spectra are the observed and noise spectra binned by 3 pixels.Some of the most prominent absorption features are indicated.The red spectrum is the best-fit ppxf model (see Section 3.2 for details).
4.01 galaxy presented in Ito et al. (2023b) is consistent with that measured from the seeing-limited imaging data from the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP; Aihara et al. 2018Aihara et al. , 2022)).We thus adopt the same approach here to measure the size of the z = 3.99 galaxy.That is to use deep optical data from HSC-SSP.We fit the i-band image, in which the galaxy is detected at S/N ∼ 20, with galfit (Peng et al. 2002(Peng et al. , 2010) ) to measure its effective radius assuming the Sérsic profile.The coadd PSF from the HSC pipeline (Bosch et al. 2018) is used as an input PSF image (FWHM= 0.72 arcsec).We run a Monte-Carlo simulation by perturbing the initial guess on the centroid, position angle, effective radius, and brightness.In each run, the Sérsic index is fixed to a randomly drawn value between 0.5 and 4. The fit with the smallest χ 2 is taken as the best estimate and ∆χ 2 < 1 as the uncertainty.We obtain r ef f = 1.04 ± 0.43 kpc.Note that we use the effective radius along the semi-major axis throughout the paper.We use stellar mass from the joint photometric and spectroscopic fitting from Section 4. Fig. 4 shows the size vs. stellar mass relation including results from the literature.The sizes of our z ∼ 4 galaxies (z = 3.99 and z = 4.01) are clearly smaller than quiescent galaxies at z 2, suggesting a significant size evolution.Note that the size of the z = 4.01 galaxy is from Ito et al. (2023b).Our measurements are in line with the literature results at z ∼ 3.5; massive quiescent galaxies at these redshifts are smaller than ∼ 1 kpc, despite their large mass.This makes them a very dense system; as discussed in Kubo et al. (2018), their surface stellar mass densities are comparable to those of globular clusters.One possible formation mechanisms for such a dense system is nuclear starburst.Their star formation histories (Section 3.3) indeed suggest that they experienced an intense burst of star formation in a recent past.

Dynamical Properties
We move on to discuss dynamical properties of the z = 4 galaxies with a focus on the stellar velocity dispersion here.We do not go very far in this dynamical analysis because the velocity dispersion is not easy to estimate robustly in young galaxies like ours and the galaxies may be rotating (Toft et al. 2017;Newman et al. 2018;D'Eugenio et al. 2023), which makes the interpretation difficult.
The observed MOSFIRE spectrum of the z = 3.99 galaxy is fit with ppxf (Cappellari & Emsellem 2004;Cappellari 2017).Simple stellar population (SSP) models from Vazdekis et al. (2010) are used in the fit.Following Tanaka et al. (2019), we exclude models older than the age of the universe at the redshift of the object as well as those with subsolar metallicities as we deal with a massive galaxy.We use an additive correction function of order 1 and no multiplicative correction.Our results are not sensitive to the choices here.We run ppxf by perturbing the each pixel of the spectrum by its uncertainty and repeat the fit.We measure σ = 305 ± 103 km s −1 .
The correlation between the stellar velocity dispersion and stellar mass is plotted in Fig. 5. Again, we include the z = 4.01 galaxy (M * ∼ 2 × 10 11 M ⊙ ).The z < 1 relations are fully consistent with each other, while z ∼ 2 galaxies seem to show a tail towards larger σ.Galaxies at z ∼ 3.5 − 4 also exhibit slightly larger σ on average compared to the z < 1 relations.The scatter is large, but this may indicate a mild evolution of the stellar velocity dispersion at high redshifts.However, as men- tioned earlier, the interpretation is tricky as they may be rotating disk galaxies (Toft et al. 2017;Newman et al. 2018).It is now feasible to measure Sersic indices of high redshift galaxies with JWST (e.g., Ito et al. 2023b), and a joint analysis of morphology and stellar dynamics will be an interesting avenue.

Star Formation Histories
The spectral fitting in the previous subsection fits the observed spectrum with combination of SSP models with different ages.We thus have relative weights between different ages for each fit.We use this information to infer the star formation histories of the z = 4 galaxies.The age grid is discrete, and we apply a Gaussian kernel to make it contiguous.The inferred star formation histories (SFHs) are shown in Fig. 6.We note that the SFHs here are based on the spectrum only.The photometric data may add further constraints, and we leave in-depth SFH analyses using both the spectra and photometry for Stellar velocity dispersion is plotted against stellar mass.The blue dashed curve is for z = 0 galaxies and is from Zahid et al. (2016).The meanings of the other symbols are the same as in Fig. 4.
all the quiescent galaxies (see the next Section for the fainter quiescent galaxies) for future work.
There are both similarities and differences between the two SFHs of the z = 3.99 and 4.01 galaxies; these two galaxies experienced a starburst, followed by a rapid decrease in their SF activities.This SFH is very similar to those observed in other massive quiescent galaxies at high redshifts (Schreiber et al. 2018a;Valentino et al. 2020) and even those inferred from the local universe (e.g., Thomas et al. 2010).It is interesting that the z = 3.99 galaxy experienced the burst earlier than the z = 4.01.In other words, the z = 3.99 galaxy is older.The z = 3.99 galaxy is in a dense region as we discuss below, while the z = 4.01 galaxy is isolated.We cannot draw any conclusion with just two galaxies, but we might be seeing the environmental dependence of galaxy formation in the sense that the galaxy formation is accelerated in high density region.

JOINT PHOTOMETRIC AND SPECTROSCOPIC FITS
We turn our attention to the faint quiescent galaxies in the over-density region.While direct redshift measurements from the spectra were not possible for them due to their low-S/N (the median S/N per pixel is 0.8), the spectra trace the continuum shape around the 4000 Å/Balmer break well.We here make an attempt to constrain their redshifts by using both the spectra and photometry.As our spectra are of low-S/N and there are residuals of the night sky lines, we first bin the spectra by 100 Å by clipping the top and bottom 10% of the flux distribution in each bin.The clipping is intended to exclude the residual of the sky lines, which is non-Gaussian and often contaminates the mean flux.We  2019).The bottom axis shows the look-back time from each galaxy.The top ticks show the corresponding redshift from z = 4 (i.e., mean observed epoch).We here use the mean redshift of the two galaxies just for clarify, and the small redshift difference of ∆z = ±0.01 from the spectroscopic redshifts does not strongly affect our discussions here.The vertical axis is in arbitrary unit and each SFH is normalized to unity at the peak.The shades show the statistical uncertainties.
then feed the binned spectra and photometry to a photometric redshift code (Tanaka 2015) to constrain their redshifts.
As a proof of concept, the technique is first applied to the confirmed z = 3.99 and z = 4.01 galaxies.This "spectro-photometric" redshift technique successfully reproduces the original spectroscopic redshifts of these two galaxies as shown in Fig. 7.The redshift probability distribution functions are sharply peaked at the spectroscopic redshifts with very small uncertainties.These strong redshift constraints come from the continuum shape traced by the MOSFIRE spectra; the spectra cover the break feature and the location of the break gives a strong redshift constraint.We perform the same exercise for the spectroscopic objects presented in Ito et al. (2023a) and confirm that we can successfully reproduce the spectroscopic redshifts.
Motivated by the excellent redshift recovery, we apply the technique to the faint quiescent galaxy candidates for which no secure redshifts could be measured from MOSFIRE.Fig. 8 shows the result; all of them show a sharp redshift spike at z = 4.There is one object with a significantly bimodal P (z) distribution (SXDS2-16609), but the 2nd peak is consistent with z = 4. Due to the location of the object in the MOSFIRE mask, we could not cover the blue part of the K-band for this object.That likely resulted in the wider P (z).In any case, it is very likely that most of these faint objects are located at z = 4.The inferred redshifts of the faint quiescent galaxies re summarized in Table 1.
The spectro-photometric fits also give fairly robust estimates of physical properties of galaxies such as stellar mass and SFR (Table 2).Fig. 9 plots SFR against stellar mass.The confirmed quiescent galaxies at z = 4 are all below the star formation main sequence, and they tend to be massive galaxies.There are less quiescent galaxies and more star-forming galaxies at lower mass.We are complete to M * ∼ 3 × 10 10 M ⊙ and the decreasing quiescent galaxies towards lower mass is likely a real trend.We will elaborate on this point in Section 5.2.

LARGE-SCALE STRUCTURE AT Z = 4
5.1.A cluster of quiescent galaxies Fig. 10 shows a zoom-in view of the over-density region.The quiescent galaxies are located in a small region on the sky (∼ 2.5 arcmin radius, which is ∼ 1 physical Mpc).This is the densest region of quiescent galaxies discovered to date at such a high redshift; the projected density of quiescent galaxies is ∼ 1.5 Mpc −2 .We refer to it as cluster for now for simplicity.Interestingly, the distribution of spectroscopically confirmed star-forming galaxies from VANDELS DR4 (Garilli et al. 2021) with secure redshift flags reveals a large-scale structure extending in the N-S direction.The z = 4 cluster is clearly part of the structure, which adds further confidence.The spectroscopically confirmed z = 4.01 galaxy is located in the northern extension of the structure, and the structure probably extends outside of the VANDELS coverage.As mentioned earlier, the close redshifts of the spectroscopic objects (z = 3.99 and z = 4.01) are not coincidence; they likely belong to the same structure.
A large number of (proto-)clusters at z 3 have been reported in the literature, but all of them used star-forming galaxies as a tracer (Miller et al. 2018;Oteo et al. 2018;Toshikawa et al. 2018).The current frontier of bona-fide clusters with red sequence is z ∼ 2 (Gobat et al. 2011;Spitler et al. 2012), and recent work has just started to imply possible higher redshift clusters (McConachie et al. 2022;Ito et al. 2023a).Our cluster is the most distant cluster of quiescent galaxies discovered to date and is an ideal laboratory to address the role of environment at the highest redshift ever probed.
A ubiquitous feature of low-redshift clusters is the red sequence on a color-magnitude diagram (Bower et al. 1992) formed by quiescent galaxies.We present the color-magnitude diagram in Fig. 11.Quiescent galaxies tend to be red as expected and there is a relatively large color scatter among quiescent galaxies.However, the quiescent galaxies in the cluster shown by the filled red circles seem to exhibit a smaller color scatter and form an early red sequence.Most of the quiescent members do not exhibit large dust extinction (Table 2) and it is likely that their red colors are due to their ages.As we are very close to the primary formation epochs of these galaxies, the observed color is sensitive to a small difference in their formation redshifts.The small color scatter implies that the cluster galaxies form (and quench) about the same time.They distribute over a 1 Mpc scale and they may not form a fully virialized system yet, but our observation here seems to suggest that the red sequence first appears during the initial collapse of a massive halo, accompanied by the roughly simultaneous quenching event in each halo.This picture is different from the pre-processing scenario suggested at much lower redshifts, in which galaxies become quiescent in low-mass groups before they fall into more massive halos (e.g., Kodama et al. 2001;Fujita 2004).Environmental effects are likely stronger at lower redshifts with more abundant intracluster medium, but any effects due to the intracluster medium are not important in our case because our cluster is likely in a pre-collapse phase (see below).Possible quenching processes may include gas-exhaustion due to starburst triggered by galaxy-galaxy interactions.Interactions may enhance AGN activities, which may then trigger AGN feedback to further suppress star formation.Recent work by Kakimoto et al. (2023) reported on a massive galaxy being quenched in a very compact group at z = 4.5 and suggested a possible role of mergers for quenching at high redshifts.Their galaxy indeed has a very close companion, and it is interesting to note that this is not the only case where a massive quiescent galaxy accompanies a close companion (e.g., Schreiber et al. 2018b).A close companion may indicate a potential role of very small-scale environment for quenching.
Star-forming galaxies, on the other hand, are separated from the red sequence and form a cloud of blue galaxies, which indicates that the bimodality of galaxy populations is already in place as early as z = 4. Overall, the segregation of galaxy properties that we observe in the local Universe has already emerged at this epoch.

The Butcher-Oemler effect at z = 4
We further discuss the environmental dependence of galaxy SFRs by extending the Butcher-Oemler effect (Butcher & Oemler 1978) to z = 4 for the first time.We first define the cluster region using a 1 Mpc aperture that encloses the concentration of the quiescent galaxies and estimate the quiescent fraction.Here, we include photometric redshift selected star-forming galaxies within the aperture as blue members.As in Fig. 2, we use galaxies  The error bars at the bottom-left corner indicate the typical fractional uncertainties including uncertainties in z phot .The filled red circles are the quiescent galaxies at z = 4 with secure spectroscopic/spectro-photometric redshifts.Recall that the quiescent galaxies are defined as those with a 1σ upper limit of sSFR being below 10 −9.5 yr −1 as indicated by the dashed line.Galaxies below SFR< 0.1M⊙ yr −1 are all shown at SFR= 0.1M⊙ yr −1 .at 3.7 < z phot < 4.3.We apply the same selection to all galaxies outside the aperture to construct a field sample for comparison.The fractions of quiescent galaxies are shown in Fig. 12.
We begin with the global trend using the entire (=cluster + field) sample.The black points nicely show that the quiescent fraction is a strong function of stellar mass in the sense that it increases towards higher mass.The fraction is as high as ∼ 30% at M 10 11 M ⊙ .On the other hand, it is very low at low mass and low-mass galaxies are predominantly star-forming galaxies.This trend is consistent with Weaver et al. (2022).In other words, quiescent galaxies are almost exclusively massive galaxies.
Turning to the cluster region, we do not have sufficient statistics to discuss the quiescent fraction as a function of stellar mass, and we instead estimate the quiescent fraction using all the cluster members.We take the mean stellar mass here and show the quiescent fraction in Fig. 12 as the filled red circle.The quiescent fraction is about 35%.For comparison, we do the same for the field galaxies and show their quiescent fraction as the blue circle.The field quiescent fraction is only 1.5% and is significantly lower than that of the cluster.While we have only a single cluster, our finding here indicates that a cluster environment exhibits a higher quiescent fraction from an early phase of its formation.In other words, the Butcher-Oemler effect might be already there at z = 4.
A larger cluster sample is clearly needed to confirm the trend, but it will be a significant challenge to identify a cluster like ours as one needs to survey a large volume and then push the current spectroscopic facilities to the limit to obtain secure redshifts.Nonetheless, it is worth the effort to search for more clusters at higher redshifts to peer deeper into the cluster formation.

SUMMARY AND DISCUSSION
We have presented the spectroscopic confirmation of a massive quiescent galaxy at z = 3.99 as well as the discovery of a significant concentration of quiescent galaxies at z = 4. Based on the deep spectroscopic observation with MOSFIRE, we have identified multiple Balmer absorption lines of the z = 3.99 quiescent galaxy.Its stellar velocity dispersion is σ = 305 ± 103 km s −1 and is slightly higher than lower redshift counterparts with similar stellar masses, indicating a mild redshift evolution.Its physical size measured from the deep imaging data reveals its compact size and dense stellar density, consistent with other quiescent galaxies observed at similar redshifts in the literature.
Interestingly, the galaxy is located in an over-density region of quiescent galaxies at similar redshifts.While we could not directly measure their spectroscopic redshifts from their spectra, the high-accuracy spectrophotometric redshifts strongly suggests that they are at the same redshifts as the z = 3.99 galaxy.All these galaxies exhibit a pronounced Balmer break, indicative of recent starburst activity followed by rapid quenching.Their UV emission is relatively weak and their SFRs inferred from their overall SEDs indeed suggest that they are quiescent.The concentration is embedded in a large-scale filament traced by the VANDELS spectroscopic galaxies, which adds further credits to the physical association.This is the first discovery of a concentration of quiescent galaxies at this high redshift and it serves as an unique laboratory to investigate the role of environment at z = 4 for the first time.The cluster galaxies form the early red sequence, suggesting that these galaxies form and quench at similar times, and the red sequence emerges during the first collapse of a cluster.In addition, the quiescent galaxy fraction is significantly higher than the field average.This may be expected because we define the over-density region around the quiescent galaxies, but we have extended the Butcher-Oemler effect to z=4 for the first time.The fact that the collapsing cluster exhibits a higher quiescent fraction suggests that we may be witnessing the environmental dependence of galaxy formation such that galaxies in higher density regions form earlier than the field galaxies.We aim to address this point further in our future paper.
Given the significant over-density, we expect that the cluster region will collapse to form the core of a massive cluster at lower redshifts.We make an attempt to roughly estimate the future halo mass of the system 11.5 +7.9 −6.0 0.9 +0.1 when they collapse.Using the halo occupation distribution modeling from Shuntov et al. (2022), we infer the halo mass of each quiescent galaxy from its stellar mass as summarized in Table 2.The individual halo mass ranges from 1 to 3 × 10 12 M ⊙ .If we assume that they all collapse to form a single halo, the total mass exceeds 10 13 M ⊙ , i.e., group-sized halo.We consider only the quiescent members, but there are star-forming galaxies with consistent photo-z's, and thus the total halo mass here is a lower limit.We find that the virial radii of the galaxies derived from their halo mass do not overlap with each other, suggesting that they are not yet in the same halo at the time of the observation.
In order to further characterize the properties of the cluster, we have searched for similar systems in the Illustris TNG-300 simulation at z = 4.01 (snapshot 21;

Nelson et al. 2019).
To be specific, we identify quiescent galaxies with M * > 10 10 M ⊙ and log sSFR < 10 −9.5 M ⊙ yr −1 in the simulation and examine their spatial distribution.Stellar mass and SFR are measured within twice the stellar halo-mass radius, and SFR is averaged over 10 Myr.Interestingly, there is no over-density region of the quiescent galaxies in TNG-300; there are only 11 massive quiescent galaxies within the 300 comoving Mpc simulation box and the closest separation between them is as large as 8 physical Mpc (3d distance).In contrast, we have 5 quiescent galaxies within 1 physical Mpc radius (projected distance).While the 3d and projected distances are different, it is clear that the observed cluster is a much denser system.We note that we perform the same analysis but considering more massive galaxies with > 5 × 10 11 M ⊙ only and confirmed there is still no counterpart in Illustris TNG-300.
We do not know yet if this is a failure of the simulation in reproducing quiescent galaxies at z ∼ 4; the number density of massive quiescent galaxies at 3 z 4 is in broad agreement between observation and simulations (Valentino et al. 2023 and references therein), although it depends strongly on how quiescent galaxies are defined.It may be that our system is so rare that the simulation box is too small to reproduce it.In any case, massive halos at high redshifts have to be rare and it is crucial to increase the cluster sample to perform further comparisons with simulations.
We expect that our spectra of the faint cluster members will be able to strongly constrain their star formation histories as the spectra cover the Balmer/4000 Å break for most of them.All of the members show a sign of suppressed star formation as discussed earlier and it would be of interest to perform extensive SED fitting to infer their primary formation epoch and quenching timescale and compare them with those of the field galaxies from the literature.Also, we can measure the physical sizes of the members and see whether there is environmental dependence of galaxy sizes at these redshifts.The environmental dependence of physical sizes of galaxies at low redshifts is somewhat controversial in the literature (see Table 1 of Yoon et al. 2017, andTable A1 of Matharu et al. 2019), and it would be interesting to look at the z = 4 system.We hope to address these points in a forthcoming paper.Also, detailed imaging/spectroscopic follow-up observations with JWST are an interesting avenue as they deliver secure spectroscopic redshifts, more accurate physical sizes, etc. ALMA observations to search for dustobscured galaxies that went unnoticed in our work will be interesting, too.This work can be extended in multiple ways, and joint analyses will hopefully lead us to deeper insights into the early cluster formation and the quenching physics in the early Universe.This work is supported by JSPS KAKENHI Grant Numbers JP23740144, 22KJ0730, JP15K17617.The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation (DNRF140).
Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration.The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community.We are most fortunate to have the opportunity to Facilities: Keck(MOSFIRE) APPENDIX A. DATA REDUCTION In this Appendix, we detail our data reduction procedure.Each A-B pair was flat-fielded and wavelength-calibrated by using common flat and wavelength reference frames.The background was removed by subtracting the frame at the B dither position from that at the A dither position followed by subtracting residuals (Kelson 2003).Finally, all A-B pairs were rectified to the same wavelength solution, while coaddition was done for each A-B pair separately by applying the dither offset at this stage.The slit mask alignment procedure was carried out every ∼ 2 hours, but there were still noticeable slit drifts in the spatial direction between the alignment procedures (Kriek et al. 2015;Larson et al. 2022).We placed a bright point source (SXDS2 19063, a quasar) to monitor the drift and the seeing variation.The drifts up to ∼ 2 pixels were corrected to locate objects in a common spatial frame.
We observed an A0V star (HIP 111996) for telluric correction twice per night with the long2pos slit mask which is designed to take an AB dither at 2 positions (Pos A and Pos C), which enables us to cover the entire wavelength range of K-band.We reduced spectra with the DRP at each position separately and extracted one-dimensional spectra by applying the optimal extraction method (Horne 1986) assuming a Gaussian spatial profile.The extracted spectra were normalized at 2.1µm and then combined.We used Pos A spectra at λ > 2.3µm, while spectra at Pos C were used for the rest of the wavelength range.This is because we found that the S/N is higher at Pos C most likely due to a misalignment of the slit geometry and preset of the telescope offsets.The telluric correction vectors were created by dividing the telluric standard spectra by a Vega template spectrum (Kurucz 1993) after convolving to the MOSFIRE resolution of R=3600.The telluric correction vector was applied for each two-dimensional spectrum closest in time to the observation time of HIP 111996.
The coaddition of two-dimensional A-B stacks was carried out by taking a weighted mean using the peak counts measured for the point source used for the drift correction.Two among 159 A-B stacks were not used since they showed too low S/N to robustly derive the traces of the source due to cloudy conditions.We also applied 3-sigma clipping to remove bad pixels and cosmic rays.Error spectra were derived by propagating error spectra associated with each A-B stack.The total integration time of the coadded spectra used in the following analysis is 15.7 hours.
Absolute flux calibration was made by using the point source common for all science exposures.We used an optimally extracted one-dimensional spectrum for the object and derived a synthetic Ks-band magnitude by convolving it with the VISTA Ks-band filter response curve for each A-B stack.We also corrected the slit loss by comparing the slit width and a Gaussian FWHM measured along the spatial direction.Then, we applied normalization factors by comparing the synthetic magnitudes with the observed VISTA Ks-band total magnitude (Mehta et al. 2018) to scale the spectra to the absolute fluxes.The spectra fully calibrated in this way are used in the main body of the paper.
Figure 1.Photometric redshifts plotted against the secure spectroscopic redshifts from VANDELS DR4.The error bars indicate the 68% confidence intervals.The red circles are spectroscopically confirmed quiescent galaxies fromSchreiber et al. (2018a),Tanaka et al. (2019), and this paper.We differentiate spectrophotometric galaxies (see Section 4) from the secure spectroscopic redshifts with the open triangles.The dashed line shows the z phot = zspec relation.

Figure 2 .
Figure2.Distribution of massive (M * > 10 10 M⊙) galaxies at z ∼ 4 in SXDS.The open red/blue circles are photo-z selected quiescent/star-forming galaxies at 3.7 < z phot < 4.3.The color contours show the density of quiescent galaxies smoothed over 1 Mpc (physical).The filled red points are galaxies in the over-density region and are the focus of the paper.There is a spectroscopically confirmed z = 4.01 galaxy in the north of the over-density region(Tanaka et al. 2019).The green square is the field field of view of our MOSFIRE observation.The bar on the bottom-left indicates 1 Mpc (physical) at z = 4.

Figure 4 .
Figure 4.Effective radius plotted against stellar mass.The dashed line is the z = 0 relation fromShen et al. (2003), and the blue shades are quiescent galaxies at z ∼ 0.7 from the LEGA-C survey(van der Wel et al. 2016;Straatman et al. 2018).Quiescent galaxies at z ∼ 2 are collected from the literature(van Dokkum et al. 2009;Onodera et al. 2012;Toft et al. 2012;Bezanson et al. 2013; van de Sande et al.  2013;Belli et al. 2017;Stockmann et al. 2020) and are indicated as the black points.Another collection from the literature of z ∼ 3.5 galaxies are shown in pink(Saracco et al. 2020;Esdaile et al. 2021;Forrest et al. 2022).The stacked object ofKubo et al. (2018) is shown as the filled triangle.The gray symbols are the JWST measurements of quiescent galaxies at z > 3 fromIto et al. (2023b).Note that they are photo-z selected galaxies (i.e., no spectroscopic confirmation yet).Our z = 4 objects are in red.The error bars show the statistical uncertainties.The solid and dotted lines in the bottom right corner show evolutionary tracks with r ef f ∝ M * and r ef f ∝ M 2 * , respectively.They represent the major and minor merger tracks(Bezanson et al. 2009;Naab et al. 2009).
Figure 5.Stellar velocity dispersion is plotted against stellar mass.The blue dashed curve is for z = 0 galaxies and is fromZahid et al. (2016).The meanings of the other symbols are the same as in Fig.4.

Figure 6 .
Figure6.Star formation histories of the z = 3.99 (red) and z = 4.01 (blue) galaxies.The latter object is fromTanaka et al. (2019).The bottom axis shows the look-back time from each galaxy.The top ticks show the corresponding redshift from z = 4 (i.e., mean observed epoch).We here use the mean redshift of the two galaxies just for clarify, and the small redshift difference of ∆z = ±0.01 from the spectroscopic redshifts does not strongly affect our discussions here.The vertical axis is in arbitrary unit and each SFH is normalized to unity at the peak.The shades show the statistical uncertainties.

Figure 7 .
Figure 7. Spectro-photometric fitting of the z = 3.99 (left) and z = 4.01 (right) galaxies.The top panel shows the observed SED.The large points are the broad-band photometry and the small points around λ ∼ 2µm are the binned MOSFIRE spectrum.The black spectrum is the best-fitting model template and the blue squares are the model fluxes.The bottom panel shows the redshift probability distribution function.It is sharply peaked at the spectroscopic redshift for both objects.

Figure 8 .
Figure8.Spectro-photopmetric fits to the faint quiescent galaxy candidates.The meanings of the symbols are the same as in Fig.7

Figure 9 .
Figure9.SFR plotted against stellar mass for galaxies at z ∼ 4. The blue dots and red open circles are star-forming and quiescent galaxies at 3.7 < z phot < 4.3, respectively.The error bars at the bottom-left corner indicate the typical fractional uncertainties including uncertainties in z phot .The filled red circles are the quiescent galaxies at z = 4 with secure spectroscopic/spectro-photometric redshifts.Recall that the quiescent galaxies are defined as those with a 1σ upper limit of sSFR being below 10 −9.5 yr −1 as indicated by the dashed line.Galaxies below SFR< 0.1M⊙ yr −1 are all shown at SFR= 0.1M⊙ yr −1 .

Figure 10 .
Figure10.Blow-up of the cluster region.The quiescent galaxies at z ∼ 4 are shown as the filled red points.Those with secure spectroscopic redshifts are indicated with the redshift labels.The green stars are spectroscopically confirmed 3.95 < zspec < 4.05 galaxies from VANDELS, and gray box indicates the region surveyed by VANDELS.The meaning of the other symbols are the same as in Fig.2.The small panel on the top shows the redshift distribution of the VANDELS galaxies.The arrow indicates z = 3.99.There is a redshift spike at z ∼ 3.98.The distribution of the green stars in the main panel remains virtually the same if we adopt a narrower redshift slice of, e.g., 3.95 < z < 4.01.

Figure 11 .
Figure11.Observed H − K color plotted against K-band magnitude.The meanings of the symbols are the same as in Fig.9, except that the filled red circles show the quiescent galaxies in the over-density region only (i.e., z = 4.01 galaxy is excluded).Objects with secure redshifts are indicated.

Figure 12 .
Figure12.Fraction of quiescent galaxies as a function of stellar mass based on galaxies 3.7 < z phot < 4.3.The black open circles are for the combined field and cluster sample to see the global trend and are plotted as a function of stellar mass.The red and blue filled circles are the average quiescent fractions of the cluster and field regions, respectively.These circles are plotted at the mean stellar mass.The vertical error bars are the Poisson errors in the quiescent fraction, while the horizontal error bards show the 68% interval of the stellar mass distribution.