A massive quiescent galaxy confirmed in a protocluster at z=3.09

We report a massive quiescent galaxy at $z_{\rm spec}=3.0922^{+0.008}_{-0.004}$ spectroscopically confirmed at a protocluster in the SSA22 field by detecting the Balmer and Ca {\footnotesize II} absorption features with multi-object spectrometer for infrared exploration (MOSFIRE) on the Keck I telescope. This is the most distant quiescent galaxy confirmed in a protocluster to date. We fit the optical to mid-infrared photometry and spectrum simultaneously with spectral energy distribution (SED) models of parametric and nonparametric star formation histories (SFH). Both models fit the observed SED well and confirm that this object is a massive quiescent galaxy with the stellar mass of $\log(\rm M_{\star}/M_{\odot}) = 11.26^{+0.03}_{-0.04}$ and $11.54^{+0.03}_{-0.00}$, and star formation rate of $\rm SFR/M_{\odot}~yr^{-1}<0.3$ and $=0.01^{+0.03}_{-0.01}$ for parametric and nonparametric models, respectively. The SFH from the former modeling is described as an instantaneous starburst while that of the latter modeling is longer-lived but both models agree with a sudden quenching of the star formation at $\sim0.6$ Gyr ago. This massive quiescent galaxy is confirmed in an extremely dense group of galaxies predicted as a progenitor of a brightest cluster galaxy formed via multiple mergers in cosmological numerical simulations. We newly find three plausible [O III]$\lambda$5007 emitters at $3.0791\leq z_{\rm spec}\leq3.0833$ happened to be detected around the target. Two of them just between the target and its nearest massive galaxy are possible evidence of their interactions. They suggest the future strong size and stellar mass evolution of this massive quiescent galaxy via mergers.


INTRODUCTION
Corresponding author: Mariko Kubo kubo@cosmos.phys.sci.ehime-u.ac.jp The deep multi-wavelength imaging surveys have uncovered wide variety in galaxy population at high redshift. Lately, massive quiescent galaxies at up to z ∼ 4 have been discovered and confirmed spectroscopically (e.g., Daddi et al. 2005;van Dokkum et al. 2008;Kriek et al. 2009;Tanaka et al. 2019;Valentino et al. 2020;. They are plausible progenitors of giant ellipticals today though their formation mechanism is a challenging problem since they have to form a stellar mass M ⋆ 10 11 M ⊙ and quench star formation in the early universe. Massive quiescent galaxies at high redshift are remarkably compact (a few to ten times smaller than local giant ellipticals) in general (e.g., Daddi et al. 2005;van Dokkum et al. 2008;Kriek et al. 2009;Kubo et al. 2018;Lustig et al. 2021). To form them, the star formation accompanied with centrally concentrated intense starburst like gas rich major merger and violent disc instability is required (e.g., Dekel et al. 2009;Burkert et al. 2010;Dekel & Burkert 2014;Zolotov et al. 2015). Lately, submillimeter galaxies (SMGs) similar in sizes, masses, and velocity dispersion with compact massive quiescent galaxies have been discovered and support this picture (e.g., Toft et al. 2014;Spilker et al. 2016;Valentino et al. 2020). It should also be considered that galaxies assembled in the early universe or when the universe was much denser tend to be compact (e.g., Williams et al. 2014;Wellons et al. 2015Wellons et al. , 2016Akhshik et al. 2021). Such massive quiescent galaxies at high redshift need strong size evolutions where evolutions of mass and size via mergers (e.g., Bezanson et al. 2009;Naab et al. 2009), expansion due to mass-loss driven by feedback (e.g., Fan et al. 2008Fan et al. , 2010, and/or evolution of stellar population ) scenarios have been proposed but not yet been constrained exactly.
Given the environmental dependance of galaxies today, we need to investigate the formation and evolution scenarios together with environments. However, massive quiescent galaxies at z > 2 in previous studies are mostly found in general fields or their environments are poorly explored, though several studies reported possible overdensities around them (Strazzullo et al. 2015;Belli et al. 2017;Schreiber et al. , 2020. To track the formation history of giant ellipticals generally in clusters of galaxies properly, we also need to study their progenitors in protoclusters. The appearance of red sequence in protoclusters at z 3 have been found by deep near-infrared (NIR) observations (e.g., Kodama et al. 2007;Zirm et al. 2008;Doherty et al. 2010;Kubo et al. 2013;Strazzullo et al. 2016;Shi et al. 2020), and massive quiescent galaxies have been confirmed spectroscopically by detecting Balmer absorption features in clusters at up to z ≈ 2 . The environmental difference of galaxy population at z ∼ 4 is shown statistically by using the deep and wide opti-cal survey of the Hyper Suprime-Cam Subaru Strategic Program survey (Toshikawa et al. 2018;Kubo et al. 2019;Ito et al. 2020). However morphologies of only a small number of massive quiescent galaxies have been studied (Zirm et al. 2008;Kubo et al. 2017, and also Strazzullo et al. 2013 andMei et al. 2015 including mature clusters at z ∼ 2), and the detailed spectral characteristics of massive quiescent galaxies have not yet been characterized in protoclusters.
We have studied a protocluster at z = 3.09 in the SSA22 field (Steidel et al. 1998) which is known as one of the most significant structures at high redshift. This structure was first identified by the overdensity of optically selected galaxies (Steidel et al. 1998;Hayashino et al. 2004;Matsuda et al. 2005;Yamada et al. 2012) and in later, further characterized with the overdensity of active galactic nuclei (AGNs) selected in X-ray (Lehmer et al. 2009a,b), SMGs (Tamura et al. 2009;Umehata et al. 2014Umehata et al. , 2015Umehata et al. , 2017Umehata et al. , 2018, and dusty starburst and passively evolving galaxies selected photometrically based on the deep NIR imaging (Kubo et al. 2013). Thus this protocluster is an ideal laboratory used to witness the transition of starburst galaxies into quiescent galaxies. Particularly, the 2 ′ x 3 ′ region containing the brightest SMG of this field was deeply observed with Atacama large millimeter/submillimeter array (ALMA) and a significant overdensity of SMGs was discovered (Umehata et al. 2015(Umehata et al. , 2017(Umehata et al. , 2018hereafter ADF22). They are clustering along the large scale (∼ 1 Mpc) Lyα filament which indicates the supply of cold gas through cosmic web (Umehata et al. 2019). This field is plausibly the central region of the SSA22 protocluster also supported from the large scale redshift distribution of galaxies (Matsuda et al. 2005;Kubo et al. 2015).
In this paper, we report the confirmation of a massive quiescent galaxy at the ADF22 field of the SSA22 protocluster by detecting the absorption features spectroscopically with multi-object spectrometer for infrared exploration (MOSFIRE; McLean et al. 2012) on the Keck I telescope. This paper is organized as follows: In Section 2, we describe the target and observation with MOS-FIRE and data analysis. In Section 3, we present the spectral energy distribution (SED) fittings and discoveries of new [O III] λ5007 emitters around this massive quiescent galaxy. In Section 4, we discuss the star formation history and future evolution scenario. Section 5 is the conclusion. In this study, we adopt cosmological parameters Ω m = 0.3, Ω Λ = 0.7 and H 0 = 70 km s −1 Mpc −1 . We assume the Chabrier (2003) Initial Mass Function (IMF). Magnitudes are expressed in the AB system.

Target
Previously, we conducted deep and wide NIR imaging observations of the SSA22 protocluster with multiobject infrared camera and spectrograph (MOIRCS; Ichikawa et al. 2006;Suzuki et al. 2008) on the Subaru telescope over a field of ≈ 112 arcmin 2 for K s ≈ 24 at 5σ level and found an overdensity of massive galaxies based on the photometric redshifts (z phot ) estimated with the SED fitting of the optical to MIR photometry (Kubo et al. 2013). We further confirmed many of these candidates as protocluster members spectroscopically (Kubo et al. 2015). Among them, we selected candidate passively evolving galaxies based on the i − K vs. K − [4.5] color following Labbé et al. (2005). These colors are similar to rest-frame U V J often used to select quiescent galaxies (e.g., Whitaker et al. 2013). We selected galaxies satisfying i − K > 3 and K − [4.5] < 0.5, and 2.6 < z phot < 3.6 as candidate quiescent galaxies. The color criterion in Kubo et al. (2013) was set to avoid the contamination of dusty starburst galaxies detected in 24 µm which distribute at the border of general color criterion for quiescent galaxies.
Our target (R.A., Dec = 22:17:37.25, +00:18:16.0, hereafter ADF22-QG1) satisfies these selection criteria. It has a total magnitude of K s,tot = 21.55. The SED is well fitted with that of a quiescent galaxy namely single burst or exponentially declining star formation history (SFH) where star formation rate (SFR) ∝ exp(−t/τ ) with τ ∼ 0.1 Gyr, and with an age of 1 Gyr at z ∼ 3.1. We show the target on a rest-frame U V J color diagram in Fig. A.1. Its spectroscopic redshift is limited to z spec = 2.7 or 3.0 − 3.15 by detecting the 4000 A break (Kubo et al. 2015). We took a high-resolution K ′ -band image of the target by using the adaptive optics AO188 and infrared camera and spectrograph (IRCS) on the Subaru telescope and found that its effective radius (r ef f ) is 1.01 ± 0.04 kpc (Kubo et al. 2017). Uniquely, it is likely in an extremely dense group of massive galaxies and SMGs identified by follow-up NIR spectroscopic observations of a 1.1 mm source found with ASTE/AzTEC (Kubo et al. 2016; here after AzTEC14 group). The left panel of Fig. 1 shows the MOIRCS K s -band image (described in § 2.3) for 20 arcsec x 30 arcsec around the target. Except for our target, there are nine galaxies at 3.0774 ≤ z spec ≤ 3.0926 confirmed by detecting the [O III] λ5007 and/or CO(3-2) emission lines (Kubo et al. 2016;Umehata et al. 2018). Assuming that they are hosted by the same halo, its expected halo mass is 10 13 M ⊙ based on the velocity distribution (Kubo et al. 2016).

Observation with MOSFIRE
We observed the target with MOSFIRE in K-band in September 2017 and 2020, and in H-band in September 2020 . The right panel of Fig. 1 shows the slit positions. The observations and data reduction are described in Umehata et al. (2019) (K-17) and in H. . Briefly, the seeing conditions were 0.7 ∼ 0.9, 0.6 and 0.6 arcsec for K-17, K-20 and H-20 runs. We adopted the 2-position mask-nod sequence with 1.5 arcsec dithers along slits with 120 (180) sec for each exposure in H(K)-band. The total exposure times at ADF22-QG1 were 2.5 and 6.5 h (3.2 h for K-17 and 3.3 h for K-20) in H and K-band, respectively. The spectra were reduced with the publicly available MOS-FIRE data reduction pipeline, MCSDRP (Steidel et al. 2014).
We extracted the one-dimensional spectrum as follows to measure the fluxes of the spectrum and photometry consistently. Here we calibrated the spectra to the H and K s -band photometry taken with MOIRCS (Section 2.4). First, H-band spectrum was re-binned to match the pixel scale with that in K-band (≈ 0 ′′ .18/pix and 2.17Å/pix) by using the magnify task of IRAF. Next, we combined the target and error spectra for ≈ 2.0 arcsec in the spatial direction by an inverse variance weighting. Then, we measured the H and K s -band fluxes of the spectra by applying the transmission curve of MOIRCS H and K s -band filters. Finally, we corrected the spectra to match these fluxes with the photometric fluxes.

Deeper K s -band image taken with MOIRCS
We newly obtained a deeper K s -band image with MOIRCS which was upgraded in 2015 (nuMOIRCS; Walawender et al. 2016;Fabricius et al. 2016) on the Subaru telescope in July 2020 (PI: H. Umehata). We collected the images with PSF FWHM of 0.3 − 0.6 arcsec for 6.6 h total exposure time. The data is reduced with the MOIRCS data reduction pipeline MCSRED 1 . We reduce the data in a standard procedure but the flux is calibrated to our previous K s -band image taken with MOIRCS (Kubo et al. 2013) which is calibrated to UKIRT Infrared Deep Sky Survey (UKIDSS). The zeropoint error is ≈ 0.05 mag relative to UKIDSS. Briefly, the final combined image has a PSF FWHM of 0.40 arcsec and the 5σ limiting magnitude at 0.80 arcsec diameter of 25.34 mag.  , the spectroscopic redshift is estimated with SLINEFIT 2 which fits a spectrum with combinations of a stellar continuum template and emission/absorption lines to find a redshift and line properties by a χ 2 minimization procedure. It returns a spectroscopic redshift z spec = 3.0922 +0.0008 −0.0004 which is consistent with the absorption features. Thus ADF22-QG1 is certainly in the SSA22 protocluster, moreover a member of the dense group of galaxies discussed in Section 4.2. This is the most distant quiescent galaxy detected of Balmer absorption features in a protocluster. Hereafter we adopt this spectroscopic redshift for ADF22-QG1.

Spectroscopic redshift
The [O II] flux of ADF22-QG1 is F [OII] = 5.7 ± 1.0 × 10 −18 erg s −1 cm −2 . The SFR from [O II] luminosity is calculated as follows; First, the SFR from Hβ luminosity is computed assuming the case B recombination value intrinsic Hα/Hβ = 2.86 (Osterbrock & Ferland 2006), and the SFR to Hα luminosity relation in Kennicutt & Evans (2012), The normal star forming galaxies (SFGs) at z ∼ 3 have log [OII]/Hβ = 0.0 − 0.6 at most according to Onodera et al. (2016). Then the SFR from [OII] is computed as, .5 which are more than ten times higher than that from the SED fitting described in later. [O II] emission lines are frequently seen in quiescent galaxies at high redshift and the SFRs from [O II] are tend to be higher than those from SED modeling (e.g., Yan et al. 2006;Lemaux et al. 2010;Belli et al. 2015; Lamastra et al. (2009), it is no wonder that ADF22-QG1 is not detected with the Chandra observation of the SSA22 protocluster with a sensitivity limit L 2−10keV ≈ 5.7 × 10 42 erg s −1 in Lehmer et al. (2009b). From these above, the [O II] emission of ADF22-QG1 may originate in a partly remaining star formation or weak AGN. Empirically, the excess [O II] of massive quiescent galaxies at high redshift is likely due to AGNs or LINERs rather than star formations (e.g., Yan et al. 2006;Lemaux et al. 2010).

SED fitting
The absorption features together with multiwavelength photometry are strong constraints on SFHs of massive quenched galaxies as demonstrated in previous studies (e.g., Belli et al. 2019;Valentino et al. 2020;Akhshik et al. 2021). Here we fit the spectrum and photometry of ADF22-QG1 simultaneously with the SED models with parametric SFHs by using FAST++ 3 , and nonpara-

Data
First, we prepare the input photometric and spectroscopic data for SED modelings.
We use the u ⋆ BV Ri ′ z ′ JHK s , 3.6, 4.5, 5.8 & 8.0 µm-band photometry measured in Kubo et al. (2013). Briefly, we con-volve the u ⋆ to K s -band images to match the PSF to a FWHM of ≈ 1.0 arcsec, and measure fluxes with a 2.0 arcsec diameter aperture. To match with them, the IRAC 3.6 − 8.0 µm photometry is applied aperture correction computed by using the K s -band image (see detail in Kubo et al. 2013). We also use the flux measured on the archival F 814W -band image taken with Advanced Camera and Spectrograph on Hubble Space Telescope (PID9760), with the same manner, and IR luminosity (at 8 − 1000 µm; L IR ) limit based on the 1.2 mm (256.98 GHz) image taken with ALMA in Cycle-2 (Umehata et al. 2018) and Cycle-5 (PID. 2017.1.01332.S, PI. H. Umehata, H. Umehata et al in prep). The L IR is calculated by assuming the average 1.2 mm flux to L IR relation for the SED library in Danielson et al. (2017). ADF22-QG1 is not detected at 1.2 mm and the 3σ limiting flux is 75 µJy corresponding to L IR ∼ 0.9 − 2.0 × 10 11 L ⊙ taking the 95% confidence interval. The input observational data for Prospector is the same as that for FAST++ but the former uses L IR while the latter uses the 1.2 mm flux.
We correct the aperture photometry and onedimensional spectrum by multiplying them with 1.21 which is the ratio of the total (Kron) flux measured on the original image and aperture flux measured on the PSF matched image in K s -band. We also correct the Galactic extinction E(B − V ) = 0.053 found from the dust extinction finding tool at NASA/IPAC INFRARED SCIENCE ARCHIVE 5 based on Schlafly & Finkbeiner (2011). The one-dimensional spectrum is re-binned with six pixels (≈ 13Å) by taking the inverse variance weighting average. The [O II] flux estimated with SLINEFIT is subtracted from the H-band photometry. The [O II] and strong OH airglow on the spectrum are excluded from the SED fitting. We use the spectrum at 1.51 − 1.80 and 1.97 − 2.19 µm avoiding the spectrum with low transmission due to the instrument and/or sky. Both transmissions are still good but we avoid using the spectrum at > 2.19 µm because the non-uniformity of the sky transmission significantly remains on the reduced spectrum.

FAST++
First, we fit the spectrum and photometry simultaneously with SED models by using FAST++ following . We fit the SED with Bruzual & Charlot (2003) stellar population synthesis models adopting the Chabrier (2003) IMF, Calzetti (2001) dust attenuation law and the solar metallicity (Z = 0.02). We fit the models with extinction value A V = 0 − 6 mag with steps of 0.1 mag and log(age yr −1 ) = 7 − 9.3 with steps of 0.05 dex. We test SFH of exponentially declining (SFR ∝ exp(−t/τ )), delayed exponentially declining (∝ t × exp(−t/τ )), where t is time from the formation, and composite SFH in . This SFH is described as a combination of, and   Adopting M ⋆,FAST++ as an upper limit value.
A parametric SED modeling can fail to reproduce a true SFH if a galaxy has a more complex SFH, e.g., starburst, sudden quenching and rejuvenation. Then we also fit the SED with nonparametric models which can handle with complex SFH following Leja et al. (2017Leja et al. ( , 2019. The fitting is performed by using Prospector ) which uses the Flexible Stellar Population Synthesis (FSPS) code (Conroy et al. 2009). As in the FAST++ SED modeling, the Chabrier (2003) IMF and the Calzetti (2001) dust attenuation law are adopted. The two-component Charlot & Fall (2000) dust attenuation model is adopted in Prospector while the uniform screen model is adopted in FAST++. It consists from birth-clouds and diffuse dust screen components. The dust attenuation is parameterized with the optical depths at 5500Å of these two components, τ 1 (birth-clouds) andτ 2 (diffuse dust). In case of the Calzetti (2001) dust attenuation law,τ 1 is set to zero and the dust attenuation is controlled only withτ 2 . We adopt the result measured adopting the solar metallicity (Z = 0.02).  We use a continuity SFH prior which fits directly for ∆ log(SFR) between adjacent time bins. Here the adopted time bins are 0 < log(age/yr) < 7.5, 7.5 < log(age/yr) < 8.0, and ∆ log(age/yr) = 0.11 binning at the range 8.0 < log(age/yr) < 9.3 (14 time bins in total). The free parameters are log(M ⋆ /M ⊙ ) (between 7 to 12), ∆ log(SFR) between adjacent time bins (scale =0.3, df =2), andτ 2 (between 0.0 to 2.0). The sampling is performed with the nested sampler dynesty (Speagle 2020) and the maximum a posteriori (MAP) values are presented as the best-fit quantities.
Though we adopt the above model, we also test the continuity SFH model with free metallicity, and prospector-α  model. The latter model adopts a Dirichlet SFH prior, the prescription of dust attenuation law from Noll et al. (2009) which allows a flexible dust attenuation law slope andτ 1 toτ 2 ratio, and adds an AGN contribution. We summarize the results of these different modeling in Appendix B. Briefly, the results do not depend on the dust attenuation law and the contribution of AGN is negligible. The model with near solar metallicity is favored for the continuity-SFH model but low metallicity is favored for the Prospectorα (Dirichlet) model, though our data has no significant spectral indices to determine the metallicity definitely. A model with lower metallicity tends to have a SFH in which the star formation quenches earlier. Thus we here adopt the continuity-SFH model with the solar metallicity (Z = 0.02) as the conservative model. The red curves in Fig. 2 and 3 show the bestfit SED, and Table 1 lists the SED parameters and SFH quantities found with FAST++. We note that there is no significant difference in the SED quantities and SFHs derived with the three SFH models.
Here we adopt the best-fit model found with the com-posite SFH modeling to compare the results directly with literature Valentino et al. 2020;. The SFH and several SFH quantities, SF R main , t 50 and t q , are computed with FAST++. The SF R main is the mean SFR during the shortest time interval over which 68 % of the star formation took place. The formation time t 50 is the time at which 50% of the total stellar mass has been formed, excluding mass loss and recycling, and the quenching time t q is the elapsed time since the SFR dropped below 10% of the SF R main . Both t 50 and t q are given as a lookback time from the observed redshift. The uncertainties of the values from the SED modeling in Table 1 are the 90% confidence interval values.
The best-fit model and its reduced χ 2 value change marginally by setting different velocity dispersion. We also tried the penalized pixel-fitting algorithm (pPXF; Cappellari 2017) but could not also obtain a robust velocity dispersion estimate. It is because all the available Balmer absorption features are partly covered by OH airglow. Here we show the best-fit parameters for velocity dispersion = 320 km s −1 where the χ 2 value minimizes (χ 2 /ν = 2.7). We also test the models with metallicities 0.2, 0.4 and 2.5 Z ⊙ but it also do not change the results significantly. A corner plot for covariance among the age, A V , Z, SFR, M ⋆ , t 50 and t q is presented in Fig. B.2, which shows that the stellar population of ADF22-QG1 is well constrained with small degeneracy.
The L IR limit for ADF22-QG1 corresponds to SFR < 9 − 21 M ⊙ yr −1 using the L IR to SFR conversion factor in Kennicutt & Evans (2012). Note that the L IR and thus SFR IR based on a single-band flux depend greatly on the adopted SED models (Santini et al. 2019;Mizukoshi et al. 2021). Scaling the 1.2 mm flux limit of our target to the model prediction in Santini et al. (2019), it can have a SFR IR upper limit 100 M ⊙ yr −1 . Even though we adopt the SFR [O II] or SFR IR at upper limit, ADF22-QG1 is classified as a quenched galaxy. We put the SFH evaluated with FAST++ in Fig. B.1 where the whole stellar mass is formed instantaneously. The formation time t 50 of ADF22-QG1 is younger (0.62 +0.09 −0.05 Gyr) than that estimated from the fitting only with the photometry ( 1 Gyr). This is expected from the spectral characteristic of ADF22-QG1 that 4000 A break is strong but Balmer absorption features are still significant. Such a SFH cannot be measured robustly without deep NIR spectroscopic observations like this study. Figure 4 compares the SED quantities of ADF22-QG1, and quenched galaxies at 2.8 < z spec < 4.0 in literature Valentino et al. 2020;. Star formation of ADF22-QG1 is suppressed well among them. ADF22-QG1 has formed later but more rapidly than other QGs. At this point, we find no significant environmental dependance in the SED quantities. Note that the figure looks lacking slowly quenched (t 50 − t q > 0.1 Gyr) galaxies however this is likely the sample bias that quiescent galaxies at z > 3 quenched earlier than the current sample are hardly surveyed completely at the current survey depth K s 24 mag ).

SED: nonparametric
The best-fit SED from nonparametric modeling of SFH is shown with blue curves in Fig. 2 and 3. It also fits the observed SED well. The SED quantities from nonparametric modeling are denoted with continuity in Table 1. Fig. B.4 shows a corner plot for covariance among the fitting parameters in the nonparametric modeling.
Both the nonparametric and parametric modeling fit the observed SED with a model of massive galaxy with suppressed star formation (sSFR continuity = 0.7 +10.3 −0.7 × 10 −14 yr −1 ). However, they find different SFHs and stellar masses. Fig. 5 shows the SFH evaluated with nonparametric modeling. It is flatter than that from FAST++ and does not take an extreme value. The SFR drops sharply after the time bin 8.80 < log(age/yr) < 8.91 or the star formation is quenched during this time bin. It agrees very well with the quenching time log(t q,FAST++ /yr) = 8.78 +0.04 −0.20 from FAST++. Because of the larger contribution from old stars, the stellar mass from nonparametric modeling is twice larger than that from FAST++.

Newly detected [O III] emitters
Adjacent to ADF22-QG1, emission lines are detected at three positions between 2.042 and 2.045 µm (Fig.  6, here after c1, c2 and c3). We extract their onedimensional spectra by combining the spectra for ≈ 1.3 and 2.2 arcsec in the spatial direction for c1, and c2 & c3, respectively, and smoothing them with three pixels in the wavelength direction. We fit the one-dimensional spectra with Gaussian profiles by a standard χ 2 minimization procedure. We derive them under the assumption that they are [O III]λ5007.
We estimate their sky positions based on the average spatial profiles in a range ≈ 17Å (8 pixels) at emission lines ( Fig. 6 right tops). In the case of c1 and c3, they show clear peaks at 1.6 and 2.3 arcsec from ADF22-QG1. Wavelength ( The peak of c2 is not so clear but likely fallen between c3 and ADF22-QG1. Their K s -band magnitudes measured with 0.80 arcsec diameter apertures are < 26.34 (< 2σ), 26.02 (2.7σ) and 25.72 (3.5σ), respectively. Table 2 summarizes their observed properties. Interestingly, the spatial locations and redshifts of c2 and c3 are just between ADF22-QG1 and the nearest massive galaxy, K15d at z spec = 3.0774 ± 0.0003 (Kubo et al. 2016).

SFH
Both the parametric and nonparametric modeling of the SFH find that ADF22-QG1 is rapidly quenched at ∼ 0.6 Gyr ago. To reproduce such a sudden quenching, strong feedbacks are needed (e.g., Belli et al. 2019;Rodríguez Montero et al. 2019;. If the SFH of ADF22-QG1 is alike with that estimated with FAST++, feedback from young massive stars is applicable for the quenching (Murray et al. 2005). The quenching via AGN feedback cannot be ignored especially for massive galaxies. ADF22-QG1 is not detected with Chandra but its [O II] can originate in a dying AGN. We will further discuss the roles of AGNs for the evolution of protocluster galaxies in an upcoming paper. It is beyond the scope of this paper but also a problem to be solved in the future that how ADF22-QG1 has maintained its quiescence for several 100 Myr although it is in a gas rich environment and surrounded by star-bursting neighbors (Umehata et al. 2019) which can induce further starburst.
Which SFH is more realistic? The SF R main,FAST++ of ADF22-QG1 is not so realistic because it is even higher than those of the brightest SMGs with SFR 5000 M ⊙ yr −1 (Riechers et al. 2013;Oteo et al. 2016;Spilker et al. 2020). Theoretically, if the SFR surface density of a SFG reaches the Eddington limit where momentum-driven wind induced by radiation pressure on dust heated by young massive stars, its star formation can be suppressed (Murray et al. 2005). The limiting SFR surface density for SMGs is ∼ 1000 M ⊙ yr −1 kpc −2 based on Andrews & Thompson (2011), and SMGs do not exceed this limit in general (Simpson et al. 2015) while the brightest SMGs are near this limit. The SFR surface density of ADF22-QG1 at SF R main,FAST++ computed by dividing it with an area with a radius of the twice r ef f at the observed redshift is 0.86 +1.67 −0.78 × 10 3 M ⊙ yr −1 kpc −2 which by far exceeds the Eddington limit SFR surface density.
Thus the SFH from nonparametric modeling is more realistic. The SFR of ADF22-QG1 before quenching is consistent with or lower than those of the observed SMGs in protoclusters at z ∼ 4 (e.g., Oteo et al. 2018;Pavesi et al. 2018;Miller et al. 2018;Long et al. 2020) and the SSA22 protocluster itself (Umehata et al. 2018). Our argument can also be true for the other massive quiescent galaxies at z > 3 in literature in which parametric SED modelings are widely used. In these studies, the SFHs are often described as a vigorous starburst similar to ADF22-QG1 computed with FAST++ (e.g., Valentino et al. 2020;. Then the most active starburst galaxies like SMGs have SFRs in concordance with the SFH for massive quiescent galaxies. However, to explain the observed properties, e.g., number densities, less-bursty SFGs on the star formation main sequence are also needed to be the major progenitors of massive quiescent galaxies (Barro et al. 2013(Barro et al. , 2017Popping et al. 2017;Gómez-Guijarro et al. 2019;Valentino et al. 2020). In addition, observed stellar masses of the brightest SMGs at z > 3 are often larger  Figure 7. The line-of-sight velocities and spatial distances of the galaxies relative to the ADF22-QG1. The black points are based on the redshifts of galaxies listed in Kubo et al. (2016) and Umehata et al. (2019) while the red points show the [O III] emitters confirmed in this study. The blue solid curve show the escape velocity for a NFW halo with a halo mass 1.1 × 10 13 M⊙ (virial radius = 170 kpc and concentration c = 5). The blue dashed curve is that scaled to projected one (the velocity and distance are divided with √ 3 and 3/2, respectively). than that of massive quiescent galaxies at 3 < z < 4 (Valentino et al. 2020). Adopting the SFH from nonparametric modeling, massive SMGs and main sequence SFGs are allowed as a progenitor of ADF22-QG1. Although we need a further robust SFH measurement in the future, application of SFHs more complex than typical parametric SFHs can correct the relation between massive quiescent galaxies and high-z SFGs.

Size and mass evolution
Here we confirm a compact massive quiescent galaxy as a secure progenitor of a giant elliptical or brightest cluster galaxy (BCG) in a cluster of galaxies today. Similarly, they reported that massive quiescent galaxies in proclusters/clusters at z 2 are more com-pact than nearby giant ellipticals in Zirm et al. (2008); Strazzullo et al. (2013); Mei et al. (2015). It needs a strong size evolution to evolve into a typical giant elliptical or BCG. Cosmological numerical simulations predict that a BCG is hierarchically formed via multiple mergers of galaxies (e.g., De Lucia & Blaizot 2007). Supporting such a scenario, ADF22-QG1 is not isolated but in a dense group of massive galaxies and SMGs. Fig. 7 shows the velocity distribution of galaxies within 100 physical kpc from ADF22-QG1. We use the redshifts derived from [O III] and/or CO(3-2) lines in Kubo et al. (2016) and Umehata et al. (2019) (black), and also newly confirmed [O III] emitters (red) in this study. The redshift of ADF22-QG1 is certainly close to the group members though it is offset from the median redshift z med = 3.087 of the group. Their velocity dispersion is σ v = 351 ± 52 km s −1 evaluated by a bootstrap resampling. According to the scaling relation based on N -body simulations in Evrard et al. (2008), where h(z) = H(z)/100 km s −1 and σ DM,15 is the normalization for a halo mass 10 15 h −1 M ⊙ . They found σ DM,15 = 1082.9 ± 4.0 km s −1 and α = 0.3361 ± 0.0026 for ΛCDM cosmology. Adopting the above σ v , the halo mass of AzTEC14 group is 1.1±0.4×10 13 M ⊙ consistent with our previous measurement (Kubo et al. 2016). The velocity offsets of all the group members are below the escape velocity for a halo characterized with the above halo mass and NFW (Navarro et al. 1997) profile. The velocity distribution and halo mass of AzTEC14 group are similar to those of SPT2349-56 which is known as an extreme overdensity of SMGs at z = 4.3 (Miller et al. 2018;Hill et al. 2020). Rotermund et al. (2021) shows that the stellar mass of the most massive member of 3.2 +2.3 −1.4 × 10 11 M ⊙ and the total stellar mass of (12.2 ± 2.8) × 10 11 M ⊙ at least for SPT2349-56. Correcting the total stellar mass of the K s -band detected galaxies of AzTEC14 group found in Kubo et al. (2016) with the newly measured stellar mass for ADF22-QG1 from the nonparametric SFH modeling, its total stellar mass is 6.7 +2.3 −0.8 × 10 11 M ⊙ at least, which is comparable or lower than SPT2349-56. On the other hand, SPT2349-56 is by far luminous at sub-mm than the SSA22 protocluster (Miller et al. 2018), i.e., the star formation in AzTEC14 group is less active. The SPT2349-56 is a close analog of AzTEC14 group but a larger stellar mass will be assembled in the former. Rennehan et al. (2020) simulated the assembly history of BCGs by performing isolated non-cosmological hydrodynamical simulations based on the observed properties of SPT2349-56. They predicted that all the group members at z = 4.3 merge into one massive galaxy by ∼ 0.5 Gyr. The spectroscopic confirmation of ADF22-QG1 in AzTEC14 group and newly confirmed [O III] emitters, which are possible evidence of interactions between ADF22-QG1 and K15d, further support such an early BCG assembly scenario via multiple mergers. We discussed the size and stellar mass growths of ADF22-QG1 via mergers based on the observed sizes and stellar masses of the AzTEC14 group members in Kubo et al. (2017). Assuming no further stellar mass and size growths in each group member, and all mergers are dry, the size and stellar mass of ADF22-QG1 can be 3 − 4 times and double from them at the observed redshift via mergers of all the group members, respectively. If ADF22-QG1 is just a progenitor of a giant elliptical, this scenario is enough. However, to grow ADF22-QG1 into a BCG, further size and stellar mass growths in each group member and/or further mergers of external galaxies are required. This argument does not change greatly by adopting the newly measured stellar mass for ADF22-QG1 in this study. Consideration of the newly confirmed SMGs in Umehata et al. (2019), the size and stellar mass growths especially for them, and stellar/AGN feedbacks in AzTEC14 group need further studies in the future. Finally we note that the candidate progenitors of the BCG of the SSA22 protocluster is not only ADF22-QG1 since there are several such extremely dense groups of galaxies (Kubo et al. 2016;Umehata et al. 2018).

CONCLUSION
We confirm a massive quiescent galaxy at the core of a protocluster at z = 3.09 in the SSA22 field. This is the most distant quiescent galaxy confirmed with Balmer absorption features in a protocluster as ever, and a securely selected giant elliptical/BCG progenitor. We fit the observed SED with both parametric and nonparametric models of SFHs. Both models agree with that our target is a massive galaxy with a suppressed star formation. The SFH found with the parametric modeling is described as a short starburst while that of the nonparametric modeling is a longer-lived SFH. The SFH found with the nonparametric modeling is more realistic given the extremely high SFR surface density at past derived from parametric modeling. On the other hand, both models support that the star formation is suddenly quenched after a starburst at ∼ 0.6 Gyr ago. To reproduce this, a strong feedback is required. This massive quiescent galaxy is confirmed as a member of an extremely dense group of massive galaxies and SMGs predicted as a progenitor of a BCG in cosmological numerical simulations. According to the simulations, such a system merge into one massive galaxy within ∼ 0.5 Gyr. We also newly confirm three plausible [O III] emitters around this quiescent galaxy. Two of them are possible evidence of the interaction between the quiescent galaxy and its nearest massive galaxy. They strongly support a hierarchical formation scenario of BCGs.
It is still unclear what is the major formation and quenching mechanism, and how to maintain the quiescence for several 100 Myr of our target in such a dense environment. Future studies of older and fainter quiescent galaxies, and stellar (gas and star formation) surface density of galaxies which may correlate with the SFH and quenching mechanism, and the roles of AGNs in protoclusters with James Webb Space Telescope (JW ST ) and Nancy Grace Roman Space Telescope (N GRST ) will enable us to discuss the typical formation scenario of cluster galaxies.

ACKNOWLEDGMENTS
We thank the anonymous referee for a number of useful suggestions. This work is supported by JSPS KAK-ENHI Grant Numbers 20K14530 (MK), 17K14252, and 20H01953 (HU). The spectroscopic data 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 observations were carried out within the framework of Subaru-Keck/Subaru-Gemini time exchange program which is operated by the National Astronomical Observatory of Japan. The K s -band image was collected with nuMOIRCS at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. We are honored and grateful for the opportunity of observing the Universe from Maunakea, which has the cultural, historical and natural significance in Hawaii. This paper makes use of the following ALMA data: Here we present the SFHs obtained with various models. First, we show the SFH and the corner plot for the main quantities evaluated with FAST++ in Fig. B.1 and Figure B.2, respectively. The corner plots presented in this paper are generated with corner.py (Foreman-Mackey 2016). In this model SFH, star formation occurs almost instantaneously. This model can explain the observed SED but results in unrealistically high maximum SFR.
We also present the nonparametric SFHs evaluated with Prospector. Here we show the SFHs evaluated adopting the models, i) Continuity SFH fixed with the solar metallicity (shown in Fig. 5 in main text), ii) Continuity SFH with free metallicity (-2.0 log Z/Z ⊙ < 0.2), iii) Prospector-α SFH fixed with the solar metallicity, and iv) Prospector-α SFH with free metallicity (-2.0 log Z/Z ⊙ < 0.5) (shown in Fig. B.3).
The continuity SFH adopted in i) and ii) models is described in section 2.5.1. The Prospector-α model is the model adopted in Leja et al. (2017). This model uses Dirichlet SFH prior, dust attenuation law prescription in Noll et al. (2009)  whereτ 2 controls the normalization of the diffuse dust, n is the diffuse dust attenuation index, k ′ (λ) is the Calzetti (2001) attenuation curve, and D(λ) is a Lorentzian-like Drude profile describing the UV dust bump. The strength of the UV dust bump is tied with the best-fit diffuse dust attenuation index following the results of Kriek & Conroy (2013). A flat prior over 0 <τ λ,2 < 4, a flat prior over 0 <τ λ,1 < 4, and a flat prior over −2.2 < n < 0.4 are adopted. This flexible dust attenuation law can handle dust attenuation law variation where Kriek & Conroy (2013) reported that galaxies with lower specific SFR tend to have larger dust attenuation indices. Here the AGN SED templates in Nenkova et al. (2008a,b) are adopted. The contribution of an AGN is controlled with the AGN luminosity as a fraction of the galaxy bolometric luminosity (f AGN ) and optical depth of AGN torus dust (τ AGN ). We show the SFH measured with ii) to iv) models in Fig. B.3. The corner plots of the parameters for ii) and iv) models are shown in Fig. B.4 and Fig. B.5, respectively. For the illustration purpose, the corner plots are shown for the modeling with the number of time bins N = 7 which give consistent parameters as those evaluated with 14 time bins. We also test Calzetti (2001) dust attenuation curve for Dirichlet prior but it does not change these parameters significantly. The AGN component in our target is negligible. From these above, the difference between the continuity SFH model and Prospector-α model mostly comes from the different SFH priors.
All i) to iv) models results in similarτ 2 and stellar mass. They are in agreement with the quiescence of the star formation in our target. In case of the continuity SFH prior, the modeling with free metallicity results in solar metallicity while that for Prospector-α results in log Z/Z ⊙ ∼ −1. As a result of the low luminosity, the SFH from iv) modeling quenched star formation earlier than those from other models. Since our observational data has no significant indices to confirm the metallicity robustly, we adopt the parameters measured adopting the continuity SFH prior and solar metallicity as conservative estimates.     Fig. B.4 but adopting the Prospector-α SFH prior. The parameters are log Z/Z⊙,τ2, z fraction which is proportional to sSFR between adjacent time bins (z fraction with n = 1 − 7), log(M⋆/M⊙), AGN fraction (fagn), the optical depth for AGN (agn tau), the ratio ofτ1 andτ2 (dust ratio), and n (dust index ).