More than 200 quasars at z > 5.6 were discovered in optical and near-infrared surveys, such as SDSS (e.g. Jiang et al. 2008;Fan et al. 2006), CFHQS (e.g. Willott et al. 2007;Willott et al. 2010), UKIDSS (e.g. Venemans et al. 2007;Mortlock et al. 2011), VIKING (e.g. Li et al. Venemans et al. 2013;Venemans et al. 2015), VST-ATLAS (Carnall et al. 2015), DES (Reed et al. 2015;Yang et al. 2019), HSC (Matsuoka et al. 2016), Pan-STARRS1 (e.g. PS1; Bañados et al. 2014;Bañados et al. 2016;Wang et al. 2019), etc (e.g. Yang et al. 2021). Galaxy formation simulations and cosmological models predict that high redshift quasars reside in the most massive dark matter halos and grow quickly in the galaxy overdensity environments formed by the hierarchical merging of galaxies (e.g. Romano-Díaz et al. 2011;Costa et al. 2014). Discovering the overdensities of galaxies associated with z ∼ 6 quasars will trace mass assembly and probe the evolution of protoclusters and large-scale structures in the early universe. Hence, it is a chance to investigate the environment where QSOs hosting SMBHs reside and to search for galaxy protoclusters close to the end of cosmic reionization.

Over the last twenty years, sub-mm surveys around AGNs and galaxies have advanced our understanding of the galaxy environment at the cosmic time. The first one was uncovered by SCUBA (Submillimeter Common-User Bolometre Array) surveys on JCMT (the James Clerk Maxwell Telescope, Smail et al. 1997). Subsequently, a series of papers have shown evidence for overdensities of submillimeter galaxies (SMGs) at z = 2 ∼ 3 in the vicinities of luminous radio galaxies (e.g., De Breuck et al. 2004;Greve et al. 2007) and dust absorbed QSOs (Stevens et al. 2004;Priddey et al. 2007). The clustering of SMGs at z ∼ 2 around luminous QSOs or galaxies indicates evolutionary scenarios in which powerful starbursts and QSOs occur in the same systems (Hickox et al. 2012). And Priddey et al. (2008) also reported deep sub-mm images with a large field of view for three z > 5 quasars in SCUBA observations. However, the observational evidence of large-scale structure close to z ∼ 6 quasars is elusive, limited by a relatively small quasar sample size and strong selection biases of quasars.

There are many different approaches for tracing overdense regions around quasars or galaxies at z ∼ 6, such as Lyman-Break Galaxies (LBGs) selected via broadband optical/NIR colors (Stiavelli et al. 2005;Morselli et al. 2014) or Lyman alpha emitters (LAEs) selected via narrow/broadband observations (e.g., Mazzucchelli et al. 2017;Ota et al. 2018). Despite the huge investment of optical or near-infrared (NIR) telescope time, the galaxy environments surrounding z ∼ 6 quasars remain unconstrained. Overdensities of galaxies were discovered with deep optical imaging around several quasars at 4 < z < 6 (e.g., the z = 6.28 quasar SDSS J1030+0524 in Stiavelli et al. 2005, and the z = 5.81 radio-loud quasar SDSS J0836+0054 in Zheng et al. 2006 andAjiki et al. 2006). However, some works failed to find any significant galaxy overdensity (e.g., LAEs surrounding PSO J215-16 at z = 5.73 in Mazzucchelli et al. 2017, and CFHQS J2329-0301 in Goto et al. 2017). Integral field spectroscopic investigations at radio/sub-mm wavelengths in quasar fields both with ALMA (Decarli et al. 2017;Willott et al. 2017;Venemans et al. 2019) and NOEMA (Wang et al. 2011;Omont et al. 2013) revealed the presence of close (< 50 kpc, < 1cMpc) companion galaxies in the immediate proximity of the first massive black holes. It only traces the small-scale clustering at close separation from the central quasars.

To test whether the companion galaxies identified around some z ∼ 6 quasars are a small-scale phenomenon or the expression of the long-sought largescale overdensities postulated by the models, in 2016-2018 we mapped the fields around 54 quasars at z ∼ 6 over a much larger field-of-view (∼15 in diameter) using SCUBA-2 (Submillimeter Common-User Bolometre Array-2) on JCMT, named as 'JCMT-SCUBA2 High rEdshift bRight quasaR surveY' ('SHERRY' survey). SCUBA-2 has a field of view of 64 sq-arcminutes (only 4.3 sq-arcmins for SCUBA) -a factor of 16 times improvement over SCUBA (Chapin et al. 2013). And it has the fastest mapping speed at 450 µm and 850 µm, which makes searching for SMGs more efficient (Chapin et al. 2013). The SHERRY survey is the first to build a large, statistical sample of z ∼ 6 quasars for sub-mm observation, including ∼ 20% of the z ∼ 6 quasars discovered to date.

In the first paper of this survey (SHERRY-I, Li et al. 2020a), we reported that about 30% (16/54) quasars were detected with a typical 850µm RMS sensitivity of 1.2 mJy beam -1 (> 3.5σ). The broadband SED for z ∼ 6 quasars is similar to that at z ∼ 2, which indicates there is probably no evolution of quasars' properties with redshift. The optical/near-infrared (NIR) spectra of these objects show 11% (6/54) of the sources have weak ultraviolet (UV) emission line features. The SHERRY survey also indicates that quasars with submm detections tend to have weaker emission lines than quasars with non-detections. The weak-line quasar at z ∼ 6 may be a young AGN-galaxy system (Li et al. 2020a), in which the broad line region is starting to develop slowly and have some unusual properties.

In this paper, we present JCMT/SCUBA-2 observations of the environment around the quasars at 5.6 < z < 7.1 and search for galaxy protoclusters close to the end of cosmic reionization. We describe the SCUBA2 observations and data reduction in Section 2. In Section 3, we introduce our method to extract SMGs candidates and calculate the number counts. We present the num-

The map depths at 850 m map / mJy ber counts results and the discussions in Section 4. Finally, we give a summary of the main results in Section 5. We adopt a Λ-model cosmology with H 0 = 71 km s -1 Mpc -1 , Ω M = 0.27 and Ω Λ = 0.73 (Spergel et al. 2007) throughout this paper.

We observed 54 quasars at 5.7 < z < 7.1, which have rest-frame 1450 Å absolute magnitudes of M 1450 ≤ -25.4. The observations were carried out with the SCUBA2 camera on the JCMT in Hawaii from August 2015 to January 2018 (Program ID: M15BI055, M16AP013, M17AP062, M17BP034). We used 'the constant velocity Daisy scan' observation mode ('CV DAISY' mode) with the field of view of 15 in diameter. This mode is for small and compact sources of order 3-arcmin or less. The most uniform and best rms region is the central one with a diameter of about 4 , corresponding to comoving 9.8 h -1 Mpc comoving. This covers typical protocluster scales at z ∼ 6 (e.g., Overzier et al. 2009).

SCUBA2 observes two bands, 450 µm and 850 µm simultaneously. The main beam sizes of SCUBA-2 are 7.9 at 450 µm and 14.9 at 850 µm; with the spatial resolution of 4 / pixel and 2 / pixel, respectively. The observations were carried out in grade 2 / grade 3 weather conditions with precipitable water vapor (PWV) of 0.83 -2.58 mm at zenith (atmospheric optical depth 0.05 < τ 225GHz < 0.12). We observed each target in four to six ∼ 30-minute scans with a total on-source time of 2-3 hours to reach the sensitivity of 1.2 mJy at 850µm . We present the details of the observations in Li et al. (2020a). The calibrated sources were observed before and after the target sources and selected from the James Clerk Maxwell Telescope (JCMT) calibrator list (Dempsey et al. 2013). Calibration sources included Mars, Uranus, Neptune, etc.

We reduced the data using the STARLINK SCUBA-2 science pipeline with the configuration file of 'dimmconfig blank field.lis', which is the recipe for processing maps containing the faint compact sources (Chapin et al. 2013). We processed each complete observation separately to produce an image and calibrated it with the flux conversion factors (FCFs) to mJy/beam. We adopted the default FCFs value of 537±24 Jy pW -1 beam -1 for 850 µm map and 491±67 Jy pW -1 beam -1 for 450 µm map (Dempsey et al. 2013). We combined all images for a given source into a single image using inverse-variance weighting. Using this recipe, we further processed the output map with a beam-match filtered of a 15 FWHM Gaussian at 850 µm band. The S/N is taken to enhance point sources, which is suitable for searching for SMGs at z ∼ 6. The final 850-and 450-µm maps have typical noise levels of 1.2 and 16.5 mJy beam -1 , shown in Figure 1.

DEC offest (arcmin)

Figure 3. JCMT-SCUBA2 850µm images of the surrounding regions around 54 z ∼ 6 quasar with a beam size of FWHM = 15 . The selected region is within 2.23 arcmins from the quasar optical position (which covers typical protocluster scales at z ∼ 6) and has one sigma uncertainty of < 1.5 mJy beam -1 for each pixel. The contour levels are +3, +4, +5, +6, +7×1.5 mJy beam -1 for each map. The star indicates the quasar optical/NIR position. The star with the blue edge is the sub-mm detected quasars. The blue circles indicate the selected SMGs (with S/N>3.5 at 850um or 450um, see Table 2). The sources marked with '#' are the SMGs used in the overdensity analysis with deboosted fluxes > 5.5 mJy.

-2

DEC offest (arcmin) The SHERRY survey fields are centered on the targeted z ∼ 6 quasars. The expected protocluster size is about ∼ 20 cMpc along the side, and the core of the protocluster is 2 -3 cMpc (Overzier et al. 2009;Chiang et al. 2017). For each map, we use the (square root of the) exposure time as a proxy for noise, which is as a function of radius (Figure 2). Note the actual noise map is produced from the pipeline, including other noise, e.g., the instrument and atmospheric noise. Considering that the noise increases rapidly at a large radius of the Daisy maps, we searched for overdensity within the central ∼2 arcmins around the quasar. This region corresponds to a scale of ∼10 cMpc in diameter, which can trace a typical protocluster (Chiang et al. 2013). Furthermore, we mask the region where the RMS noise level is > 1.5 mJy beam -1 based on the error map to ensure the best/uniform sensitivity in the SCUBA-2 images. The error difference between the center and the edge is less than ∼ 15%. The source extraction is processed in this effective region, shown in Figure 3.

To extract the sources in the maps, we adopted the 'top-down' algorithm to find the source with the highest SNR in the two bands maps for each quasar field (e.g., Chen et al. 2013;Geach et al. 2017). SMGs at high redshifts are unresolved in SCUBA-2 observations; thus, the flux density of the source is the peak value of the pixels in the beam. We first identified the pixel's position with the highest SNR for each source. Then we subtracted a scaled point spread function (PSF) from the image at the source position. We generated the Gaussian PSF model with an FWHM of 14.9 arcsec at 850 µm map and 7.9 arcsec at 450 µm . During each iteration, we recorded the coordinates, flux densities, noise values, and SNRs of each > 3.0σ detection of candidate sources. To create a unique source list, we removed multiple detections if the secondary peaks were within one beam area of a higher SNR one. We iterated this process until there were no more sources with SNR of greater significance than the threshold in the image.

After identifying the > 3.0σ candidate sources in the two band maps, we cross-checked the sources between 850µm and 450µm maps to find the counterpart in the radius of 15 (the beam size at 850µm map). The final catalog of submm source candidates is generated as those satisfying > 3.5σ at least in one band. The flux density of its counterpart on another band is also listed in the catalog, if it has. The 850µm images of each of the 54 quasar fields are shown in Figure 3. The final SMGs selected as described above are marked as blue circles, and listed in Table 2 of Appendix A. We also carried out MCMC deboosted process for each individual source's flux for repeating 50000 times, and calculated the deboosted flux and its error.

In this section, we determine the pure galaxy number counts in the quasar fields. To do so, we first need to estimate the number of spurious sources contaminating the counts and estimate completeness levels, particularly at low SNR. These completeness and flux deboosting calculations are all processed for each field individually. In this section, we use SDSS J0100+2802 as an example to show how it works. Note the number counts result in this paper is only based on 850 µm maps due to its better sensitivity.

We examined the negative noise statistics to calculate the contamination rate of false detections. For example, we plot the S/N histograms of the pixels in the effective area of the signal maps at quasar nearby field for SDSS J0100+2802 (z = 6.3258±0.0010) in 850µm and 450µm bands, in Figure 5. The red line represents the Gaussian fitting of the noise statistics. The excess above this line on the positive side represents true detections assuming that the noise is symmetric. On the other hand, the negative tail is from the spurious sources due to the matched-filter PSF (Chapin et al. 2013). We need to estimate the distribution of negative to obtain the pure number counts.

We performed this check on each field individually. We first constructed an inverted flux map and then applied the source extraction procedure described in Section 3.1 to them. Nine detections at σ > 4 were found in the inverted maps at 850µm band. We estimated the false detection rate above 4σ significance to be 11.4% (9/79). We also estimated the false detection rate at 3 < σ < 4 to be 57.6% (57/99). Considering the effects of spurious sources, we calculate the number counts with the flux bin beginning from 4σ. Then the impact of the spurious source will significantly diminish, and the conclusions below are not affected.

To describe how we make some corrections to calculate the number counts, we still use the field of quasar SDSS J0100+2802 (z = 6.3258 ± 0.0010) as an example in the following. The completeness is the rate at which a source is expected to be detected in the image (Hatsukade et al. 2016). We used the Monte Carlo simulation approach to estimate completeness. The completeness test is only in the central best/uniform-sensitivity regions. We injected artificial sources of 0 -10 mJy/beam at random positions in the maps > 8.0 away from ≥ 3σ peaks to avoid the contribution from nearby sources. The flux bin step is 0.5 mJy for 850 µm map. Then we ran the source extraction on the maps to see if we could recover the injected sources. When the output source is above 3σ, the source is considered to be recovered. This procedure is repeated 5000 times for each flux bin. We fit the completeness as a function of f(S in )=[1 + erf((S in -A)/B)]/ 2 (Hatsukade et al. 2016), as shown in Figure 6. Here S in is the input flux density. The best-fitted parameters are A = 4.67 ± 0.01, B = 2.05 ± 0.02. Based on this Monte-Carlo experiment, we estimate a 50% completeness level for 4.7 mJy point sources and 80% for 5.9 mJy point sources in 850µm map.

With the Monte Carlo simulations, we also considered the effect that flux densities of low S/N sources are boosted (Murdoch et al. 1973;Hatsukade et al. 2016). More specifically, we calculated the ratio between the flux of the recovered sources and those of the injected artificial sources. We then estimate the mean ratio in 0.5 mJy flux bins. To illustrate this process, we plot the mean flux boosting value for the field of J0100 in Figure 6 (right). The flux boosting is more pronounced for low S/N sources. We fit the ratio between input flux density and output flux density as a function of f(S in ) =1+Aexp(-BS C in ) (Hatsukade et al. 2016). The best fitted parameters are A = 1.74 ± 0.27, B = 0.11 ± 0.04 and C = 1.95 ± 0.22. At S/N = 4 (input flux = 4.8 mJy), the estimated mean flux boosting factor is 1.19 at 850 µm map. This is consistent with the results of Monte Carlo simulations in previous SCUBA studies (Eales et al. 2000;Scott et al. 2002;Cowie et al. 2002;Casey et al. 2013;Chen et al. 2013), such as the boosting factor of 1.25 at 850µm map in Cowie et al. (2002). Potential contamination for the SMG number counts is the gravitational lensing effects of foreground objects (e.g., Chen et al. 2011). If some of the detected submillimeter sources are lensed background sources, they may increase the uncertainties of the source positions and also boost the measured fluxes. (Knudsen et al. 2008) 4. RESULTS AND DISCUSSION

To obtain an unbiased estimate of the number counts for each SMG candidate, we statistically correct the observed flux density using the Monte Carlo simulation for each field, see Figure 6. Then, the deboosted flux density is amended for completeness, as described in Section 3.2.2. Note we also removed the quasar itself in the analysis of number counts.

There are two number counts to describe the source abundance: the cumulative number counts and the differential number counts. To study the general environment of z ∼ 6 quasars, we built a combined area to calculate the differential number counts, which have independent errors for each count measurement, making it straightforward to do fitting. The combined area includes all 54 quasar nearby fields and totals ∼ 620 arcmin 2 , corresponding to ∼ 3700 cMpc 2 (∼ 76 pMpc 2 ) at z ∼ 6. It represents the overall environmental property of quasar fields at z ∼ 6. Another number count is the cumulative number count. It applies to the small sample size with a visual underlying counts shape. This method can be used in the analysis of individual quasar nearby fields.

The differential number counts of the combined area are calculated by dividing each > 3.5σ detection by the total combined area and adding up all measurements lying with each flux bin, then dividing by the flux width. Following previous studies, we fit the differential number counts with a Schechter function of the form (Casey et al. 2013;Geach et al. 2017):

We find the best fit parameters are N 0 ∼ 3000deg -2 , S 0 ∼ 1.2 mJy and α ∼ 1.5. We show the corrected differential and cumulative number counts for the effective area of the 54 quasars, together with the fit results in Figure 7. We list the corrected data points in Table 1. Note that our fitting is restricted to sources brighter than 5.5 mJy. The first set of uncertainties is Poisson errors (Gehrels 1986). And the second reflects the field variance that we used a bootstrapping approach. The individual field was selected randomly as many times as the size of a given sample, allowing for replacements, to create a new sample to calculate the counts. And each source is also corrected deboosting and completeness. This process was repeated 1000 times, and then we fitted a normal distribution to these 1000 realizations.

The standard deviation of this distribution was the field variance. These variances are of comparable magnitude to the Poisson errors.

We then compared our results with the previous blank field survey -SCUBA2 Cosmology Legacy Survey (S2CLS, ∼5 deg 2 ) at 850 µm. To a 4σ detection (∼ 5.5 mJy), the SMGs overdensity (δ SMG = ρ SMG / < ρ SMG > -1) is 0.68 (Gehrels 1986). Points in this work have been corrected for flux boosting and incompleteness (see figure 6). of 0.17 degree 2 , exceeding (∼ 1.68 ×) the source counts of blank sky area. And to a 5σ detection, the SMG overdensity is a factor of 0.40 +0.27 -0.23 (±0.24). If the galaxy excess is mainly due to overdensities at a narrow redshift range around z ∼ 6, the SMG overdensity near quasar redshifts must be higher.

To investigate which quasars reside in an overdense environment and search for protoclusters at z ∼ 6, we also calculated the completeness and boosting corrected SMG cumulative number counts N(>S) individually for each quasar field. We compared the results with the number counts in the blank field and gave a rough definition of 'overdensity'. The lower limit of the SMGs cumulative number counts to 5.5 mJy (∼ 4σ) is above the value in the blank field. With this criterion, 13 quasar fields show clear sub-mm overdensities with the overdensity of δ SM G ∼ 1.5 -5.4. We label their overdensity in Figure 3. Here the calculation of the overdensity is based on the sources with deboosted fluxes brighter than 5.5 mJy. There are 55 sources from the whole survey (marked with '#' in Figure 3). 24 of these lie in the 13 over-densities fields, with four of the over-densities comprise just one source above this flux limit.

We also combined the overdense fields around these 13 quasars and compared them with the blank field and those around other lower redshift galaxies. For comparison we show the cumulative number counts of SMGs, in fields centered on (i) luminous high-redshift AGN: z ∼ 5 optically selected quasars (Priddey et al. 2008); z ∼ 4 radio galaxies (Stevens et al. 2003); (ii) z ∼ 2 protoclusters (Zhang et al. 2022); (iii) z ∼ 1 clusters (Best 2002;Webb et al. 2005) and (iv) z ∼ 6 optically selected quasars (in this paper) in Figure 8. The black line shows the blank field number counts of SMGs from the SCUBA2 survey of the S2CLS field in Geach et al. (2017). The combined field of these 13 quasars exceeds the blank field counts with the overdensity to 5.5 mJy of δ SMG ∼ 2.46 +0.64 -0.55 (±0.25) in the regions of ∼ 150 arcmin 2 . It is consistent with the lower redshift cluster and the surrounding regions of quasars and radio galaxies. But to the bright end, the excess is not significant, e.g., to 7.5 mJy, the overdensity decreases to 1.00 +1.08 -0.74 (±0.73). Compared with the lower redshift galaxies/quasar fields or cluster fields, the combined overdense area shows a lower overdensity. That is because the data from the literature are case studies that intentionally target overdense regions (including significant selected bias). Still, our data points represent the average environmental properties around z ∼ 6 quasars, and overdensity is diluted.

Therefore, considering the abundant SMGs candidates in the faint-end, we select the fields where the lower limit of the SMGs cumulative number counts to 6.5 mJy (∼ 5σ) is above the value in the blank field. There are three (3/13) z ∼ 6 fields that show a significant excess over the blank field: P007+04 (z = 6.0), J0100+2801 (z = 6.3) and P308-21 (z = 6.2), shown in Figure 8. They are more likely similar to the lower redshift overdense fields/clusters with the overdensity of 2.1 -5.4 at 5.5 mJy. However, Champagne et al. 2018 analyzed 1.2 mm ALMA dust continuum maps for 35 bright quasars at z ∼ 6. They reported there is no evidence of 0.1 -1mJy millimeter continuum sources in the environments of P007+04 (z = 6.0) and P308-21 (z = 6.2). The overdensity is -0.07±0.56, which implies no detected overdensity of dusty galaxies around these quasars. This conclusion is based on a small ALMA field of view of about 25 , corresponding to an angular diameter distance of ∼ 1 cMpc at z ∼ 6. It could miss the overdensity on scales larger than one cMpc.

Some previous works used other tracers, like LAEs and LBGs, to probe the large-scale structures. But the question of whether the galaxy environment traced by LAEs or LBGs of z ∼ 6 quasar is overdense or not is still far from clear. Morselli et al. (2014) showed that the environment traced by LBGs around four z ∼ 6 quasars is overdense in a wide area of ∼ 575 arcmin 2 , ∼3100 cMpc 2 . Conversely, Willott et al. (2005) reported three SDSS quasars at z ∼ 6 that have no clear signs for LBG overdensity in an area of ∼ 30 arcmin 2 (∼170 cMpc 2 at z ∼ 6).

At z = 2 -6, some observations reported that SMG brightness appears to be correlated with the highdensity region of LAEs and LBGs (e.g. Pope et al. 2006;Tamura et al. 2009;Matsuda et al. 2011;Hodge et al. 2013;Umehata et al. 2017;Harikane et al. 2019).

Here we compare some individual sources that have previous environment studies traced by LAEs and LBGs on a similar scale. Becker et al. (2018) reported that J0148+0600 field has a highly significant deficit of LAEs. It is roughly one-quarter of background density with ten cMpc (similar size in our survey). But in our millimeter survey, the environment of J0148+0600 has an overdensity of ∼ 2.1 (∼ 3 × the background field).

Conversely, Mazzucchelli et al. (2017) reported no evidence for an overdensity of LAEs in an area of ∼ 37 arcmin 2 centered on the quasar PSO J215-16. We also find no evidence of SMG overdensity in the field of PSO J215-16 in our SHERRY survey. At a similar scale, the different findings may be ascribed to different tracers. These LAEs or LBGs trace the young, dust-poor starforming galaxies. Nevertheless, some observations and

Our SHERRY survey covers 54 z ∼ 6 quasars (20/54 detected at the sub-mm band, see Li et al. 2020b). We identified 170 SMGs (> 3.5σ) at 450 and 850 µm, in an area of ∼620 arcmin 2 . SHERRY represents the first systematic study to probe the environment of z ∼ 6 quasars in the sub-mm bands. For these SMGs lie z ∼ 6, their FIR luminosities are 2.2 -6.4 × 10 12 L , and their star formation rates are ∼ 400 -1200 M yr -1 (assuming a modified black body SEDs with a typical SMG dust temperature of 35 K, an emissivity index of β = 1.6 and a Salpeter IMF).

We found that 13/54 fields around z ∼ 6 quasars exhibit excess above SMG counts compared to a blank field. However, we caution that the other 41/54 quasars could also reside in the overdense fields, because our survey is not sensitive to the fainter SMGs with L < 2.0 × 10 12 L . Furthermore, even in an actual dense region of galaxies, the expected number of massive galaxies will not be very large based on the luminosity function.

On the other hand, our SHERRY survey does not provide an accurate redshift determination. Therefore, the SMG overdensity may be diluted over the large cosmological volume. Meyer et al. (2022) presented ALMA observations for three quasar-SMG fields (P231-20 at z = 6.59, P308-21 z = 6.24 and J0305-3150 z = 6.61) selected from our SCUBA2 survey to search for [Cii] emission. 12 out of 17 SCUBA2-selected sources are detected with ALMA in the continuum. ALMA detections have offset from the center of pointing due to the large beam of SCUBA2 imaging (13 at 850µm , e.g., of the order of ALMA field of view, Figure 1 in Meyer et al. 2022). They also search for the [C ii] emission lines in six of the continuum-detected sources: four ones in the field of P231-20 (> 5σ) and two ones around P308-21 (3-4σ). In the field of PJ231-20, they confirmed a foreground cluster at z ∼ 2.4. But unfortunately, none of the sources have a photo-z consistent with the background z ∼ 6 quasars. In the future, more high-resolution imaging and spectroscopic observations are required for these JCMT-selected protocluster candidates. If we confirm that a significant fraction of the SCUBA-detected galaxies are indeed at the redshift of the quasars (e.g., P007+04, J0100+2802, etc.), we will have striking evidence that some high-z quasars reside in highly overdense regions.

In this paper, we have presented a SCUBA-2 sub-mm survey 'SHERRY' of 54 z ∼ 6 QSOs to study their environment and to search for galaxy protoclusters at the early epoch. The large field of view of JCMT enabled us to identify 170 SMGs (> 3.5σ) at 450 and 850 µm . SMGs are selected within a distance of ∼2 arcmins away from the central quasars and have the RMS noise level of ∼ 1.2 mJy beam -1 in the 850 µm band. This region corresponds to 10 cMpc in diameter, which can trace the typical galaxy protocluster (Chiang et al. 2013). If the SMGs are near the quasar redshifts, their FIR luminosities are 2.2 -6.4 × 10 12 L and star formation rates are ∼ 400 -1200 M yr -1 . The SMGs discovered in SHERRY provide us with essential candidates for follow-up spectroscopic observations (e.g., [C II ]) to secure redshifts and study the first massive large-scale structures in the universe.

We also calculated the SMGs differential and cumulative number counts after correcting for completeness and flux boosting in a combined area of ∼ 620 arcmin 2 , which corresponds to ∼ 3700 cMpc 2 (∼ 76 pMpc 2 ) at z ∼ 6. To a 4σ detection (at ∼ 5.5 mJy at 850µm ), the SMGs overdensity is 0.68 +0.21 -0.19 (±0.19), exceeding (∼ 1.68 ×) the source counts of blank sky surveys (S2CLS; Geach et al. 2017). We found 13/54 quasars show SMG overdensities (at ∼ 5.5 mJy) in their fields, and the others do not. The combined field of these 13 quasars exceeds the blank field counts with the overdensity to 5.5 mJy of δ SMG ∼ 2.46 +0.64 -0.55 (±0.25) in the regions of ∼ 150 arcmin 2 ; but to the bright end (e.g., 7.5 mJy), the overdensity is not significant. Among them, three (3/13) z ∼ 6 fields show a significant excess over the blank field, with the overdensity of 2.1 -5.4 at 5.5 mJy.

We also compared our work and some previous environment studies traced by LAEs and LBGs on a similar scale. These LAEs or LBGs trace the young, dust-poor star-forming galaxies, while the SMGs trace massive, metal-enriched, and highly obscured objects. They may be correlated in the high-density region; therefore, more follow-up spectroscopic works to search for LAEs with JWST, HST, VLT, and ALMA. are required for these fields at a large scale. Here we present SMGs and SMG candidates in the quasar fields at z ∼ 6. The full catalog includes 170 SMGs selected from 450 and 850 µm images. Column (1) contains the SMG IDs. SMGs are detected above 3.5σ at least in one band, marked as 'ObjectShortName-SMG+index'. Columns (4) and ( 7) are the SMG peak flux density in 850µm and 450µm bands of SCUBA-2 images; Columns (5) and ( 8) are signal-to-noise ratio in each band, respectively. Columns (6) and ( 9) are the deboosted flux density for 850µm and 450µm detections calculated from the simulation. Column (10) is the angular distance between SMGs and their central quasar.

P007+04-SMG1 00h28m07.4s +04d57m41.7s 10.79±1.08 9.95 10.87±1.78

00h28m05.2s +04d56m53.7s 7.71±1.12 6.86 7.77±1.79

00h28m07.9s +04d57m17.7s 6.84±1.10 6.24 6.88±1.76

00h28m05.7s +04d56m29.7s 4.60±1.17

P007+04-SMG6 00h28m08.7s +04d59m01.7s 4.18±1.17 3.56 3.54±1.21

P007+04-SMG7 00h28m02.3s +04d59m01.7s 4.30±1.22 3.51 3.70±1.28

J0100+2802-SMG1 01h00m13.9s +28d00m25.8s 9.94±1.47 6.77 10.00±2.08

J0100+2802-SMG2 01h00m09.0s +28d04m13.8s 8.00±1.32 6.05 8.03±1.95

J0100+2802-SMG3 01h00m14.2s +28d02m53.8s 5.48±1.12 4.88 5.18±1.50

J0100+2802-SMG4 01h00m21.4s +28d02m21.8s 5.99±1.40 4.29 5.81±1.79

J0100+2802-SMG5 01h00m17.5s +28d03m17.8s 4.11±1.16 3.53 3.37±1.14

J0148+0600-SMG1 01h48m37.8s +05d58m28.0s 6.83±1.49 4.57 6.76±2.05

J0148+0600-SMG2 01h48m40.0s +05d59m36.0s 5.43±1.27 4.27 4.98±1.55

(10) J0148+0600-SMG4 01h48m42.9s +06d01m08.0s 5.86±1.43 4.11 5.56±1.79

J0148+0600-SMG5 01h48m30.8s +06d00m12.0s 5.23±1.42 3.67 4.72±1.61

J0148+0600-SMG6 01h48m41.0s +06d01m20.0s 4.79±1.32 3.64 4.12±1.39

P036+03-AGN 02h26m02.8s +03d02m59.4s 5.47±1.01 5.42 5.31±1.48

02h26m02.6s +03d04m19.4s 11.32±1.08 10.49 11.39±1.77

P036+03-SMG2 02h25m58.0s +03d03m47.4s 6.33±1.06 6.00 6.33±1.68

J0305-3150-SMG1 03h05m22.3s -31d49m47.8s 6.50±1.30 5.00 6.49±1.83

P055-00-SMG1 03h41m44.2s -00d48m28.6s 6.38±1.31 4.86 6.21±1.79

P055-00-SMG5 03h41m38.1s -00d46m36.6s 4.66±1.33 3.51 3.97±1.36

P056-16-SMG1 03h46m58.0s -16d29m44.8s 5.85±1.05 5.55 5.81±1.51

P060+24-SMG1 04h02m08.4s +24d51m32.5s 6.90±1.28 5.39 6.83±1.86

P060+24-SMG2 04h02m07.2s +24d50m08.5s 6.78±1.46 4.65 6.69±1.97

P060+24-SMG3 04h02m19.8s +24d52m20.5s 4.88±1.34 3.63 4.32±1.46

P089-15-AGN 05h59m46.4s -15d35m00.1s 3.56±1.17 3.05 2.60±0.96

P089-15-SMG1 05h59m42.8s -15d34m56.1s 6.28±1.17

J1152+0055-SMG2 11h52m23.4s +00d56m28.6s 4.16±1.13 3.69 3.48±1.17

.7s +04d20m47.6s 6.88±1.35 5.09 6.79±1.92

P217-07-AGN 14h31m40.7s -07d24m47.5s 6.03±1.17