JWST’s First Glimpse of a z > 2 Forming Cluster Reveals a Top-heavy Stellar Mass Function

Clusters and their progenitors (protoclusters) at z ∼ 2 − 4, the peak epoch of star formation, are ideal laboratories to study the formation process of both the clusters themselves and their member galaxies. However, a complete census of their member galaxies has been challenging due to observational difficulties. Here we present new JWST/NIRCam observations targeting the distant cluster CLJ1001 at z = 2.51 from the COSMOS-Web program, which, in combination with previous narrowband imaging targeting Hα emitters and deep millimeter surveys of CO emitters, provide a complete view of massive galaxy assembly in CLJ1001. In particular, JWST reveals a population of massive, extremely red cluster members in the long-wavelength bands that were invisible in previous Hubble Space Telescope (HST)/F160W imaging (HST-dark members). Based on this highly complete spectroscopic sample of member galaxies, we show that the spatial distribution of galaxies in CLJ1001 exhibits a strong central concentration, with the central galaxy density already resembling that of low-z clusters. Moreover, we reveal a “top-heavy” stellar mass function for the star-forming galaxies (SFGs), with an overabundance of massive SFGs piled up in the cluster core. These features strongly suggest that CLJ1001 is caught in a rapid transition, with many of its massive SFGs likely soon becoming quiescent. In the context of cluster formation, these findings suggest that the earliest clusters form from the inside out and top to bottom, with the massive galaxies in the core assembling first, followed by the less massive ones in the outskirts.

One of the most prominent features of local galaxy clusters is a high concentration of massive quiescent galaxies (QGs) in their cores (Dressler 1980).The formation process of these massive galaxies and, in particular, the role of dense environments in shaping their formation and evolution remains unclear.Galactic archeology suggests that most of their stars were arXiv:2403.05248v3[astro-ph.GA] 30 May 2024 formed at z ∼ 2 − 4 (e.g.Thomas et al. 2010), which makes (proto)clusters at these redshifts ideal laboratories to investigate their formation process.During the last decade, a number of spectroscopically confirmed (proto)clusters have been found at these redshifts (e.g.Pentericci et al. 2000;Kajisawa et al. 2006;Gobat et al. 2011;Hayashi et al. 2012;Overzier 2016;Wang et al. 2016;Noirot et al. 2018;Oteo et al. 2018;Zhou et al. 2020Zhou et al. , 2024)).Studies of these early formed structures and their member galaxies have provided important insights into the formation process of both clusters and cluster galaxies (e.g.Strazzullo et al. 2016;Afanasiev et al. 2023;Mei et al. 2023;Pérez-Martínez et al. 2023).However, it is still an open question how these early formed (proto)clusters evolve into low-redshift, more mature clusters and how the dense environments impact the formation and quenching of their member galaxies.
The major difficulty in studying galaxy formation in these high-z (proto)clusters is obtaining an unbiased census of their member galaxies.In most cases, only a small fraction of member galaxies at any given mass are spectroscopically confirmed.Without a representative and unbiased sample of cluster members, it is challenging to gauge the actual relevance of any environmental dependence.So far, various efforts have been made to obtain an unbiased census of member galaxies in a few structures at z ≳ 2, which are based on deep near-infrared (NIR) spectroscopy (e.g.Gobat et al. 2013;Balogh et al. 2021), deep CO spectroscopy (Wang et al. 2018;Jin et al. 2021), or narrowband imaging (Hayashi et al. 2016;Shimakawa et al. 2018a,b;Zheng et al. 2021).However, the prevalence of dusty starbursts and crowded distribution of galaxies in these structures often require combining several methods in order to fully and unbiasedly cover the member galaxy population.
In this Letter, we combine the deep-and highresolution images from JWST/NIRCam, narrowband imaging from Subaru/MOIRCS, and CO spectroscopy from the Atacama Large Millimeter/submillimeter Array (ALMA)/NOEMA/Very Large Array (VLA) to provide a complete census of member galaxies in the z spec = 2.51 cluster CLJ1001 (hereafter, J1001; Wang et al. 2016).Previous studies of J1001 have found that it has enhanced star formation rate (SFR), short gas depletion time and strong evidence that the properties of member galaxies are affected by the extreme environment (Wang et al. 2016(Wang et al. , 2018;;Gómez-Guijarro et al. 2019;Xiao et al. 2022;Xu et al. 2023).Here we further show a population of Hubble Space Telescope (HST)-undetected red member galaxies with bluer companions, which shows that even at z ∼ 2.5, the high-resolution NIR-to-mid-infrared images from JWST are necessary to uncover the com-plete massive galaxy population of (proto)clusters.In addition, we find that J1001 has a highly concentrated density profile and a top-heavy stellar mass function (SMF), which support an inside-out and top-to-bottom evolutionary path.
The Letter is organized as follows: Section 2 presents the data set used for this work; Section 3 shows the methods used to select member galaxies and obtain their physical properties; Section 4 presents the main results, including a population of HST-dark cluster members, the density profile, and stellar mass function.In this work, we assume the cosmological model with 3 and Ω Λ = 0.7.The AB system (Oke & Gunn 1983) is used for the magnitudes throughout the Letter, and a Kroupa (2001) initial mass function is adopted for the estimation of stellar mass.

JWST/NIRcam images
Cluster J1001 has been observed by JWST/NIRCam with four filters (F115W, F150W, F277W, and F444W) as part of the COSMOS-Web survey (Casey et al. 2023) in December 2023.We reduce the raw JWST/NIRCam images from the Mikulski Archive for Space Telescopes (MAST)1 with the JWST Calibration Pipeline v1.12.5 (Bushouse et al. 2023).We run both stages 1 and 2 of the pipeline with default parameters.Then, we perform additional subtraction of the stripe-like 1/f noises (e.g.Bagley et al. 2023), background structure, and some remaining cosmic rays (e.g.Rieke et al. 2023) for each "cal" file output by stage 2.During stage 3, we use a reference catalog based on previous HST/Wide Field Camera 3 (WFC3) F160W image (Xu et al. 2023) for astrometry calibration.

Other multiwavelength data
The deep HST/WFC3 F125W and F160W images of J1001 have been obtained during Project 14750 (PI: T. Wang) with details shown in Xu et al. (2023).These mosaics from HST have been resampled to match the pixel size of JWST/NIRCam images with Montage (Berriman et al. 2003;Good & Berriman 2019).In addition, to identify Hα emitters (HAEs) within J1001, narrowband imaging with the NB23002 filter on Subaru/MOIRCS has been conducted (PI: T. Kodama).The upper panel of Figure 1 shows the transmission curve of the NB2300 filter.According to the FWHM of this narrowband filter, it can be used to identify HAEs at 2.483 < z < 2.519 combined with the deep K s imaging from UltraVista.The membership selection of J1001 with this high accuracy is essential due to the presence of a number of protoclusters at similar redshifts as J1001 in the same field (Casey et al. 2015;Cucciati et al. 2018;Huang et al. 2022).Based on this narrowband observation, the survey area of this work is taken as the maximum square area (3.48 ′ × 3.48 ′ , or 1684kpc × 1684kpc) with narrowband coverage, which is shown as the blue square in Figure 2 (left).
CO spectroscopy in the millimeter toward J1001 has been conducted through a few programs with ALMA, NOEMA, and VLA.Most of them have been reported in previous studies (Wang et al. 2016(Wang et al. , 2018;;Champagne et al. 2021;Xiao et al. 2022).New CO observations include a high-resolution CO(1-0) observation with VLA and new CO(3-2) observations with ALMA toward a larger area than that reported in Xiao et al. (2022).Details of these observations will be described in a forthcoming paper (Xu et al., in preparation).Here we only use the redshift information as derived from the CO spectroscopy, totaling 23 CO(3-2) emitters including 11 CO(1-0) emitters, to select member galaxies.
Additionally, because J1001 is located in the COS-MOS field with abundant multiwavelength data, we further retrieve archival multiband images from U to IRAC channel 4 of the cluster area and perform source photometry on them.Detailed information on these images is given in Appendix A.

Catalog construction
Based on the multiband images described in Section 2, we build a multiwavelength photometric catalog for CLJ1001.We first convolve the F814W, F115W, F150W, and F277W images to match the point spread function (PSF) of F444W for color fidelity.Next, we perform source detection by running Source-Extractor v2.25.0 (Bertin & Arnouts 1996, 2010) on the convolved F277W image, which is the deepest JWST/NIRCam image of J1001, and most sensitive to galaxies at z ∼ 2.5.Then, we measure the total flux for each F277Wdetected source in each band with Source-Extractor and T-PHOT v2.0 (Merlin et al. 2015(Merlin et al. , 2016)); details of the photometric methods are described in Appendix B.

Sample selection
Based on the photometric catalog, we have consistent color measurements for every F277W-detected source in The red stars, blue circles, and gray points are HAEs, spectroscopically confirmed members and interlopers, respectively.The magenta solid line shows the limitation given by Equation 1 with Σ = 2.The magenta dashed line corresponds to the color cut given by Equation 2.
the broad K s and narrow NB2300 band, which are used to search for HAEs in J1001.The K s and NB2300 fluxes are both measured by T-PHOT with the same configuration of the code.Following Shimakawa et al. (2018a), the selection criteria of HAEs are where f NB is the flux density in the NB2300 narrowband, σ Ks and σ NB are the 1σ limiting flux density at the broadband K s and narrowband NB2300, respectively; Σ is the confidence level in σ, the color cut 0.35 (Shi-makawa et al. 2018a) corresponds to a rest-frame equivalent width (EW) limit of 45 Å for the Hα line; and the color term 0.13 is applied since narrowband NB2300 is at the red end of the broadband K s , which may lead to an overestimation of the Hα line flux for red sources.In addition, we also require that the signal-to-noise ratio (SNR) of f NB should be larger than 2 (SNR NB > 2) for all HAE candidates.To reject emitters of other emission lines (e.g.[OIII] emitters at z ∼ 3.6), we fit the photometric redshift (z phot ) of each source using code EAZY (Brammer et al. 2008; the minimum χ 2 redshifts are used) with agn blue sfhz 134 templates and apply 5% systematic error floor.To calculate the uncertainty on the photometric redshifts, we compare the spectroscopic and photometric redshifts using the 28 spectroscopic members.We find two outliers with |z phot − z spec | > 0.5 (both are heavily blended), and the standard deviation of the other members is 0.043.Then, we only keep those line emitters with 2.374 < z phot < 2.638, i.e., within the 3σ standard deviation of |z phot − z spec | as HAEs.These HAEs can all be considered to have 2.483 < z < 2.519 according to the FWHM range of the NB2300 filter.
Figure 1 shows the color-magnitude diagram used to select HAEs.Here, 64 HAEs are selected by the criteria with Σ > 2, 24 of which are spectroscopically confirmed members with either Hα detections by VLT/KMOS (Wang et al. 2016) or the CO spectroscopy mentioned in section 2.2.In addition to the 64 HAEs, we further include four non-HAE spectroscopic members.All of these 68 narrowband or spectroscopic members can be selected as SFGs with their rest-frame UVJ color (Williams et al. 2009;Carnall et al. 2018).The completeness of this sample of SFGs reaches ∼ 80%(50%) at log(M * /M ⊙ ) = 9.5(9.1)(see Appendix D).Lastly, we also include six UVJ-selected quiescent members with 2.374 < z eazy < 2.638 (three out of six quiescent members are close neighbors of confirmed members).This yields a total number of 74 cluster members.
To derive stellar masses of these 74 member galaxies, we use the code BAGPIPES (Carnall et al. 2018) within a narrow redshift range [2.481, 2.531] centered at z = 2.506, with stellar population synthesis model in 2016 version of Bruzual & Charlot (2003) (age∈ [0.1, 10] Gyr and metallicity log Z/Z ⊙ ∈ [0, 2.5]), a delayed star formation history (timescale τ ∈ [0.3, 10] Gyr), the dust attenuation law from Calzetti et al. (2000) for young and old populations separately (divided by 0.01 Gyr, A V ∈ [0, 5]), and nebular emission (Byler et al. 2017) (ionization parameter log U ∈ [−5, −2]).For members of J1001, the typical uncertainty of stellar mass is 0.06 dex, which includes the uncertainty of all other parameters and the photometric uncertainty.This has been significantly suppressed by the deep JWST/NIRCam F150W, F277W, and F444W images since they directly trace the stellar masses at z ∼ 2.5.Two examples of the best-fit SED are shown in Appendix C.

A population of red, massive, and HST-dark cluster members
The spatial distribution of our selected members is shown in Figure 2 (left), with all member galaxies marked as open circles.Among these member galaxies, we find four massive red members, ∼ 16% of the population at M ⋆ > 10 10.15 M ⊙ , which have been missed by previous blind source detections on the HST/WFC3 F160W image.The right panels of Figure 2 show the cutouts of these four HST-dark members.All of them have been spectroscopically confirmed via CO(3-2) emission by ALMA (Xu et al. in preparation).Considering the fluxes within a small aperture (d = 0.32"), these HST-dark members have mag F150W,aper ∼ 27 and mag F150W,aper −mag F444W,aper ∼ 3.However, their total fluxes at H band can be much higher (mag F150W,total ∼ 24), which means that these HST-dark members are extremely extended at short wavelengths.In this case, since they are often associated with one or more brighter (and bluer) companions, blind source detection on the F160W image cannot identify them and might recognize them as extended structures (e.g.tails) of their brighter companions.This suggests that the merger rate of galaxies and stellar mass function in J1001-like forming (proto)clusters can be underestimated without the JWST observation.

A highly concentrated galaxy density profile
We quantify the spatial distribution of member galaxies in J1001 via their projected density profile (Figure 3).To make a fair comparison with the results at z ∼ 1 (van der Burg et al. 2014) and z ∼ 0.15 (van der Burg et al. 2015), we also only use galaxies with M * > 10 10 M ⊙ to derive the density profile and adopt a similar IMF. Figure 3 shows that the central density of J1001 can be even higher than clusters at lower redshifts, while its density in the outer region is much lower than mature clusters.According to the result from the Millennium Simulations (Chiang et al. 2013), a cluster with M 200,J1001 = 10 13.9±0.2M ⊙ at z ∼ 2.5 will have M 200 = 6.0 +3.5  −2.2 × 10 14 M ⊙ at z ∼ 1.In this case, J1001  can be considered to be the predecessor of those z ∼ 1 clusters (van der Burg et al. 2014) within the uncertainties.This, combined with the results in Figure 3 supports an inside-out formation scenario (van der Burg et al. 2015) up to z ∼ 2.5, which means the massive central galaxies in clusters are formed ahead of the smaller galaxies in the outskirts.
We then fit the data points in Figure 3 with the projected Navarro-Frenk-White (NFW) profile (Navarro et al. 1995;Bartelmann 1996)  Figure 3 (right) compares the c number of J1001 with that of other clusters at lower redshifts (Budzynski et al. 2012;Strazzullo et al. 2013;Annunziatella et al. 2014;van der Burg et al. 2014van der Burg et al. , 2015)), as well as massive dark matter halos (M halo = 10 14 M ⊙ ) from numerical simulations (Duffy et al. 2008).Previous studies of low-redshift clusters show that the concentration of member galaxy distribution decreases with time below z ∼ 1.5 (van der  (2014,2015).Following van der Burg et al. ( 2015), the original results with dimensionless units are converted to have physical units by assuming M200 = 3 × 10 14 M⊙ at z ∼ 1 and M200 = 9 × 10 14 M⊙ at z ∼ 0.15.Field level is given by galaxies with 2.483 < z phot < 2.519 from the COSMOS2020 catalog (Weaver et al. 2022).The vertical dotted line shows the R200 of J1001.Middle panel: similar to the left panel, but the number density profiles are shown.Right panel: The redshift evolution of the concentration of cluster-galaxy number density profiles from observations (filled circles) and of the density profile of massive (M h = 10 14 M⊙) dark matter halos (dotted line) based on simulations from Duffy et al. (2008).Burg et al. 2015;Ahad et al. 2021), and we further push this increasing trend of concentration with redshift for cluster galaxies to z ∼ 2.5.These observed concentrations of galaxies can be much higher than dark matter halos because the massive galaxies (or subhalos) tend to be concentrated in the center due to dynamical friction, which makes the distribution of galaxies deviate from the host halo (e.g.Han et al. 2016).We note that the surveys shown here have different fitting ranges and virial radii, which may also affect the concentration, and we do not compare c mass since it is too sensitive to the presence of massive central galaxies.

4.3.
A "Top-heavy" stellar mass function of star-forming member galaxies The combination of deep JWST/NIRCam imaging, Subaru/MOIRCS narrowband imaging, and ALMA/NOEMA/VLA CO line surveys yields a complete census of star-forming (SF) members at M * ≳ 10 9.5 M ⊙ (see Appendix D) in J1001.Based on this sample, we determine the SMF of SF members in J1001 (Figure 4).Instead of using the data points in Figure 4, we perform the maximum-likelihood fitting to avoid the arbitrary binning procedure.We fit the SMF with both the single Schechter function  After Quenching J1001 (SFGs,original) J1001 (SFGs,after quenching) J1001 (QGs,after quenching) z 1.2 clusters (vdB+20) (SFGs) z 1.2 clusters (vdB+20) (QGs) Figure 5. Left panel: similar to Figure 4, while all the massive central galaxies within 50 kpc are merged into a quiescent BCG.The green dotted line is a copy of the blue line in Figure 4.The dashed line is the SMF of clusters at z ∼ 1.2 (van der Burg et al. 2020) within 1Mpc of their center (comparable with the size of our survey).Right panel: the SF and quiescent SMF of cluster J1001 with all galaxies to be quenched within 0.5Gyr considered as QGs.The blue and red dashed lines show the SF and quiescent SMF of cluster galaxies at z ∼ 1.2, respectively.
which is a single gamma distribution of the term 10 log M −log M * , and also the double Schechter function which is the mixture of two single Schechter functions and is similar to the mixtures of gammas discussed by Young et al. (2019) Compared with the SMF shape of other protoclusters (Shimakawa et al. 2018a,b;Edward et al. 2023), the SMF of SFGs within J1001 shows a clear top-heavy feature around M * ∼ 10 10.5−11.0M ⊙ .This means that J1001 has an excess of massive SFGs, most of which are concentrated in the cluster core (Figure 2).This is consistent with the high concentration of J1001 shown in Figure 3. Low-redshift clusters only exhibit a topheavy structure when both the QGs and SF members are combined (e.g.van der Burg et al. 2013Burg et al. , 2018Burg et al. , 2020;;Annunziatella et al. 2014).A direct indication is that the massive SF members in J1001, at the time of the observation, likely include a large population of progenitors of massive QGs.Another possibility is that many of these massive SFGs may merge into a single brightest cluster galaxy (BCG), considering that many of them are located within a small volume.
We illustrate these two scenarios on the possible evolutionary path of the SMF in J1001 in Figure 5. Firstly, consider that most of the massive SFGs in J1001 in the central area might soon be merged into a BCG, which is currently still missing for J1001.Specifically, we assume all galaxies within 50 kpc of the cluster center will be merged into a quiescent BCG, which can then be excluded from the SMF of SFGs.In this case, as shown in the left panel of Figure 5, the top-heavy feature will be obviously suppressed with ∆(BIC) decreasing from 11.02 to 1.01.
Secondly, we consider the effect of quenching on the evolution of the SMF of both SFGs and QGs.We collect the gas depletion time t dep ≡ M gas /SFR for members with CO spectroscopy (see section 2.2) and assume that the SFGs with t dep < 0.5Gyr will soon become QGs.After this transformation, we derive the yielding SMF of SFGs and QGs in J1001, as shown in the right panel of Figure 5.We find that the top-heavy feature in the original SMF of SFGs in J1001 is strongly suppressed with ∆(BIC) decreasing to 0.96.Meanwhile, a rapid rise for the SMF of QGs in J1001 is observed, which surpasses that of SFGs at the massive end and matches QGs in low-redshift clusters.
The above analysis shows that both mergers and quenching can help eliminate the "top-heavy" feature in the SMF of SFGs in J1001.In practice, both processes might take place simultaneously.Moreover, by comparison with the SMF in z ∼ 1.2 clusters from the GOGREEN program (van der Burg et al. 2020), we show that the massive end of J1001 is already as abundant as the clusters at z ∼ 1.2, while the number of low-mass members with log M * = 9−10 is much smaller.We have briefly discussed that J1001 can be approximately considered as the predecessor of clusters observed at z ∼ 1 (section 4.2).In addition, the dearth of small members in J1001 is unlikely to be an observational effect since our sample should be complete at log M * > 9.5 (Appendix D), and all the expected incompleteness has already been corrected.This means that the massive cluster galaxies are already in place at z ∼ 2.5, suggesting a top-to-bottom growth of these early formed clusters.

DISCUSSION AND CONCLUSIONS
Based on the deep multiwavelength observations including JWST/NIRCam and narrowband imaging with Subaru/MOIRCS, we have performed a complete census of the member galaxies in cluster J1001 down to M ⋆ ≳ 10 9 M ⊙ at z = 2.51.In this Letter, we focus on the global properties of this sample of member galaxies, including their density profile and stellar mass function, aiming to constrain the evolutionary state of J1001 in the context of cluster formation at high redshifts.Our main findings are listed as follows: 1. JWST/NIRCam observations reveal a population of red and massive cluster members, which have been missed from previous deep HST/F160W imaging (HST-dark).In addition to their faintness, most of them have a close, bluer HST companion.This inhibits blind detections of the red sources in HST images, which only become detectable at longer wavelength by JWST/NIRCam with better sensitivity and resolution.The prevalence of the massive, HST-dark cluster members suggests that JWST is necessary to obtain a complete census of member galaxies, even at z ∼ 2.5.Consequently, previous estimates based on HST or groundbased observations on the massive end of the mass function, as well as the merger rate of member galaxies, might be underestimated.
2. We find that the spatial distribution of member galaxies in J1001 is highly concentrated.By performing NFW profiling fitting of the stellar mass and number density profiles, we extend previous measurements on the concentration of cluster member galaxy distribution to z ∼ 2.5.In contrast to the cosmic evolution of dark matter halos, we find that cluster galaxy distribution at higher redshifts exhibits a higher concentration.The central stellar density of J1001 is already comparable to or even higher than more massive low-redshift clusters, while the outskirts lie below that of low-redshift clusters.This strongly suggests an inside-out formation scenario, at least for these early formed clusters.
3. Based on the mass-complete sample of star-forming members, we show that the stellar mass function of J1001 shows a prominent "top-heavy" feature, with overabundant massive SF members compared to a single Schechter function.The total number and stellar density of these massive SFG members are comparable to that of massive SFG and QGs combined in clusters at lower redshifts.Together with the low quiescent fraction in J1001, these findings suggest a minimal role of preprocessing, and most of the QGs should be quenched only after they were accreted onto the cluster.
These findings provide novel insights into the formation processes of clusters and their member galaxies at z ∼ 2.5.Based on a complete census of member galaxies, we show that the galaxy density profile and stellar mass functions are powerful tools to characterize the evolutionary state of galaxy (proto)clusters.The centrally concentrated galaxy density profile and the "topheavy" stellar mass function of the SF members in J1001 indicate that these early formed clusters and their member galaxies grow in an inside-out and top-bottom fashion: most of the massive galaxies in the core are assembled first, and less massive galaxies in the outskirts are formed later.Future studies with much larger samples of (proto)clusters at z ∼ 2 − 4 will provide a more comprehensive understanding of clusters and their member galaxy formation at their peak formation epoch (Zhou et al. 2024).(2014); and the Spitzer/IRAC channel 1, 3, and 4 images from the Spitzer Large Area Survey with Hyper-Suprime-Cam (SPLASH) (Steinhardt et al. 2014).

B. THE METHODS OF PHOTOMETRY
During the construction of our multiwavelength photometric catalog, we combine different methods of photometry for different bands since our data cover a wide resolution range.Firstly, for the high-resolution images from HST/ACS and JWST, we run Source-Extractor v2.25.0 (Bertin & Arnouts 1996, 2010) in the dual image mode to obtain the Kron aperture photometry (Kron 1980; FLUX AUTO from Source-Extractor), during which the fluxes are measured on the convolved images with the resolution of the JWST/NIRCam F444W image.For the HST/WFC3 images whose resolutions are slightly larger than the F444W image, we get Kron aperture photometry without performing a PSF match on them.Instead, we further convolve the F150W image to match the PSF size of HST/WFC3 images and take the ratio of the F150W fluxes after and before this convolution f 150,F444W resolution /f 150,WFC3 resolution as the correction factor for WFC3 fluxes.We note that these direct results given by Kron aperture photometry can miss some fluxes in the outskirts of each source, especially since a nonnegligible fraction (∼ 10%) of light can be scattered to large radii by the PSF of JWST/NIRCam.In this case, to obtain the accurate total flux, we estimate the ratio of missed light by performing Kron aperture photometry on the PSF and correct the Kron fluxes based on this ratio.
For other images with poorer resolution than HST and JWST, we use T-PHOT v2.0 (Merlin et al. 2015(Merlin et al. , 2016) ) to obtain prior-based deblended photometry based on the source positions and morphologies detected in the JWST/NIRCam images.T-PHOT is a deconfusion photometry software that convolves the cutouts of each source from a high-resolution image to produce low-resolution priors and then fits the flux of the same source in a low-resolution image with these priors.For images from ground-based telescopes and Spitzer, we perform T-PHOT photometry using the image of the closest JWST/NIRCam band as high resolution prior to obtaining deblended total flux for each F277W-detected source.
When combining the results measured by different methods, Merlin et al. (2021) have shown that the results of T-PHOT photometry do not have systematic bias compared with Kron aperture photometry.In addition, we also perform a similar test with our data: We directly perform source detection and Kron aperture photometry on the K s -band image from UltraVISTA (McCracken et al. 2012) and cross-match our F277Wdetected catalog with this blind K s -band catalog.

C. EXAMPLES OF THE BEST-FIT SED
Two examples of the best-fit SED given by Bagpipes are shown in Figure C2.The first one is representative of massive members that dominate the top-heavy feature, while the second one is representative of small members around the limiting depth of this work.

D. SAMPLE COMPLETENESS
We present the derivation of selection completeness with Monte Carlo simulation in this appendix.Firstly, following Shimakawa et al. (2018b), selection completeness is considered as a function of narrowband magnitude.At each given narrowband magnitude, we embed 500 point sources into the real narrowband and K s -band image.The K s fluxes of these sources are determined by their narrowband magnitude according to the best-fit power-law relation shown in the lower panel of Figure D3 with a random 0.12 dex scatter.This tight relation be- tween the flux excess and narrowband magnitude is fitted with the 64 real HAEs from observation.Next, we perform T-PHOT photometry on the simulated narrowband and K s -band images, during which the simulated point sources are also embedded into the F277W image to provide the high-resolution priors.We then test the recovery rate of the membership selection procedure described in Section 3. Next, we also test the selection completeness as a function of stellar mass.To determine the K s band fluxes for sources at a given mass, we collect the galaxies at 2.4 < z < 2.6 within the ultradeep stripes of the Ultra-VISTA survey (McCracken et al. 2012) from the COS-MOS2020 catalog (Weaver et al. 2022), and fit the stellar mass-K s flux relation at each given mass with a Gaussian distribution.Then, the K s -band fluxes of simulated sources can be randomly generated using this best-fit Gaussian distribution.As for narrowband fluxes, we firstly generate the SFR for each source based on its stellar mass by taking the SF main sequence at z = 2.5 (Popesso et al. 2023) with a 0.3 dex scatter.This SFR can be converted to a Hα flux based on the Equation 2of Kennicutt (1998) with a reduction of 0.24 dex for IMF calibration.Since the redshift of J1001 is close to the edge given by the NB2300 filter (Figure 1), we randomly give the redshift of each simulated source using a Gaussian distribution given by fitting the distribution of available z spec .We then calculate the wavelength of the Hα line at this redshift and multiply the Hα flux with the relative transmission of the NB2300 filter at this wavelength.Next, in order to obtain the narrowband flux of each simulated source, we adopt the dust attenuation-stellar mass relation (Garn & Best 2010;Shapley et al. 2022) with a 0.4 dex scatter, the contamination from [N II ] ( Shimakawa et al. 2018b), and a color-term variation 0.1.Once the K s -band and narrowband fluxes of every simulated source are assigned, we can follow the steps from the last paragraph to test and fit the recovery rate at each stellar mass.The best-fit result is also shown in Figure D3 (upper) and can be written as (D2) As a robustness test, the SMFs of SFGs with two different incompleteness-correction procedures are compared in Figure D4.We find the top-heavy shape of SMF is robust under different incompleteness-correction procedures.These two SMFs are almost identical at the high-mass end (M * > 10 10 M ⊙ ) and only show a small difference around the limiting depth (M * ∼ 10 9 M ⊙ ) without affecting the top-heavy shape.Since the bestfit parameters in Equation D2 might be influenced by the uncertainty (e.g.environmental dependence and the bias from observation) of the stellar mass-K s flux relation, SF main sequence, and dust attenuation, we use Equation D1 for the incompleteness correction in Section 4. We caution that the completeness discussed here is targeted at the sample of HAEs.This tends to slightly underestimate the completeness of the full SF sample since we also include four additional non-HAE spectroscopic members (mostly at M * = 10 10.0−10.5 M ⊙ ).

E. THE SCHECHTER FITTING OF SMFS
In order to avoid the uncertainty caused by the arbitrary binning procedure, the SMFs are fitted with maximum-likelihood estimation (MLE).We perform the MLE fitting by minimizing the negative log likelihood (e.g.Marshall et al. 1983;Hill et al. 2022)  sidered as the mixture of two single Schechter functions, which is similar to the mixture of gamma distributions (Young et al. 2019), we use the Expectation Maximization algorithm (Dempster et al. 1977) to determine the contribution of each member to each Schechter component.In addition, when we are fitting the SMF of SFGs, the power-law slope of the double Schechter function is fixed at α 1 = −1.5 (e.g.Shimakawa et al. 2018a,b) because it cannot be well constrained by our limited mass range and sample size.To obtain the uncertainty of the best-fit parameters, we randomly determine the stellar mass of each galaxy according to the mass uncertainty given by Bagpipes, as well as its weight considering the Poisson error.We repeat this procedure to perform the MLE for 100 times and then take the standard deviation of the 100 best-fit values as the uncertainty.All of the best-fit parameters and their uncertainties are listed in Table E1.

Figure 1 .
Figure1.Upper panel: the relative transmission curve of the NB2300 filter.The redshift range of HAEs is given by the FWHM range of this filter.Lower panel: the colormagnitude diagram used to select HAEs.The red stars, blue circles, and gray points are HAEs, spectroscopically confirmed members and interlopers, respectively.The magenta solid line shows the limitation given by Equation 1 with Σ = 2.The magenta dashed line corresponds to the color cut given by Equation2.

Figure 2 .
Figure 2. Left: the red, green, and blue composite color image of J1001.The R, G, and B channels of the image correspond respectively to the NIRCam/F444W, F277W, and F115W bands.The large blue square indicates the region covered by Subaru/MOIRCS narrowband imaging targeting HAEs.All member galaxies are marked with open circles, with massive HAEs, less massive HAEs, non-HAE SF members, and quiescent members denoted in yellow, cyan, green, and red, respectively.Right: JWST/NIRCam RGB and HST/F160W images of the four HST-dark red members (circles) and their companions (crosses).
using the "curve fit" algorithm from Scipy v1.8.1 (Gommers et al. 2022).Data points outside the virial radius R 200 =369 +61 −53 kpc (Wang et al. 2016) are not considered since they can be affected by subhalo structures.According to this NFW fitting, the concentration of J1001 is c mass (≡ R 200 /r s ) = 93 +∞ −40 for stellar mass density profile and c number = 37 +∞ −27 for the number density profile, where +∞ represents a single power law.

Figure 3 .
Figure3.Left panel: the projected stellar mass density profile of J1001 compared with the best-fit projected NFW profiles at lower redshifts given by van derBurg et al. (2014Burg et al. ( , 2015)).Following van derBurg et al. (2015), the original results with dimensionless units are converted to have physical units by assuming M200 = 3 × 10 14 M⊙ at z ∼ 1 and M200 = 9 × 10 14 M⊙ at z ∼ 0.15.Field level is given by galaxies with 2.483 < z phot < 2.519 from the COSMOS2020 catalog(Weaver et al. 2022).The vertical dotted line shows the R200 of J1001.Middle panel: similar to the left panel, but the number density profiles are shown.Right panel: The redshift evolution of the concentration of cluster-galaxy number density profiles from observations (filled circles) and of the density profile of massive (M h = 10 14 M⊙) dark matter halos (dotted line) based on simulations fromDuffy et al. (2008).

Figure 4 .
Figure 4.SMF of cluster J1001 compared with results in other protoclusters at z ∼ 2. The blue and red solid lines show the SMF of SFGs and QGs in cluster J1001 at z = 2.506.The blue dashed line is the SMF of HAEs in J1001, which provides fair comparisons withShimakawa  et al. (2018a,b).The purple dashed line is the stacked SMF of SFGs in 14 (proto)clusters between z = 2.0 − 2.5 in the COSMOS field(Edward et al. 2023), the green dotted line is the SMF of HAEs in cluster USS 1558 at z = 2.5(Shimakawa et al. 2018a), the brown dashed-dotted line is the SMF of HAEs in cluster PKS 1138 at z = 2.2(Shimakawa et al. 2018b).All the three SMFs of other (proto)clusters have been renormalized to match the normalization of J1001.
APPENDIX A. ARCHIVAL IMAGES In this work, the archival multiband images from U to IRAC channel 4 of the COSMOS field are used, including the U band from The COSMOS-WIRCam Near-Infrared Imaging Survey taken with the Canada-France-Hawaii Telescope (CFHT) (McCracken et al. 2010); B and IB427 from the COSMOS-20 survey taken with Subaru/Suprime-Cam (Taniguchi et al. 2015); g, r, i, and z from the third public data release (PDR3) of the Hyper Suprime-Cam Subaru Strategic Program (Subaru/HSC-SSP) (Aihara et al. 2022); the F814W band from the Cosmic Evolution Survey taken with the Advanced Camera for Surveys on HST (HST/ACS; Scoville et al. 2007); Y , J, H, and K s band from the fourth data release of the UltraVISTA nearinfrared imaging survey (McCracken et al. 2012); an archival HST/WFC3 F110W image from Negrello et al.
shows the comparison of the K s fluxes measured by T-PHOT and Kron apertures.We find that the two fluxes are consistent for the unblended sources with a median deviation −0.9%.Meanwhile, the T-PHOT fluxes of blended sources can be much smaller and more reliable since the Kron results can be contaminated by nearby sources.

Figure B1 .
Figure B1.The comparison of the Ks-band fluxes measured with T-PHOT and the Kron method.The blue and red points show the unblended and blended sources, respectively.
2, and fit the resulting selection completeness with a complementary error function.The best-fit result is C(NB) = 0.488 × Erfc[1.372× (NB − 23.565)] (D1) This best-fit error function used for completeness correction is shown in Figure D3 (upper).
where M * ,min = 10 9.2 M ⊙ and M * ,max = 10 12.0 M ⊙ are the lower and upper limits of the fitting range, respectively, and C(NB) is the incompleteness-correction factor given by Equation D1 in Appendix D. During the fitting, since the double Schechter function can be con-

Figure C2 .
Figure C2.Left panels: examples of the best-fit SED of two member galaxies from Bagpipes.Right panels: the probability density distribution functions of the photometric redshift from EAZY and stellar mass from Bagpipes.

Figure D3 .
Figure D3.Upper panel: the selection completeness.Blue and red points are the simulated selection completeness as a function of narrowband magnitude and stellar mass, respectively.And the blue and red lines are the completeness function fitted by the complementary error function.Lower panel: the narrowband exceeds fluxes as a function of narrowband magnitude.The red solid line is the result of linear fitting, while the dotted lines show the scatter of 0.12 dex.

Table E1 .
Best-fit parameters of the stellar mass functions.Note.-(a) ∆(BIC) ≡ BIC single − BIC double , ∆(BIC) > 0 means that the double Schechter function can describe the sample better than the single Schechter function, and the typical threshold for a valid comparison is ∆(BIC) > 8 ∼ 10 (Kass & Raftery 1995).(b) Fixed parameters are enclosed in square brackets.