Physical Properties of a Coma-analog Protocluster at z = 6.5

We present evidence for the discovery of a protocluster of starburst galaxies (Lyα emitters, or LAEs) near the end of the epoch of reionization. The recent trend in the search for high-redshift protoclusters focuses on utilizing bias tracers, such as luminous starburst galaxies, as signposts of overdensities. Thus, we conducted a photometric selection of LAE candidates around a pair of spatially close, luminous LAEs at z = 6.5 in the Subaru/XMM-Newton Deep Survey field, using OSIRIS in its imaging mode at the 10.4 m Gran Telescopio Canarias in La Palma, Spain. The spectroscopic follow-up was done with OSIRIS in its multiobject spectroscopy capability. We have spectroscopically confirmed 10 LAEs via their recognizable Lyα emission feature. The redshifts of these LAEs shed light on their 3D distributions within the observing window defined by the photometric selection. We have derived the galaxy number density contrast of δ gal = 3.18 − 1.99 + 3.47 , which led to the expected mass of the overdensity of 8.40 − 1.39 + 2.98 × 10 14 M ⊙ . We also found evidence for the presence of a virialized core with M 200 = 4.06 − 1.90 + 2.77 × 10 13 M ⊙ within this overdensity. Based on the extended Press–Schechter formalism, this overdensity would continue to grow in the linear regime and collapse to form a galaxy cluster at z coll = 0.84 − 0.43 + 0.57 . By the time this protocluster reaches z = 0, it will be a massive cluster of galaxies with mass 1.54 − 0.69 + 1.12 × 10 15 M ⊙ , comparable to the Coma cluster. Thus, our careful analysis has pointed to evidence that this protocluster would evolve into a Coma-analog cluster in the present-day universe.


Introduction
Observations of the most massive gravitationally bound structures such as galaxy clusters play a crucial role in our understanding of the large-scale structure formation of the universe. Galaxy clusters at different redshifts trace the evolution of regions with enhanced dark matter density, which provides valuable information on the cosmology-dependent matter density fluctuations. Studying galaxy clusters and protoclusters at the highest redshifts possible puts a meaningful constraint on the cosmological models of the universe.
While galaxy groups and clusters at low to intermediate redshifts (z≈1) have been thoroughly studied (e.g., Ellis et al. 1997;Stanford et al. 1998;Carlberg et al. 2001;Blakeslee et al. 2003;Eke et al. 2004;Halliday et al. 2004;Holden et al. 2005), only recently have the frontier facilities and instruments enabled us to venture into the realm of the high-redshift (highz) universe. For instance, a Virgo-like cluster with M vir ≈10 14.4 M e at z∼2.1 with multiple M halo ∼10 13 M e groups was discovered in the COSMOS field using a mediumband filter targeting the Balmer-break feature of the member galaxies (Spitler et al. 2012) and later spectroscopically confirmed (Yuan et al. 2014). Moreover, the mass and formation history of galaxy clusters are the objects of interest in many studies, such as testing cosmological models, constraining structure formation scenarios, deducing physical properties of dark matter, and constraining the nature of dark energy (e.g., Zwicky 1939; Sunyaev & Zeldovich 1972;Press & Schechter 1974;White & Rees 1978;Efstathiou & Eastwood 1981;Davis et al. 1985). For example, the massive high-z galaxy cluster observed by Gonzalez et al. (2015) with M 200 =(1.1±0.2)×10 15 M e is supposedly one of the five most massive galaxies at z1.19 (according to ΛCDM cosmology).
Many early studies, numerical simulations, and analytical calculations were devoted to studying the growth and formation of galaxy clusters: first, they originate from small gravitational instabilities and rapidly grow in both mass and size in the expanding universe; then, finally, they collapse and virialize (e.g., van Albada 1960van Albada , 1961Aarseth 1963;Peebles 1970;Icke 1973). On the observational side, the overdensities of galaxies have been detected at very high redshifts, even though they are not massive or mature enough to be classified as galaxy clusters. Many studies have discovered and identified the high-z overdensity of galaxies in the process of their formation as protoclusters (e.g., Steidel et al. 1998;Shimasaku et al. 2003;Ouchi et al. 2005;Wang et al. 2005). For example, Cai et al. (2017) found a rare overdense Lyα absorber at z=2.32±0.02 along with 20 Lyα emitters (LAEs), BOSS1441, in SDSS-III, by tracing the absorption features from approximately 80,000 quasi-stellar object (QSO) sight lines.
However, at very high z, many protoclusters have also been discovered in large galaxy surveys with various detection methodologies (e.g., Steidel et al. 1998;Shimasaku et al. 2003;Ouchi et al. 2005Ouchi et al. , 2018Wang et al. 2005;Overzier et al. 2006a;Toshikawa et al. 2012Toshikawa et al. , 2014Toshikawa et al. , 2016Lee et al. 2014;Higuchi et al. 2018). For instance, Toshikawa et al. (2012) discovered a massive protocluster of LAEs at z∼6 in Subaru Deep Field (SDF) using the i'-dropout technique (see also Toshikawa et al. 2014). One way to search for a protocluster at such high redshift besides conducting an extensive galaxy survey is to follow up on signposts of an overdensity (Overzier 2016). These are massive and biased toward dense environments, such as radio galaxies, QSOs, submillimeter galaxies, and starburst galaxies (e.g., Le Fevre et al. 1996;Carilli et al. 1997;Pentericci et al. 1997;De Breuck et al. 2002Barr et al. 2004;Venemans et al. 2004Venemans et al. , 2007Overzier et al. 2006b;Utsumi et al. 2010;Capak et al. 2011;Cucciati et al. 2014). Specifically, Davies et al. (2014) have proven that star-forming galaxies can be efficient signposts for high-z overdensities by discovering the >20σ clustering of 36 star-forming galaxies, a QSO, and a submillimeter galaxy (SMG) at z∼1.8 around LESS J033336.8-274401 in the COSMOS field. Capak et al. (2011) found a massive protocluster in the COSMOS field at the redshift of 5.3 using bright quasars and starburst galaxies (especially, Lyman-break galaxies, LBGs) as signposts. Jiang et al. (2018) have shown that using luminous LAEs as bias tracers is effective in searching for a massive protocluster at z=5.7. The recent success in searching for an overdense region at very high redshift using bias tracers has led us to adopt this method to look for a protocluster at z=6.5.
To look for the first site of large-scale structure formation and galaxy cluster assembly, we have pierced through the epoch of reionization by using the high-z active star-forming galaxies, called LAEs, as signposts of the overdensity. Even though they are fainter than quasars, LAEs prove to be more suitable for probing the faint end of the high-z star-forming galaxy luminosity function (e.g., Kashikawa et al. 2011;Konno et al. 2018). The visibilities of LAEs can also be slightly enhanced when they are in groups or clusters, due to the ionized bubbles in the intergalactic medium (IGM; Miralda-Escudé 1998; Dayal et al. 2009;Dayal & Ferrara 2011;Mortlock et al. 2011;Hutter et al. 2015). Chanchaiworawit et al. (2017) found 45 fainter LAE candidates clustered around the two spectroscopically confirmed massive LAEs at z=6.5 in the Subaru/XMM-Newton Deep Survey field (SXDS; Furusawa et al. 2008;Ouchi et al. 2010) by conducting OSIRIS medium-band photometric observations at the Gran Telescopio Canarias (GTC). We have also confirmed 10 LAEs via GTC/OSIRIS spectroscopic follow-ups (Calvi et al. 2018). In this paper, we discuss the clustering properties, the mass of the protocluster, and the predicted cluster mass at the present day inferred from the statistical analysis of photometric and spectroscopic results as pieces of evidence supporting the discovery of this protocluster. The calculations throughout this work adopt a ΛCDM cosmology with Ω M =0.3, Ω Λ =0.7, and h=0.7, consistent with the latest Planck cosmological parameters (Planck Collaboration et al. 2016). The magnitudes used in this work and the photometric data paper (Chanchaiworawit et al. 2017) are all in the AB system (Oke & Gunn 1983).

Observations
The observations of this overdense region of LAEs at z=6.5 were carried out in two phases: photometric selections and spectroscopic confirmations. The photometric selections were obtained from the observations on OSIRIS in its imaging mode at the 10.4 m GTC, using the three reddest medium-band filters from the SHARDS program (Pérez-González et al. 2013), namely F883w35, F913w25, and F941w33 (henceforth blue, green, and red bands). The sensitivities reached by the three bands were 26.54, 26.56, and 25.84 magnitude, respectively (Chanchaiworawit et al. 2017 Chanchaiworawit et al. (2017). There were 45 fainter LAE candidates selected from the photometric observations in addition to the two spectroscopically confirmed LAEs in the field (Ouchi et al. 2008(Ouchi et al. , 2010. The candidates were grouped into three classes based on the reliability of their photometry. Class I, II, and III LAE candidates were those that exhibit high signal-to-noise ratio (S/N?5) and compact circular shape (Φ∼1″), lower S/N∼5 or somewhat resemble a compact circular shape, and very noisy and spurious-like detections, respectively. Only class I and II were considered for our final set of LAE candidates (Chanchaiworawit et al. 2017). We only selected the LAE candidates that were in the area of the field of view (FOV) with high completeness level and low percentage of spurious sources (e.g., noise spikes, cosmic-ray residuals, fringe patterns of sky background), or the high-completeness and low-contamination regions as discussed in detail by Chanchaiworawit et al. (2017).
In semesters 2016B and 2017B, we conducted a spectroscopic follow-up of 16 LAE candidates and one spectroscopically confirmed LAE using the multiobject spectroscopy (MOS) capability of OSIRIS at the GTC. The spectroscopic observations spanned 36 hr and reached the sensitivity in flux density of 5×10 −18 erg s −1 cm −2 Å −1 . We defined four criteria to gauge the success of the spectroscopic confirmation of each observed LAE: (1) the detected emission line shows asymmetry resembling the P-Cygni profile in the 1D spectrum, (2) the observed size of the Lyα emission feature is comparable to the seeing during a typical spectroscopic night (1″), (3) the emission line is located far from prominent sky emission lines or their wings, and (4) the spectroscopically measured Lyα luminosity must be consistent within the uncertainty to the photometrically estimated Lyα luminosity during the photometric selection phase. The spectroscopically detected LAEs were graded using these criteria. Those that satisfy at least three or more criteria were given grades "B" and "A," respectively. We extracted the physical properties, such as redshift and Lyα flux, from the observed emission line by fitting a skew-Gaussian profile (O'Hagan & Leonard 1976;Mudholkar & Hutson 2000;Azzalini 2013) to the 1D spectrum. The skew-Gaussian profile was selected to mimic the asymmetry of the Lyα line at high redshift. Further details on the spectroscopic observation of the LAE candidates at z=6.5 and its results have been discussed in Calvi et al. (2018). We have obtained the following crucial information from the spectroscopic observations: (1) Ten LAEs show the profiles of Lyα emission that satisfied the spectroscopic confirmation criteria as presented in Calvi et al. (2018) from the total of 16 new LAE candidates observed. This success rate in spectroscopic confirmation of 10 16 is in excellent agreement with the previously estimated spectroscopic success rate of 2 3 as adopted by Chanchaiworawit et al. (2017). (2) The extracted 1D spectra of all of the observed LAEs and LAE candidates provide accurate redshift measurements, required for pinpointing the position of the object along the line of sight (LOS). (3) The fitted spectra of Lyα emission lines constrain the Lyα luminosities and the star formation rates (SFRs) of the LAEs.
However, when it comes to determining the Lyα luminosity function, we decided to use the photometric estimations of Lyα luminosities for the following reasons: (1) the S/N of the emission line in spectroscopic data was limited by contaminations from night-sky OH glows, and (2) there were only 10 LAEs with reliable spectroscopic measurements of the Lyα luminosities (Chanchaiworawit et al. 2017;Calvi et al. 2018).
The analysis of the luminosity functions of Lyα emission derived from this particular region as compared to the one derived from the 1 deg 2 SXDS and SDF (Furusawa et al. 2008;Ouchi et al. 2008Ouchi et al. , 2010 at z=6.5 showed at least a factor of 2× overdensity (Chanchaiworawit et al. 2017). The result suggests the existence of a protocluster in this region. Furthermore, the spectroscopic redshift measurements of the Lyα lines provided us with some insight into the underlying redshift distributions of these LAEs. With some handle on the spatial distributions of LAEs in 3D space, we can now reassess the survey area and effective volume, gaining a better constraint on the level of overdensity of LAEs at z=6.5.

Survey Area
The previously determined survey depth was derived from the FWHM of the F913 filter function. The depth was equivalent to the LOS of z=6.4 to z=6.6, while the surface area was derived from the full FOV area (Chanchaiworawit et al. 2017). This approach might overestimate the effective survey volume for the following reasons. First, the LAEs may not be uniformly distributed along the LOS from z=6.4 to 6.6. Second, the LAE candidates were only detected in highcompleteness and low-contamination regions of the FOV, as shown in Figure 1. Based on the results in Chanchaiworawit et al. (2017), the medium-band imaging mode of OSIRIS, which operates off axis, causes the wavelength variation of the filter's central wavelength across the FOV (see also Méndez-Abreu et al. 2011;Pérez-González et al. 2013). This variation introduces differential sky emission and, thus, creates the rising gradient of background noise from left to right as seen in the right panel of Figure 1. In the same figure, the overplotted distribution of LAE candidates also confirms the conclusion that only approximately one-half of the FOV is deep enough to detect such faint high-redshift LAEs.
The revised survey area was derived from the total area with high visibility (high completeness and low contamination) of LAEs. This area corresponds to the probability of the real detection in the F913 band of 0.5 (P(LAE)0.50) in the faintest available magnitude bin (i.e., 26.2<m F913 26.6), as demonstrated in Chanchaiworawit et al. (2017). The revised survey area as shown in Figure 1 is 48% of the total FOV area and contains all of the LAEs and LAE candidates regardless of their F913 magnitudes.

Survey Depth
In the effort to determine the appropriate LOS depth of the survey, we simulated the visibility of Lyα emission lines at different redshifts in the range of z=6.3-6.7. Under the assumption of the photon-noise-limited regime with sky emission features as the primary source of noise, we simulated photometric S/N as a function of redshift. The black dashed line and the light brown shaded region represent the S/N of the average LAE in photometric selection as a function of redshifts and 1σ uncertainty. The result was highly dependent on the shape of the F913 filter function, as demonstrated in Figure 2. The coadded LAE from Calvi et al. (2018) was the basis for the model of an average Lyα emission line at z≈6.5. This model LAE was built from averaging the flux profiles of the 10 new spectroscopically confirmed LAEs (excluding the brightest one). It corresponds to a Lyα luminosity of ∼2.5×10 42 erg s −1 and F913 magnitude of 26.4. Thus, this model is truly the representative of the faintest magnitude bin in our photometric survey. The S/N in the photometric selection of an average LAE as a function of redshift shows the enhanced probability of detection around z=6.42-6.49. However, all of the confirmed LAEs in this field are located beyond z=6.49 (Calvi et al. 2018). This finding rules out the possibility of a uniform distribution of LAEs along the LOS. The contribution of the LAEs in the z=6.41-6.49 range should be less than 6 16 (approximately one-third) or the ratio of unsuccessful confirmation to the total number of LAE candidates observed spectroscopically. When taking into account the redshift measurements of the spectroscopically confirmed LAEs (Calvi et al. 2018), we concluded that at least 2 3 of the LAEs in this field reside within the range between z=6.49 and 6.62. However, we need a better description of the underlying redshift distribution of the LAEs in this field.

Redshift Distribution
The spectroscopic redshift measurements of the LAEs and LAE candidates from the spectroscopic follow-up were not a complete representation of the underlying redshift distribution. However, these redshift measurements were drawn from an incomplete but uniformly distributed sample of LAEs on the FOV. Thus, we can use this information to constrain the underlying redshift distribution. We constructed a redshift histogram of the total 14 LAEs in our field. The sample comprised two previously confirmed LAEs (Ouchi et al. 2010), 10 new individually confirmed LAEs from our observations, and two spectroscopically confirmed LAEs from Higuchi et al. (2018). Among our 10 spectroscopically confirmed LAEs, there are three LAEs in common with the spectroscopically confirmed sample by Higuchi et al. (2018). We performed the Kolmogorov-Smirnov test (KS test) with various models for redshift distribution, as shown in Figure 3. The tested distribution models consisted mainly of Gaussian (normal) and uniform distribution functions. The boundaries of the uniform distribution function are at z=6.40 and z=6.62 (model 7). Model 1 was derived from a Gaussian distribution with the median and standard deviation (σ) equal to the weighted average and the standard deviation around the mean of the visibility function above the threshold of S/N4.0 as in Figure 2. Similarly, model 2 was derived from a Gaussian distribution with the median and standard deviation equal to the median and standard deviation around the median of the visibility function above the threshold of S/N4.0. Models 3, 4, and 5 were derived from Gaussian distributions with the median equal to the median redshift of the spectroscopic redshifts of LAEs and LAE candidates in our field. The standard deviations of models 3, 4, and 5 are the standard deviation of the spectroscopic redshifts of our LAEs and LAE candidates, and the standard deviations around the mean and median of the visibility function above the threshold of S/ N4.0, respectively. Model 6 was generally similar to model 1 and model 4, but the median was shifted to the position of average redshift of the large-scale (FOV∼0.1 deg 2 ) overdensity around our field, as reported by Higuchi et al. (2018). All of the fitted models and the redshift histogram of the LAEs are listed in Table 1 and shown in Figure 3.  . Redshift histogram of all 14 spectroscopically observed LAEs (two confirmed massive LAEs, 10 recently confirmed LAEs from this work, and two spectroscopically confirmed LAEs from Higuchi et al. 2018). The KS test was performed to find the best distribution function for the observed distribution along the LOS of the LAEs. The KS test statistics and parameters are listed in Table 1.  Note.
(1) Models were used in the KS test; (2) type of distribution functions of the models; (3) location of median redshifts of the models (the lower bound for the uniform distribution); (4) standard deviations of normal distributions (width in case of uniform distribution); (5) KS test statistics; (6) p values derived from the test statistics, representing the probability of the fitted functions being drawn from the same original distribution.
We used the one-sample KS test to find the probabilities that the spectroscopic redshift distribution of our LAEs and LAE candidates was drawn from the listed distribution models. We found that we could not reject the null hypothesis for models 3, 4, 5, and 6 with 90% confidence level, as shown in Table 1. Model 5 showed the highest p value. However, we could not ignore the other models that satisfied the null hypothesis, due to the small sample size of the redshift measurements. Furthermore, to avoid being biased toward our spectroscopic sample, we chose model 6 (derived from the suggested large-scale structure in Higuchi et al. 2018) to be the representative underlying redshift distribution of LAEs in this field. Model 6 also satisfied the requirement that at least ∼67% of the LAEs in this field are located in the range z=6.49-6.62, showing ∼85% of the area under the curve within this redshift range.

Effective Survey Volume
The revised effective survey volume corresponds to the derived LAEs' redshift distribution and their visibility across the FOV as discussed in Sections 2.2 and 2.4. The resultant comoving volume is 9242±1427 cMpc 3 , just slightly below one-third of the previously determined volume from the width of the green filter and full FOV of the survey. This volume contains 85% of the LAEs according to the best-fit redshift distribution discussed in Section 2.4. However, there is a caveat in this approach, since it may disregard the void surrounding the overdensity (e.g., Hayashino et al. 2004), or some faint LAEs farther away from the proximity of the overdensity are less enhanced in their visibility. Therefore, with this revised survey volume and the number of LAE candidates bounded by this volume (which is only slightly fewer than the total number of LAE candidates), the level of overdensity arising from this volume becomes at least 3× that of the surrounding area of the same epoch. This overdensity was derived from the comparison between our LAEs in this survey volume and the LAEs found by Ouchi et al. (2010) over their entire volume surveyed in the 1 deg 2 SXDS field. Note that the results of the entire SXDS field revealed a mild overdensity but were still consistent within 3σ of the cosmic variant, as concluded by Ouchi et al. (2010). Moreover, our observations focus on the subfield with the densest Lyα luminosity function from among five SXDS subfields (Ouchi et al. 2010). Thus, this result strengthens the evidence for the existence of an overdense structure, which could potentially be a protocluster in this field.

Density Contrasts
With the revised effective survey volume, the observed number density of LAEs becomes 1.69 0.77 We express the level of overdensity of LAEs in the field in terms of a density contrast, δ, as shown in Equation (1) (see Carroll et al. 1992;Steidel et al. 1998;Weinberg et al. 2013): where the subscript "i" represents the type of density used in the calculation, such as δ gal for galaxy number density contrast, δ m for matter density contrast, and δ l for linear perturbation density contrast. We compute the number density contrast δ gal with our revised observed number density of LAEs and the average observed number density of LAEs found in the extended SXDS and SDF surveys (Furusawa et al. 2008;Ouchi et al. 2008Ouchi et al. , 2010. We obtain δ gal =3.18 1.99 3.47 -+ for the LAEs in this field. Conveniently, matter density contrast is linked to the observed galaxy number density contrast as expressed in Equation (2) where C is the redshift space distortion factor of the collapse structure and b is the bias parameter (Kaiser 1987;Steidel et al. 1998;Overzier 2016). The number density contrast for each type of galaxy has different associated bias levels depending on the probability of finding that type of galaxy in a dense dark matter halo. Ouchi et al. (2010) give the bias parameter for LAEs at z∼6.6 to be 3.7. However, δ m and C are interconnected and affected by one another, as also shown in Equation (3): Here, f (z) is a function of redshift and depends on the cosmological model (Linder 2005;Overzier 2016 where Ω M (z) is the matter density parameter as a function of redshift, z. Thus, the value of C decreases as δ m grows, resulting in the survey volume appearing to be more compact.
But, as C decreases, δ m slightly decreases, as shown in Equation (2) (Steidel et al. 1998;Overzier 2016). Therefore, we simultaneously solve for both parameters by finding the intersection between Equations (2) and (3)  The redshift space distortion factor, C, affects the true volume of space that we sample in our observations as well. To get a handle on the mass of the overdense structure in this field, we compute the true survey volume by taking into account the effect of space distortion, expressed in Equation (6): where V survey and V true are the effective survey volume as observed with this overdensity and the distortion-corrected survey volume, respectively. Then, the mass of the overdensity is the comoving matter density of the field ( 1 m r d á ñ + ( )) times the distortion-corrected comoving survey volume, V true , as expressed in Equation (7): The derived mass of the overdensity, M survey , is 8.40 1.39 2.98 -+ ×10 14 M e . This mass is the estimated virial mass of the protocluster at the time of its collapse. The estimated mass of the overdensity is also in excellent agreement with the mass of the overdensity from other studies at similar redshift range. For example, Toshikawa et al. (2012Toshikawa et al. ( , 2014 give the range of the mass of the confirmed protocluster at z∼6 to be 5-10×10 14 M e . However, the estimation here is close to a lower limit for the present-day cluster mass for the following reasons: (1) only 85% of the expected number of LAEs reside within the volume and are used in the calculations; (2) the true survey volume is more compact than initially derived in the photometric selection phase and comparable to the typical volume of cMpc h 15 3 3 ( ) used in protocluster searches (Chiang et al. 2013;Higuchi et al. 2018); and (3) this estimation does not account for the possible major mergers with nearby groups or overdensity surrounding the protocluster, as presented by Higuchi et al. (2018) that there could be up to two other overdense structures within the distance of ∼7.5 cMpc from this protocluster (but outside the FOV of our survey).

Halo Masses of the LAEs
Another intuitive method to estimate the mass of the protocluster is summing up the halo mass of each LAE. This way of mass estimation would face a higher uncertainty from the completeness and contamination corrections. On the other hand, the total mass of the protocluster derived from summing up the halo masses of LAEs would miss the portions of matter associated with other types of galaxies (e.g., dwarfs, dusty starbursts, dust-obscured star-forming galaxies) or dark matter without prominent stellar sources, and diffuse IGM. Another source of uncertainty lies in the estimation of the individual halo mass itself. Nevertheless, the halo mass can only be inferred from the intrinsic Lyα luminosity. We need the value of the Lyα escape fraction, f esc,Lyα , at the redshift of interest to get a handle on the intrinsic Lyα luminosity, as demonstrated in Equation (8) et al. 2018), and f esc,Lyα is the escape fraction of Lyα, while the average optical depth can be expressed as a function of the LAE's halo mass as in Equation (9). Moreover, the best representative set of conditions that matched the observables of LAEs at z>6 was selected by Inoue et al. (2018)  Note that in some cases the normal distribution function of τ α (Equation (10)) yields a negative value, which is unphysical. When this occurs, our matching routine redraws a new τ α value until a nonnegative value is picked. Since there are approximately one-third of LAEs and LAE candidates in our field with Lyα luminosities obtained spectroscopically, we use the photometrically estimated Lyα luminosities and their uncertainties in the calculations of halo masses and escape fractions. The derived and best-match halo mass and Lyα escape fraction of each object are listed in Table 2.
The total halo mass of all LAEs and LAE candidates is shown in Equation (14) (Furusawa et al. 2008;Ouchi et al. 2010;Higuchi et al. 2018), and the suggested redshift evolution of the Lyα escape fraction by Hayes et al. (2011). The results from the SXDS field put the Lyα escape fraction for z∼6.6 at 0.30±0.18 (Ouchi et al. 2010), which is well within 1σ uncertainty from our derived average Lyα escape fraction of LAEs at z=6.5 as shown in Equation (15). However, there is one caveat to keep in mind: the initial distribution of the optical depth of Lyα photons, τ α,10 , has to be assumed based on the simulations by Inoue et al. (2018) in order for us to successfully solve for the hosting halo mass and Lyα escape fraction with the limited data we have (i.e., only the observed Lyα luminosity).

A Core-like Structure
We investigated the 3D structure of this overdense region of LAEs to find any further evidence of a protocluster. The information about the locations of LAEs and LAE candidates along the LOS is limited. Nevertheless, we have derived a probable underlying redshift distribution of the LAEs in this region, as discussed in Section 2.4. We apply this redshift distribution to fill in the blank and create a 3D grid of probability distribution for the LAE candidates without the spectroscopic redshifts.
The three dimensions of this information cube are R.A., decl., and LOS. The grid size is designed to be 55 pkpc in all directions, which should comfortably house a halo of a high-z LAE. Next, we create the cube for LAE halo mass, by multiplying the probability density as a function of xyz location with the probability of real detection of a particular LAE candidate, taking into account the spectroscopic success rate and each LAE halo mass as expressed in Equations (16) and (17): where S(x, y, mag) is the spurious fraction as a function of position on the FOV and F913 magnitude of an LAE candidate, P′(LAE) is the probability of the particular LAE candidate being a real detection, p(x, y, z) is the probability density of an LAE or LAE candidate belonging to a certain grid as a function of XYZ coordinates, M h í is an individual hosting the halo mass of an LAE, and Cube m (x, y, z) is the resulting halo mass falling within that grid. Then, the probability of detection being real becomes P′(LAE)=1 for confirmed LAEs. Furthermore, the probability density of finding an LAE at a specific location p(x, y, z) comes in two different forms. The probability densities of finding an LAE in the grid with the coordinates (x, y, z) for the spectroscopically confirmed LAEs and the photometrically selected LAE candidates are respectively expressed as normal distributions with the forms shown in Equations (18) and (19): The median and variance of the normal distribution are the previously determined median redshift of this overdense region and its variance as discussed in Section 2.4. The probability density goes to zero if the xy location does not match the location of the LAE or LAE candidate on the plane of the sky.
We sum up all of the probability-modified LAE halo mass grids that fall within the characteristic linking radius from the reference grid point. The characteristic linking radius at z∼6.5 is 0.45 pMpc (Ouchi et al. 2010). Next, we divide the sum of halo masses by the spherical volume with a radius of 0.45 pMpc. The calculated halo mass density is assigned to the reference grid point. The process is repeated for all grid points. The results of these simulations are shown in Figure 4, with each slice of Δz=0.01 thick and showing the contours of the average density along the LOS.
From Figure 4, the 6.50<z<6.51 slice shows the densest region of LAEs. Furthermore, to the left and right sides of this slice, similar patterns appear in lower intensities. This could be a hint of the front and back infalling of matter surrounding the center of a protocluster. However, this pattern is somewhat enhanced by the artifacts of a spherical volume in the calculations of halo mass densities. Nevertheless, the thickness of each slice is considerably larger than the characteristic linking radius of 0.45 pMpc. The total LAE halo mass within this core is M 1.41 10 0.35 0.29 12 -+  . Therefore, this core-like structure, revealed by our simulations of LAE halo mass density, provides another piece of evidence for a protocluster at z=6.5.

Is It Truly a Coma-analog Protocluster?
This dense core at the center of the observed overdensity harbors the two bright LAEs observed by Ouchi et al. (2010; hereafter, Ouchi's pair). These two galaxies are separated by only 305 pkpc and have spectroscopic redshifts of 6.511 and 6.519. Their relative velocity along the LOS is 2400 km s −1 . These two galaxies could be part of a virialized (or currently virializing) cluster around z=6.5 with σ LOS ∼400 km s −1 .

Harboring a Relaxed Core?
We pursue this line of argument for a relaxed core as the main structure of the protocluster by looking from various angles. From the first angle, we present the possibility of a pair of LAEs being members of the virialized core. On the other hand, we attempt to trace the radial density profile of the protocluster to see whether there is any LAE falling within the boundary of the caustic profile under the assumption of a virialized structure.

Friends-of-friends Analysis
To assess the possibility of a virialized core in this overdense region, we rely on the friends-of-friends (FoF) algorithm to find any grouping or clustering of confirmed LAEs in our sample. The FOF algorithm is a widely used procedure to look for virialized groups or clusters of galaxies in an observed or simulated sample (e.g., Huchra & Geller 1982). According to this technique, two galaxies belong to the same structure if their projected proper separation, Δr, is less than the "linking length" (two times the assumed virialized radius of the group), while their relative LOS velocity, Δv los , is less than the "linking velocity" (six times the assumed LOS velocity dispersion of the group).
A virialized halo should satisfy Equation (20): where f is a structural factor, which equals 0.7 for the Navarro-Frenk-White (NFW) dark matter density profile (Navarro et al.  Konno et al. (2018), and the fitted relation in Equation (12). The spherical subvolume with radius 450 pkpc is used to calculate the halo mass density. 1996,1997), and G is the gravitational constant. The radius of the defined virialized halo, or R 200 , is the largest extent that yields a matter density 200 times the average matter density of that epoch ( z 1 3 r á ñ + ( ) ), as shown in Equation (21) where r á ñ is the comoving matter density of the universe. With an assumption of an isotropic velocity distribution in 3D space, the LOS component of the velocity dispersion, σ los , is 3 s .
With a given σ los , the virialized radius of that epoch, R 200 , can be calculated as expressed in Equation (22) We have carried out the FoF analysis for the most probable case of a virialized core. Using the epicenter between the two Ouchi LAEs as the center of the group, we find the minimum σ los that could bound the pair within R 200 =150 pkpc from the group's center to be 572 km s −1 . This gives M 200 of 4.87×10 13 M e from solving Equation (21). We find that this minimum virialized mass would contain not only the Ouchi pair (C1-01 and C1-02), but also another newly confirmed LAE, C1-15, as well. The derived M 200 is indeed less massive than the mass of the whole overdensity.
Nonetheless, the calculated mass is still quite massive for a virialized halo at z=6.5. The could be the direct effect of the assumption of virialization, which may not hold at such an early stage of galaxy assembly, and also the degeneracy between LOS depth and LOS velocity of an LAE derived from the spectroscopic redshift. The derivation of the radial density profile and the caustic analysis in Sections 4.1.2 and 4.1.3 could shed some light on this problem.

Density Profile of the Protocluster
The derivation in the previous section assumed that the differences in redshift of the LAEs are caused by their LOS component velocities from the infalling toward the core of the overdensity. However, these differences in redshift can also be treated as the different LOS depths as well. Even though the derived virial mass from the previous section is quite large for that epoch, we cannot discard the possibility that the overdensity harbors a virialized core. Thus, the differences in redshift could not be purely geometrical effects, but rather the combinations of both different LOS depths and orbital motions around the core.
In this section, we assume that the distribution in redshift space is purely geometric and derive the radial density profile of the overdensity. The derived radial density profile will be used to perform a caustic analysis as a means of determining whether or not the overdensity or a part of it is virialized. The only requirement to use this caustic analysis to determine the level of virialization of the structure is to obtain a reasonable constraint on the matter density profile. Even though we cannot obtain a precise distribution in 3D for all LAEs and LAE candidates, the derived probability distribution of redshifts as discussed in Section 3.4 can help us secure the matter density profile with little challenge.
First, we find the best-fit matter density profile. However, the derived LAE halo mass density maps in Figure 4 only trace the mass of dark matter halos directly associated with the LAEs and not others (non-LAE galaxies, diffuse dark matter without stellar components, IGM, and so on). We cannot use them as a basis for constructing the matter density profile.
We go around this issue by following the procedures in Section 3.4 to construct a data cube of LAE number density. Next, we obtain the probability-normalized number counts of LAEs bounded by the assumed spherical volume. The spherical volume used here is smaller to give a higher frequency in the sampling of the number density profile, with R sphere =0.30 pMpc. From this number density data cube, we construct a radial profile of LAE number density by averaging the values of number density grids with the same distance from the center of the overdensity. The center of the overdensity is assigned from the densest grid point in the data cube.
The radial profile of LAE number density, or 1+δ gal as a function of distance from the center, is plotted in Figure 5 along with the available overdensity profiles of protoclusters at z=5 in the "Millennium simulation" by Chiang et al. (2013). Our radial profile of LAE number density is distinguishably denser than the one from star-forming galaxies of the Coma-clusteranalog protocluster within the inner 2 cMpc. Our number density profile is consistent within the uncertainty of the ones from star-forming galaxies in all types of protoclusters at z=5 within R2 cMpc. The result suggests the inner 2 cMpc of our protocluster harbors a dense core-like structure, as demonstrated in Section 3.4. However, our data set is plagued with large uncertainty and limited information on the 3D distributions of the LAEs. Therefore, we cannot draw any decisive conclusions from this number density profile alone.
Then, we transform the radial number density profile to the matter overdensity profile using the previously determined bias parameter, b, and redshift space distortion factor, C, as demonstrated in Section 3.1. The physical matter density profile corresponding to the observed overdensity as a function of physical radial distance from the center is illustrated in Figure 6. We find the best-fit matter density via χ 2 fitting of the Figure 5. Overdensity profile as a function of comoving radius from the center of the overdense structure. The light blue solid line with error bars represents the overdensity profile of this protocluster of LAEs at z=6.5. The red, green, and blue sets of lines represent the seed of Coma-like, Virgo-like, and Fornaxlike clusters at z=5 based on the Millennium simulation, respectively. The solid, dashed, and dotted lines trace SFR1 M e yr −1 , M10 10 M e , and M10 9 M e galaxies, respectively (Chiang et al. 2013).
Einasto profile (Einasto 1965;Navarro et al. 2010;Dutton & Macciò 2014) with the form described in Equation (23): where the parameter, r −2 , is the radius at which the logarithmic slope of the density equals −2, and α is the curvature parameter of the profile. We used the Einasto profile to fit the observed radial matter density profile of the overdensity to avoid the diverging cusp at r=0 from the NFW profile (Navarro et al. 1996(Navarro et al. , 1997. We also include the profiles that give sufficient goodness of fit within 68% confidence level as well. Thus, we have the ability to (1) alleviate the effect of the circular argument of obtaining the caustic profile from the 3D distribution of the LAEs and use it to determine the membership of the same set of LAEs and (2) provide a range of possible radial matter density profiles that could yield the observed distribution of the LAEs under the limited sample size that we possess. The best-fit Einasto parameters are shown in Equation (24) These fitted matter density profiles are used to derive the radial gravitational potential profiles as a function of distance from the center of the overdensity in pMpc. The gravitational potential profiles yield the radial profiles of velocity dispersion, which are discussed in Section 4.1.3.

Caustic Analysis
We follow the standard treatment for derivation of the LOS velocity dispersion as a function of tangential (sky projection) radius from the center of the virialized structure (e.g., Kent & Gunn 1982;Kent & Sargent 1983;Diaferio 1999;Geller et al. 1999;Łokas & Mamon 2003;Rines et al. 2013). We first express the gravitational potential of this overdense structure in terms of escape velocity, as shown in Equation (25): where Φ(r) is the radial gravitational potential, and V esc represents the escape velocity at a particular radius from the center of the structure. From the previously obtained ρ(r), we calculate the radial profile of gravitational potential from the Poisson equation for gravitational potential of an isolated sphere (Serra et al. 2011), as demonstrated in Equation (26): where ρ(r) is the matter density as a function of radius from the center of the structure, or the peak of the overdensity. The escape velocity is linked to its LOS component with a function of the anisotropy parameter, β (Binney & Tremaine 1987), as shown in Equation (27): where g(β(r)) is some function of the anisotropy parameter, β.
The value of g(β(r)), itself, can be expressed in a simple relation as discussed in Serra et al. (2011) and also in Equation (28): where β(r) is the anisotropy parameters as a function of radius from the center of the cluster or the overdense structure. The value of β takes a form based on the ratio between the orbital component and the total magnitude of the velocity dispersion at a certain radius from the center of the overdensity, as shown in Equation (29): where V θ (r) and V r (r) are the orbital component and the total magnitude of the velocity dispersion as a function of radius. The value of β differs based on the assumption. For instance, β=0 for an isotropic orbits (i.e., V V V r = = f q | | | | | |), β=1 for purely radial orbits, and b = -¥ for purely circular orbits (Binney & Tremaine 1987;Merritt 1987). To test the virialization level of such high-z overdensity, we adopt the β value that reflects the early stage of a virialized structure. In this situation, the velocity distribution should be close to isotropic with a hint of radial orbits (starting to infall to the center of the overdensity). We select the anisotropy parameter 0.18 b á ñ = (Serra et al. 2011). The approximation of constant anisotropic level has been tested to work well for the inner part of the galaxy clusters, where β(r) is fairly constant and has small positive values. With this b á ñ, we can derive the LOS velocity as a function of tangential radius by solving for V esc,| | in Equation (27).
The derived caustic profile of our protocluster as a function of projected radius in physical space is presented in Figure 7 for both escape velocities and LOS components of the escape velocities (3σ los ) as functions of the projected radius. The result demonstrates that there could be up to three confirmed LAEs with a high probability of being virialized (fit within the boundaries of the caustic profile) within the radius of 0.25 pMpc. This size of the possible virialized region is comparable to the size of the core-like structure we have demonstrated in Section 3.4, which strengthens the evidence for the protocluster.
Furthermore, the total caustic mass can be derived by integrating the product of the best-fit density profile and shell volume, as demonstrated in Equation (30): where ρ(r), in this case, represents the best-fit radial density profile. The obtained caustic mass is in good agreement with the total mass derived from the density contrast within 1σ uncertainty, as shown in Equation (31)

( )
Considering that this integrated mass is dominated by the main core-like structure with only three confirmed LAEs as its potential members (C1-01, C1-02, and C1-15), it is understandable that it is a bit less than the estimated mass of the overdensity. The zero-velocity radius (Kent & Gunn 1982), or the radius in which the influence of a group or cluster starts to dominate the cosmological flow, is expressed in Equation (32): where t z is the age of universe at the redshift, z. The derived zero-velocity radius is R 0.98 pMpc . This caustic mass represents the central structure of the protocluster gravitationally attracting other less massive clumps to form a more massive cluster in the future (as suggested in Figure 4) , which is in good agreement with Ouchi et al. (2010) as well.
However, the derived caustic mass does not represent the mass that may be fully virialized, but rather the total mass within the sphere of influence of this dense central structure of the overdensity. To obtain the virialized mass, or M 200 , we calculate the average matter density of this structure with different scales around the densest part of the overdensity. The largest comoving radius that yields the average matter density larger than or equal to 200× the present-day average matter density of the universe, m,0 r , is the virial radius, =-+  . The derived virial radius suggests that not all of the three LAEs thought to be virialized are full members of this virialized core. However, with the level of uncertainty and the results of the caustic analysis, we cannot reject the possibility that these three LAEs are on their way to being virialized either.

Growth of the Protocluster
At this point, we have obtained the different estimations of mass associated with this overdensity of its putative relaxed core from various methods. Now, we attempt to predict its evolution (i.e., the redshift of collapse of the overdensity and the growth of the relaxed core until z=0). This derivation allows us to assess the masses of the protocluster or cluster at different stages of formation and check the consistency of our derivations.

Redshift of Collapse of the Overdensity
First, we estimate the redshift of collapse of the overdensity using the linear theory of density perturbations. To do so, we transform the matter density contrast, δ m , into the linear density contrast, δ L . The relation between the two quantities is given by Bernardeau (1994): . The cosmological growth factor, g(z), will amplify this overdensity as it evolves through redshifts (Carroll et al. 1992;Weinberg et al. 2013). The parameter g(z) is expressed as a simple function of cosmological parameters based on the flat ΛCDM cosmology, as shown in Equation (34) where Ω M (z) and Ω Λ (z) are redshift-dependent matter and dark energy density parameters, respectively. Then the growth of linear density contrast is dictated by the normalized cosmological growth factor, D(z), as expressed in Equations (35) and (36): where z coll is the redshift of collapse of the whole overdensity, and δ c is the critical density contrast to trigger the collapse. We solve Equation (36) for z coll , where the redshift-dependent Figure 7. LOS velocity of the LAEs as measured from Δz (with respect to the central redshift of the overdense structure) in the inner 2 pMpc radius. Each LAE is represented by a green circle with ±1σ uncertainty of its V los . The caustic profiles are the orange-shaded regions representing the 68% confidence intervals of 3σ los derived from the fitted Einasto profiles, while the light blue hatched regions show 68% confidence intervals of V esc from the overdensity. The top right subpanel is the zoom-in of the inner 1 pMpc region.
The value of δ c is almost constant with very weak dependence on redshift from z=6.5 to z=0. In general, one can use the value of δ c =1.69. The resulting redshift of collapse of the detected overdensity at z=6.5 (e.g., Steidel et al. 1998;Weinberg et al. 2013 The detailed evolution of the protocluster can be followed using the extended Press-Schechter (EPS) formalism (Bond et al. 1991) starting from some initial relaxed seed. Following Manrique et al. (2003), one can integrate the instantaneous merger rate of a halo of mass M at a cosmic time t (Bond et al. 1991;Lacey & Cole 1993) over the mass of the accreted objects up to the threshold mass for major mergers (equal to one-third of the total halo mass). This yields the halo instantaneous accretion rate. Integrating over this accretion rate, we can obtain the expected mass increase over any period of time. The assumption that the halo grows by simple accretion (minor mergers) is justified, because the estimated virial mass of this overdense structure is already so massive that the probability of it merging with a similarly massive seed over the rest of its evolution (a major merger) is negligible.
Such an evolution can be followed, starting from the possible relaxed core at the center of the overdensity at z=6.5. In Figure 8, we show the result of growing the relaxed core at z=6.5. For comparison, we also plot the expected evolution of the most massive progenitor for current halos of the Coma cluster regime according to Chiang et al. (2013). At the derived redshift of collapse, the predicted virial mass according to the EPS formalism is ∼8×10 14 M e , which is in excellent agreement with the total mass of the overdensity as well. The result is a piece of evidence for the existence of the massive, virialized or at least in the process of being virialized, core of a protocluster (M 200 ∼4×10 13 M e ) at z=6.5, which will grow through simple, nonmajor merger accretion as predicted by the EPS model. Moreover, the predicted mass of the cluster at z=0 equals M 1.54 10 0.69 1.12 15 -+  , comparable to the mass of the Coma cluster. The result further strengthens our proposition that this overdensity at z=6.5 is a Coma-analog protocluster.
At very high redshift, z>2, the difference in masses of the seed or the most massive progenitor shows signs of ∼2σ tension between our derivation and the results of the Millennium II simulation, as shown in Figure 8 (Chiang et al. 2013). One possible explanation for this tension is that the assumption of virialization of a group of LAEs at the center of the overdensity may not hold at such high redshift, because there is a degeneracy in transforming differences in redshift to LOS depths and LOS velocity dispersion. Also, the virialized group of LAEs may be sparser than previously thought. However, from the caustic analysis in Section 4.1.3, we cannot discard the probability of these three LAEs being the members of a virialized structure.
Nevertheless, the 2σ consistency of the Millennium II simulation and our work suggests that, if it exists, this virialized group of LAEs at the center of the overdensity at z=6.5 is indeed very massive and quite rare. At the lower redshift, from z=2 down to z=0, our derivation of the cluster mass and the results of Millennium II simulations agree very well (within 1σ level). Thus, we are confident that this overdensity would produce a Coma-analog protocluster, which will evolve into one of the most massive clusters (M∼10 15 M e ) by z=0. Furthermore, with the level of number density contrast observed in this field, recent computational studies have put the mass of the resultant cluster to at least many times 10 14 M e (e.g., Toshikawa et al. 2014) and most likely at the level of the Coma cluster or ∼10 15 M e (e.g., Higuchi et al. 2018). Moreover, the wide FOV survey of LAEs at z≈6.6 by Higuchi et al. (2018) also shows two other dense regions seemingly interconnected to our field, with the projected separation of ∼1 pMpc. This finding could bring up the total mass of the protocluster and, ultimately, the predicted mass of the cluster at z=0. Therefore, with all of the pieces of evidence presented here, we are confident that this protocluster will become a massive Coma-like cluster at z=0 with mass on the order of 10 15 M e .

Conclusion
We have conducted both photometric and spectroscopic observations of the overdense region of LAEs at z=6.5. The photometric phase with OSIRIS in its imaging mode on the 10.4 m GTC has revealed 45 fainter LAEs or LAE candidates (Chanchaiworawit et al. 2017). The spectroscopic phase with MOS mode of OSIRIS at the GTC followed up on 17 LAEs and LAE candidates (one confirmed LAE from Ouchi et al. 2008, 2010, and16 LAE candidates). The spectroscopic follow-up has confirmed 10 fainter LAEs with sufficient S/N of Lyα emission lines, as demonstrated in Calvi et al. (2018).
First, the new spectroscopic results have suggested that the previously determined survey volume was overestimated, especially in the LOS dimension. The measured spectroscopic redshifts of the LAEs suggested a much tighter distribution along the LOS. We have revised the survey volume according to the effective survey area, or 48% of the full FOV, and the Figure 8. The most massive progenitor mass as a function of redshift of Comalike clusters is demonstrated as the black dashed line, with 68% and 95% confidence intervals displayed as opaque and transparent blue hatched regions, respectively (Chiang et al. 2013). The overplotted black solid line represents the evolution track of the possibly virialized seed with mass of 4.06×10 13 M e at z=6.5, with 68% and 95% confidence intervals displayed as opaque and transparent orange regions. derived underlying distribution along the redshift space of the LAEs, which was centered at z=6.537 with standard deviation σ z =0.048. With the revised survey volume, we have found the LAE number density contrast of δ gal = . The caustic analysis has also pointed out that there is a nonnegligible probability of a group of three LAEs at the center of this overdense structure being virialized. While the derived virial radius is R 200 =1.06±0.20 cMpc, the virial mass is M 4.06 10 1.90 2.78 13 -+  . This mass estimate is indeed quite massive for a protocluster core at such high redshift. Nevertheless, the virial mass is still consistent within the 2σ level from the Millennium II simulation results at z=6.5. We have found that the evolution of this seed can grow into a Coma-type cluster at the present-day universe with mass ) . This cluster mass at z=0 is consistent with that of a very massive cluster we have seen in the local universe. Therefore, we are confident that this observed overdense region of LAEs harbors a seed of a massive cluster of galaxies, in other words, being a Comaanalog protocluster at z=6.5.