Deep Submillimetre and Radio Observations in the SSA22 Field. II. Sub-millimetre source catalogue and number counts

We present the deepest 850 $\mu$m map of the SSA22 field to date, utilizing a combination of new and archival observations taken with SCUBA-2, mounted at the James Clerk Maxwell Telescope (JCMT). The mapped area covers an effective region of approximately 0.34 deg$^2$, achieving a boundary sensitivity of 2 mJy beam$^{-1}$, with the deepest central coverage reaching a depth of $\sigma_\text{rms}$ $\sim$ 0.79 mJy beam$^{-1}$, the confusion noise is estimated to be $\sim$ 0.43 mJy beam$^{-1}$. A catalogue of 850 $\mu$m sources in the SSA22 field is generated, identifying 390 sources with single-to-noise ratios above 3.5, out of which 92 sources exceed 5$\sigma$. The derived intrinsic number counts at 850 $\mu$m are found to be in excellent agreement with published surveys. Interestingly, the SSA22 number counts also exhibit an upturn in the brighter flux region, likely attributed to local emitters or lensing objects within the field. On the scale of $\sim$ 0.3 deg$^2$, the 850 $\mu$m number counts are unaffected by cosmic variance and align with the blank field. In the deep region ($\sigma_\text{rms}$ $\leqslant$ 1 mJy), the counts for fluxes below 8 mJy are consistent with the blank field, and the excess in the brighter regime is not significant. Due to the limited number of very bright sources and the insubstantial cosmic variance in our field, we attribute the fluctuations in the number counts primarily to Poisson noise. The SCUBA-2 850 $\mu$m detection in the SSA22 field does not exhibit indications of overdensity.


INTRODUCTION
The exploration of the formation of the first stars/galaxies and the evolution of galaxies is crucial for our understanding of the universe.In the 1980s, the Infrared Astronomical Satellite (IRAS) made a significant discovery, revealing that the bolometric luminosities of some galaxies were predominantly dominated by far-infrared emission (Sanders & Mirabel 1996;Simpson et al. 2019).Subsequently, the space-based Cosmic Background Explorer (COBE), aimed at studying the cosmic background radiation, found a substantial portion of the background radiation was emitted in the far-infrared and submillimetre regime (Puget et al. 1996;Fixsen et al. 1998;Hauser et al. 1998;Casey et al. 2014;Zavala et al. 2017).These observations of galaxies suggested that the interstellar dust absorbed ultraviolet radiation produced during intense star formation and reradiated it at far-infrared wavelengths (Hildebrand 1983;Dole et al. 2006;Simpson et al. 2019).Consequently, due to spectral redshift (Garratt et al. 2023), this radiation was detected in the submillimetre band, leading to the identification of a population of dusty star-forming galaxies known as submillimetre galaxies (SMGs).
Despite severe water vapour absorption in the far-infrared band, some submillimetre atmospheric windows remain survived, enabling the observation of SMGs.Additionally, submillimetre observations ★ E-mail:ypao@pmo.ac.cn benefit from a strong negative K correction, resulting in luminosities of submillimetre galaxies remaining nearly constant regardless of their redshift.Specifically, for the 850 μm detection, the radiation from thermal dust approaches the Rayleigh-Jeans limit, and as redshift increases, the rest-frame band shifts toward the peak of the spectral energy distribution (SED) of galaxy, compensating for cosmological dimming (Blain 2002;Dale & Helou 2002;Casey et al. 2014;Geach et al. 2017;Garratt et al. 2023).For instance, submillimetre galaxies exhibit similar flux densities whether they are located at a redshift of 0.5 or 10 (Blain & Longair 1993).Therefore, the 850 μm observation offers a significant advantage for tracing back to the epoch of cosmic reionization.
In the late 1990s, the Submillimetre Common User Bolometer Array (SCUBA) mounted on the 15 m James Clerk Maxwell Telescope (JCMT) at Mauna Kea captured the first glimpses of submillimetre galaxies (Smail et al. 1997;Barger et al. 1998;Hughes et al. 1998;Lilly et al. 1999).Subsequently, the new generation of the bolometer array camera, SCUBA-2, advanced our understanding of the properties and physical processes of SMGs.The redshift distribution of SMGs selected at 850 μm exhibit a peak ∼ 2 -3 (Chapman et al. 2005;Wardlow et al. 2011;Chen et al. 2013a;Zavala et al. 2017).It is essential to note that the redshift distribution of SMGs depends on the observation wavelength, with galaxies selected in longer wavebands tending to have higher redshifts (Chapin et al. 2009;Chen et al. 2013a;Roseboom et al. 2013;Zavala et al. 2014; Wang et al. 2017).Regarding their properties, these galaxies are heavily dusty and have masses on the order of ∼ 10 10−11 M ⊙ (Swinbank et al. 2004;Hainline et al. 2011;Michałowski et al. 2012), and contain large amounts of gas (Frayer et al. 1998;Greve et al. 2005;Ivison et al. 2011;Thomson et al. 2012).SMGs are luminous in the infrared, with luminosities on the order of ∼ 10 12−13 L ⊙ (Barger et al. 1998;Zavala et al. 2017;Hyun et al. 2023), and their star formation rates usually range between ∼ 100 -1000 M ⊙ yr −1 (Magnelli et al. 2012;Swinbank et al. 2014;Barger et al. 2014).Despite these observations, the exact formation mechanisms of SMGs remain unclear.Leading hypotheses include starbursts triggered by major merger phases of gas discs in galaxies (Baugh et al. 2005;Ivison et al. 2012) and galaxy interactions (Swinbank et al. 2010;Kartaltepe et al. 2012;Chen et al. 2015;Wang et al. 2017).Cosmological simulations have also inspired other formation mechanisms, such as long-period gas accretion (Narayanan et al. 2015) and disc instabilities (Lacey et al. 2016).
X-ray detections of submillimetre galaxies selected at 850 μm revealed that approximately 20 per cent of these galaxies host active galactic nuclei (AGNs) (Alexander et al. 2005;Almaini et al. 2005;Pope et al. 2008;Wang et al. 2017), some SMGs had been found to reside in quasars, which are further associated with supermassive black holes (SMBHs) (Jones et al. 2017;Alexander et al. 2008;Wang et al. 2013).Due to their nature and characteristics, SMGs were often considered as progenitors of massive elliptical galaxies in the local universe (Eales et al. 1999;Lilly et al. 1999;Genzel et al. 2003;Simpson et al. 2017).Observations had suggested that these galaxies were located in overdense regions (Dressler 1980;Menéndez-Delmestre et al. 2013;Dannerbauer et al. 2014;Arrigoni Battaia et al. 2018).While they appeared to trace the large-scale structure in the highz universe, their individual connection to clusters or protoclusters remained unclear (Umehata et al. 2014).
The SSA22 field, originally one of the four target fields of the deep imaging and spectroscopic survey, derived its name from 'Small Selected Area (SSA) at 22 h ' (Lilly et al. 1991;Cowie et al. 1994).Steidel et al. (1998) reported a highly significant concentration of galaxies at a redshift of z ∼ 3.09 in the SSA22 field, traced by Lymanbreak galaxies (LBGs).Subsequent studies of Lyman  emitters (LAEs) supported the existence of an overdensity region in this field, with the local density being 6 times the average, indicating a robust 'protocluster' (Steidel et al. 2000;Hayashino et al. 2004;Tamura et al. 2009;Yamada et al. 2012;Umehata et al. 2014).Umehata et al. (2019) also discovered a large-scale filamentary structure in the core of the SSA22 protocluster.Several studies supported the synchronous formation of SMGs and LAEs within the same cosmic structure, implying that SMGs might be also preferentially formed in dense regions (Tamura et al. 2009;Umehata et al. 2015), where star formation is fuelled and black holes are grown (Lehmer et al. 2009;Umehata et al. 2019).A population of extended LAEs was referred to as Lyman  blobs (LABs), which might preferentially reside in dense environments at high redshift, hinting at an association with massive galaxy formation (Matsuda et al. 2004;Ao et al. 2017).The origin of LABs remained unclear and might be diverse, but some observations suggested that they were associated with large-scale gas outflows driven by AGN or intense starbursts, and filamentary LABs might be linked to cold accretion streams from the surrounding intergalactic medium (IGM) (Matsuda et al. 2011;Geach et al. 2014;Ao et al. 2015Ao et al. , 2020)).
The paper is organized as follows.In Section 2, we provide a detailed description of our observations and data reduction process in the SSA22 field.Section 3 presents the SSA22 coverage map and outlines the source extraction procedure.We elaborate on the Monte Carlo simulation approach used to correct the extracted information of the submillimetre sources and present our source catalogue in Section 4. Finally, we report the number counts for the SSA22 field and discuss the results in Section 5.The article employs the cosmological parameters: H 0 = 67.8km −1 s −1 Mpc −1 , Ω Λ = 0.692, Ω  = 0.308.

Observations
A subregion of the SSA22 field was observed with SCUBA-2 on the JCMT between 2015 April and June (program ID: M15AI91; Ao et al. 2017).The submillimetre transmission quality heavily relied on the atmospheric extinction, quantified by the precipitable water vapour (PWV), which was obtained from a 183 GHz water vapour monitor (WVM) mounted on JCMT (Holland et al. 2013).The PWV was directly converted to the zenith opacity at 225 GHz ( 225 ), a crucial factor for extinction correction of the observation frequency (Dempsey et al. 2013;Mairs et al. 2021).The weather conditions during the observation ranged from 0.04 to 0.08, with an integral time of 22 per cent categorized as Grade 1 for  225 < 0.05 (PWV < 0.83 mm), and the remaining categorized as Grade 2 for 0.05 <  225 < 0.08 (0.83 mm < PWV < 1.58 mm).The quality of our data was also influenced by the observing elevation.To ensure an adequate exposure time for the map, an elevation restriction of under 70 • 1 was adopted.The elevation of our observations ranged from 44 • to 69 • .
The new observation map was centred at  = 22 h 17 m 31.s 70,  = +00 • 17 ′ 50.′′ 0 (J2000), with a total on-source mapping time of 19 h.To achieve a larger coverage area with a smoothly increasing noise gradient, we employed a rotating PONG pattern.In this observation mode, the telescope scanned across a predetermined sky region and 'bounced' off the edge.At the end of each 'bounce', the map was rotated, and the PONG pattern repeated the motion (Geach et al. 2017).This scanning method was available in various sizes (e.g., 900 ′′ , 1800 ′′ , 2700 ′′ , 3600 ′′ , and 7200 ′′ ) to meet different coverage requirements.For our observations, we used the PONG-900 observing strategy with a scan velocity of 280 ′′ s −1 , which repeated 11 times during an approximate 40-minute integration time 2 , covering a field of diameter 15 arcmin.The SSA22 field had been previously mapped as one of the target fields of the SCUBA-2 Cosmology Legacy Survey (S2CLS; Geach et al. 2017).The earlier observation (program ID: MJLSC02) adopted the PONG-1800 scanning mode at the centre of the SSA22 field, coordinated at  = 22 h 17 m 36.s 30,  = +00 • 19 ′ 22. ′′ 7 (J2000).It took 72 hours between 2012 September and 2013 December, with a sky opacity ( 225 ) ranging from 0.05 to 0.114.Our new observation map achieved an average 1 depth of 1.1 mJy beam −1 , which is comparable to the SSA22 archive map.Combining our new observations with the archival data from the Canadian Astronomy Data Centre (CADC), with a total exposure time of 91 h, allowed us to generate the deepest submillimetre 850 μm map of the SSA22 field to date.
The SCUBA-2 records information simultaneously from two wavebands: 450 μm and 850 μm.However, the 450 μm data are less stable for analysis compared to the 850 μm data due to its higher 1 https://www.eaobservatory.org/jcmt/instrumentation/continuum/scuba-2/observing-modes/#High_elevation_ constraints 2 https://www.eaobservatory.org/jcmt/instrumentation/continuum/scuba-2/observing-modes/ sensitivity to transmission fluctuations.Observations at 450 μm requires much drier conditions as the atmosphere is more opaque at this wavelength, and the background radiation in the sky is more severe than at 850 μm.Furthermore, the telescope surface is not optimized for 450 μm (Coppin et al. 2006), leading to a beam pattern that is not well described by a two-dimensional Gaussian and is very sensitive to even small JCMT dish deformations (Knudsen et al. 2008).As a result, the data quality of 450 μm is generally inferior to that of 850 μm, and 450 μm data are typically used as a supplement to the 850 μm data.For this study, we solely present the 850 μm data, catalogue, and results of the SSA22 field.

Data reduction
During the observation, SCUBA-2 records information from the astrophysical signal, atmosphere, ambient emission, and noise in the 5120 bolometers distributed across four sub-arrays.The raw scan data is a time-varying signal stream typically split into many subscans of about 30 seconds, which are saved in the STARLINK NDF format with a '.sdf' extension (Holland et al. 2013).
The SCUBA-2 850 μm data were reduced using the Dynamical Iterative Map Maker (DIMM) of the Sub-millimeter User Reduction Facility (SMURF) (Chapin et al. 2013) tool package in the STARLINK 2021A.We specifically used the 'makemap' command and the configuration file 'dimmconfig_blank_field.lis' that was designed to detect extremely low significant sources, and we set the pixel size to 1 arcsec.The map-maker started with pre-processing, independently cleaning the raw data for each observation.All segments of individual scans, approximately 40 minutes each, were connected to create continuous time-series.The time-streams were then subjected to a flat-field correction (Dempsey et al. 2013) that bracketed each observation and calibrated the data in units of pW using the associated fast-flat scan, and were re-sampled to match the desired map pixel size of 1 arcsec.Next, the spikes in the data were identified and removed if the bolometer values were higher than 10 times the local noise level within a box size of 50 timeslices (Chapin et al. 2013).These spikes had high amplitude and short duration in each bolometer.Sudden steps in the baseline were corrected, and any gaps were filled.Finally, the baseline of each time-stream was subtracted using polynomial fitting.To eliminate the low-frequency noise, a high-pass filter was applied, corresponding to 200 arcsec at 850 μm.This step was taken to avoid convergence problems caused by the high-pass filter during the iteration stage (Thomas et al. 2021).
When the cleaning of time-stream data was completed, the data reduction entered the iteration stage.Firstly, the common-mode signal, primarily dominated by variations in atmospheric emission and temperature, was removed by averaging values of working bolometers of the sub-array at each time-step (Chapin et al. 2013).Secondly, an extinction correction was applied based on the opacity monitored by the WVM.Next, the time-series was re-gridded to celestial coordinate, and the astronomical signal mode was produced and subtracted from the time-series data, which was then reprojected back to the time domain.In the final step, the noise mode for each bolometer was measured from the data with the removed signal mode, and this noise was used to weight each map in the mosaic step (Chapin et al. 2013;Geach et al. 2017;Thomas et al. 2021).The iterative program continued looping until the maximum iteration number of 4 was reached.
After obtaining a set of individual reduced maps from the map-maker, several additional post-reduction steps needed to be performed.We used the Pipeline for Combining and Analysing Reduced Data (PICARD) of STARLINK and the recipe CAL-IBRATE_SCUBA2_DATA to calibrate the data from pW to mJy beam −1 .Since all observations were completed before 2016 November 19, a new flux calibration factor (FCF) of 525000 mJy beam −1 pW −1 was applied simultaneously to the flux and rms layer of each map, as recommended in the new version 2021A of STARLINK (Thomas et al. 2021).This FCF had an uncertainty of 7 per cent (Mairs et al. 2021).The calibration method of beam FCF, also known as peak FCF, integrated the total flux density of a source distributed over the beam region into one pixel.We noticed that a single scan was skipped by the program because it only had 25 bolometers, and it was ignored in the remainder of the processes.Alternatively, if we used the advanced science pipeline ORAC-DR and the recipe REDUCE_SCAN_FAINT_POINT_SOURCES, the program would called the same make-map configuration parameters and FCFs by default to reduce the raw data, and the same resultant scan maps were output.
The PICARD recipe MOSAIC_JCMT_IMAGES was used to mosaic all calibrated maps pixel-by-pixel weighted by inverse-variance (Simpson et al. 2019).The default WCSMOSAIC approach was chosen for this purpose.Next, we applied a script to remove any pixels beyond 4 from the median of a box region centred on each pixel.This step is necessary to prevent bad pixels from forming bright spots after matched-filtering and making them difficult to distinguish from the detected source.Each pixel value of the mosaicked map was sampled from multiple bolometers, the root-mean-square (rms) of the bolometer values contributing to each pixel was the noise (Thomas et al. 2021).It was contributed by the instrument and atmosphere, but which we would later call as local instrument noise (Simpson et al. 2019).
The commonly used matched-filtering technique was employed to optimize source identification and characterization, and to suppress any residual large-scale noise that might not have been entirely rejected in the preceding filtering step (Thomas et al. 2021).For this purpose, we requested the recipe SCUBA2_MATCHED_FILTER.The matched-filtering procedure involved two steps: first, the mosaic was convolved with a wide Gaussian function with a full width at half maximum (FWHM) of 30 ′′ to obtain a smooth version of the mosaicked map, which was then subtracted from the original map to remove any residual large-scale structure.Next, the image was convolved with the point spread function (PSF).For a more detailed visualization of this process, refer to Chen et al. (2013a, fig.1).Finally, after applying the matched-filtering, we obtained the SSA22 field 850 μm map.The flux calibration check would be presented in Section 4.6.2.

Coverage map
We define the effective area as the total area used for source detection and subsequent analysis, where the local instrument noise (i.e.rms) is less than 2 mJy beam −1 , corresponding to approximately 0.34 deg 2 .The flux map of our 850 μm mosaic is shown in Fig. 1, the white circles represent significant sources with fluxes higher than 5 times rms, which are later used to construct the empirical point spread function.The white dashed lines in the figure depict the distribution of local instrument noise, due to the observation mode, the exposure time gradually decreases and the instrument noise gradually increases outward.Within the deepest coverage of the map, approximately 5 arcmin in radius, the rms level reaches a uniform value of ∼ 0.79 mJy beam −1 .The contours illustrate the variation in noise levels at 0.9, 1.0, 1.2, 1.4, 1.6, and 1.8 mJy, smoothly increasing from the centre outwards.The thick white dashed line delineates the region with  rms ⩽ 1 mJy, referred to as the 'deep region' in subsequent chapters.The white circles represent the 92 sources detected above 5, with a circle radius equal to the FWHM of the 850 μm empirical PSF, which measures 14 ′′ .On the other hand, the pixel distribution of the science map extends beyond the Gaussian noise on both the left and right sides.The right side reflects the submillimetre signal, while the left side results from the matched-filter process.
In Fig. 2, we present the pixel value histogram of the signal-tonoise ratio (S/N) map, which exhibits excess in both positive and negative directions relative to Gaussian noise (indicated by the grey dashed line or light grey shadow).The long tail of the S/N distribution, extending almost to 20, is shaped by the emission of astrophysical objects.The negative excess is contributed by the negative rings (see Section 3.2.1)around the source peaks introduced by the matched-filter in the final data reduction step.

Modelling PSF
The SCUBA-2 instrumental response to 850 μm light from distant galaxies can be expressed by a two-component Gaussian function, with an expected effective FWHM of ∼ 12. ′′ 6 (Mairs et al. 2021).However, during the matched-filtering reduction, the shape of the instrumental PSF undergoes changes.As a result, it becomes necessary to reconstruct the PSF to accurately extract the filtered submillimetre sources.
To reconstruct the PSF, we first identified all sources with significance greater than 5 over the effective area and extracted the pixel values within a radius of 40 ′′ of each source.Subsequently, we stacked the pixel values of every source and obtained the median point spread function.The reshaped profile of the 850 μm sources, with a FWHM of approximately 14 ′′ , is depicted in Fig. 3.This profile exhibits a circular depression resulting from the application of the matched-filter to the map and instrumental PSF.Notably, the width of the empirical PSF is slightly narrower than previously reported due to updates in the matched-filter (Mairs et al. 2021).

Source extraction
Dusty submillimetre galaxies in the early universe typically occupy compact regions (with extents of a few kiloparsecs; Umehata et al. 2017).Due to their small angular sizes compared to the SCUBA-2 average beam size (approximately 130 kpc at the typical redshift of ∼ 2-3 of SMGs; Hayward et al. 2013a;Lovell et al. 2021), these distant dust emission populations appear as unresolved point sources in SCUBA-2 survey maps (Casey et al. 2013;Shim et al. 2020;Zhang et al. 2022).In addition, the data reduction pipeline uses the default units of mJy beam −1 , meaning that the pixel value of each source peak in the map represents its total measured submillimetre flux density (see Section 2.2 for context on flux calibration).
The matched-filter fully optimizes the science map for source extraction and enhances the clarity of each source.For detection and source identification, we adopted a simple yet common 'top-down' iterative algorithm, following the approach used in Wang et al. (2004Wang et al. ( , 2017)).Firstly, we located the most prominent peak in the signal-tonoise ratio map and subtracted 5 per cent (referred to as 'CLEAN gain'; Wang et al. 2017) of the peak flux scaled PSF from a circular region with a radius of 40 ′′ , centred on this peak in the flux density map.We then recorded the noise level and the subtracted flux density, along with the peak pixel position, which served as the location of the detected source.This process was repeated, identifying the most significant peak, subtracting part of the peak flux, and recording the source information until the signal-to-noise ratio peak reached a predetermined detection threshold, set to 3 to ensure the adequate detection of potential sources.When identifying a peak within a radius of 7 ′′ (approximately half FWHM of the PSF) of a previously detected source, we considered them as a single source and continued subtracting at the original extracted location.Once the iterative procedure was completed, we combined each piece of cleaned flux density with the remaining flux below the threshold at the identified location to obtain the total flux density for each source.The final catalogue consisted of 390 sources selected from the primary catalogue, with a cut limit of 3.5.While the choice of 3.5 was somewhat arbitrary, it was commonly used in the literature, and the total false sources accounted for 17 per cent (see Fig. 8).
In the literature, the commonly adopted proportion for peak subtraction was 100 per cent.However, this method sometimes produced fake sources (Zhang et al. 2022), which was also observed in our tests.The 'CLEAN' method, involving a gradual deduction of peak flux, effectively separated somewhat blended sources.Comparing to the catalogue from Geach et al. (2017), this subtraction approach allowed us to distinguish two pairs of close submillimetre sources (with separations less than the FWHM).Besides, our method identified four additional pairs of close sources in the 3.5 cut-off catalogue, which could be valuable for follow-up observations and the search for counterparts.

THE MONTE CARLO SIMULATION
Because of the random measurement errors of the instrument during the reception of submillimetre photons, a certain luminosity of light can be 'scattered' to higher or lower values, resulting in a Gaussian or Poisson distribution.Moreover, within the detection capability of telescope, there are always more faint sources than bright ones.This situation leads to a higher number of events of dimmer sources being erroneously assigned higher flux densities during observations than the opposite case, presenting a persistent risk of flux overestimation.
The phenomenon of flux overestimation, also known as 'flux boosting', is influenced by both the Eddington and Malmquist biases discussed above.Also, source blending of multiplicity can contribute to the boosting effect.On the other hand, some of the faintest sources may be undetected due to being obscured by random noise, while some noise peaks may be misidentified for the sources because of random fluctuates of noise.Moreover, the relatively large FWHM of the SCUBA-2 beam, coupled with the pointing accuracy of JCMT of around 1-2 arcsec, introduces a considerable level of uncertainty in the positional accuracy of detected sources.
Given the challenges in accurately recognizing and extracting sources, conducting Monte Carlo simulations becomes essential to obtain a robust catalogue and reliable number counts, and facilitate follow-up observations of individual sources.

Simulation Method
First, we needed a jackknife map.We randomly divided the total of 134 calibrated scan maps into two groups.Each group consisted of half of the new data and half of the archived data, ensuring a balanced distribution of exposure time between the two subsets.Subsequently, we synthesised two half-mosaics separately and performed a subtraction between them based on the world coordinate system.The pixel values in the resulting subtracted map were then scaled using the formula √  1 ×  2 /( 1 +  2 ), where  1 and  2 represented the exposure times of the corresponding pixels in the two half-mosaics (Chen et al. 2013a).After applying additional steps to remove bad pixels and use matched-filter, the resulting map was designated as the jackknife map, also known as the true noise map (Chen et al. 2013a;Wang et al. 2017).Since the astrophysical signals from the two half-mosaics were eliminated during the subtraction process, the flux layer of the jackknife map contained only Gaussian instrument noise (as indicated by the grey-filled portion in Fig. 2).Thus, this map provided a sophisticated estimation of the true random noise in the science flux density map.The jackknife rms map showed a high level of consistency with the science rms map (Wang et al. 2017).
In the second step, we conducted a Monte Carlo simulation by injecting millions of artificial sources into our jackknife true noise map.As the fundamental model, we adopted the differential counts from Geach et al. (2017), represented by a Schechter function: Considering that the observed flux densities of the brightest and faintest sources detected in the science 850 μm map were approximately 3 and 15 mJy beam −1 , respectively, we set the flux density of the injected artificial sources to range from 1 to 20 mJy.Hayward et al. (2013b) indicated that this faintest flux setting below the detection limits of almost single-dish surveys.We divided this injected flux range into 20 bins in logarithmic space and calculated the theoretical number of sources in each flux density bin according to the counts model.For simulated flux densities greater than ∼ 13 mJy, the   small map area resulted in theoretical source numbers less than one.
To ensure sufficient sources in the bright flux density bin and obtain accurate statistics, we applied Poisson randomisation to determine the actual number of simulated sources in each bin.The flux density of each artificial source was then randomly and uniformly sampled from the corresponding flux bin.
For injecting the artificial sources into the jackknife map, we utilized the PSF outline of the filtered sources obtained in Section 3.2.1.Each source was placed randomly at any pixel in the map, regardless of any clustering effects (Simpson et al. 2019), and anywhere within the pixel.After injection, we performed the source extraction program described in Section 3.2.2 on the simulation map, repeating the injection and extraction procedure 5000 times for a robust statistical analysis.
Finally, we matched the injection and extraction source catalogues.In each simulation, if a source with a significance of ⩾ 3.5 from the extracted catalogue was found within a radius of 0.75 times the FWHM (i.e.0.75 × 14 arcsec ∼ 11 arcsec) of an injected source, we considered the injected source as 'recovered'.In cases where an extraction matched multiple injections, we retained the brightest component for the match.Once an injected source was successfully matched, it was excluded from further matching attempts.Sources from the injection catalogue that remained unmatched were considered to be submerged in noise and not extracted, while unmatched sources from the extraction catalogue were regarded as false sources resulting from Gaussian noise.

Flux boosting
To estimate the flux boosting, we compared the extracted and injected flux densities of each matched source, which depended on both the observed flux density and the local instrument noise.As sources with the same injected flux were affected by different local noise perturbations, their observed fluxes could vary, in turn, an observed flux has a corresponding intrinsic flux density probability distribution, denoted as p( true ).
Following an empirical method (Geach et al. 2017), we obtained p( true ) for each two-dimensional bin ( obs , ) by calculating the histogram of flux density of injected sources falling into that bin.This allowed us to estimate the intrinsic (deboosted) flux density of the detected sources in the real map.The median value of the injected flux distribution was used as the deboosted flux density of the source, and the standard deviation was considered as the uncertainty in the deboosting process (Shim et al. 2020;Zhang et al. 2022).An instance of p( true ) is shown in Fig. 4, where the centre of the injected flux density distribution gets closer to the observed flux as the observed flux increases, indicating a decreasing flux boosting effect with increasing observed flux density.Similarly, as the local noise decreases, the flux boosting becomes weaker.
Once the matching procedure of the 5000 simulated input and output catalogues was completed, we aggregated all two-dimensional bins and the relevant median results of the injected flux density distribution to create a two-dimensional image (Fig. 5, panel (b)), which was used for deboosting the sources extracted from the SSA22 850 μm map.Meanwhile, we summarized the statistical results regarding the injected source number, completeness, and positional uncertainty using a flux interval of 0.4 mJy and an rms interval of 0.04 mJy.
We used the meshed values of the panel as a look-up table to estimate the correction for each source, based on their observed flux density and rms.By locating the corresponding ( obs , ) bin for observed flux density and instrument noise of a source, we derived the deboosted flux and its uncertainty, which, combined with the instrument noise and confusion noise, provided the total uncertainty of the source in the real catalogue (Shim et al. 2020).
The average flux boosting effect could be described by a simple power-law function of the measured signal-to-noise ratios of sources (Geach et al. 2017;Simpson et al. 2019;Shim et al. 2020): In Fig. 6, the average flux boosting reaches an amplification factor of 1.63 for 3.5 sources, and the boosting curve decreases to about 1.04 at a S/N level of 15, with the average boosting factor for all sources being about 1.4.It is essential to note that this power-law curve is an approximation, and the flux boosting of each source depends on its observed flux density and local instrument noise, rather than being directly derived from the above equation (3).We conducted a simple comparison of flux boosting with results from literature.In the region where the S/N is significant, our flux boosting is roughly consistent with Zhang et al. (2022), slightly higher than Geach et al. (2017), Shim et al. (2022), and Garratt et al. (2023), and slightly lower than Hyun et al. (2023).We attribute these differences primarily to variations in the uniformity and coverage areas of surveys.Additionally, the presence of denser regions, where faint objects below the detection threshold crowd into the instrument beam, can lead to a higher boosting curve.The higher boosting values in the low S/N section reported by Shim et al. (2022) and Garratt et al. (2023) can likely be attributed to the higher noise levels of their maps.Differences in simulation performance may also contribute to the distinct variations in the boosting results (Scott et al. 2008(Scott et al. , 2010;;Wang et al. 2017;Zavala et al. 2017).

Completeness
Completeness is defined as the ratio of the number of recovered sources to the number of injected sources in the input flux density and instrument noise bin.To estimate the completeness correction for each source in the science map and the randomly sampled flux density of sources in Section 5.1, we used the two-dimensional meshed image (Fig. 5, image (c)).From the figure, it is evident that the efficiency of extracting simulated sources rapidly decreases when the input flux density is lower (< 8 mJy) due to the increase in rms, indicating that noise significantly influences the detection rate of sources in the low flux region.
The 50 per cent completeness limit corresponds to 2.6 mJy in the deepest region of the map and 6.2 mJy in the shallowest region.Similarly, the 80 per cent completeness limit corresponds to 3.4 and 8 mJy, respectively.It is essential to note that these corresponding flux values are the deboosted (intrinsic) flux densities of the sources.The average completeness shown in Fig. 7 has a simple relationship with the input flux density, achieving 50 per cent and 80 per cent completeness when the input flux density is 3.7 and 5.1 mJy, respectively.

Positional uncertainty
We also determined the positional uncertainties of the detected sources using a similar method.The positional uncertainty of a source is correlated with its injected flux density and local noise.The median positional errors are shown in graph (d) of Fig. 5 and can be described by a mathematical form similar to boosting: The positional offset gently decreases from 3.3 to 0.73 arcsec within the observed S/N range of 3.5 to 15, the position uncertainty for all sources is about 2.5 arcsec on average.The positional accuracy is usually influenced by the mapping depth, uniformity, and the survey coverage area.When using the positions of the SMGs detected in the 850 μm map for follow-up observations in other wavebands, caution is advised due to both the telescope pointing accuracy and the position uncertainty of the source extraction, which can reach about a few arcsec.

False detection rate
As mentioned earlier, fluctuations of Gaussian noise could occasionally be misidentified as sources.To measure the false detection rate (FDR), we generated 40 jackknife true noise maps, and compared the number of fake sources detected in these maps to the number of extracted sources in the real map within each S/N bin.
The signal-to-noise ratio pixel histogram (Fig. 2) illustrates that noise still contributes significantly at 3.5, accounting for up to 43 per cent of the detected sources (see Fig. 8).At 4.5, the differential false detection rate is about 4 per cent, and when the source significance reaches 5, the FDR drops to approximately 1 per cent.Therefore, caution is warranted when considering sources with low detection significance.Despite the presence of a large number of false sources at low detection significance, we chose an extraction threshold of 3.5 to achieve a more comprehensive source catalogue, as most of the sources were genuine.The cumulative false detection rate for sources with S/N > 3.5 was 17 per cent, meaning that out of all 390 sources, about 66 sources were false.Imitating the approach of Geach et al. (2017), we fitted the differential FDR with a piecewise linear function in logarithmic space: then we applied the fitting function to determine the differential FDR correction for each observed source.
An alternative approach to evaluate the FDR involved directly dividing the number of unmatched sources in the output catalogues of simulation by the number of sources extracted in the science 850 μm map within each S/N bin (Wang et al. 2017).The FDR results obtained from both methods were highly consistent since they were based on the same underlying principle: the unmatched sources primarily originated from Gaussian noise fluctuations in the jackknife map.

Confusion noise
Confusion noise is a persistent challenge in astronomical measurements (Scheuer 1957;Condon 1974;Takeuchi & Ishii 2004), arising from the combined effects of thermal emission from Galactic dust and crowding by unresolved faint extragalactic sources (Helou & Beichman 1990).In the context of this paper, cirrus confusion is negligible, and thus, it will not be considered further.Nevertheless, confusion noise remains an unmitigated factor, impervious to longer exposure times or increased detector sensitivity (Helou & Beichman 1990), affecting both position and photometry measurements (Hogg 2001).While not always rigorously accounted for in the SCUBA-2 signal map, we provide a simple estimation of its impact.
A commonly used rule of thumb for estimating confusion noise is that the confusion limit is reached when the source density reaches approximately one source per 20-30 resolution elements (Hogg 2001;Dole et al. 2003).Here, a resolution element is the beam area defined as Ω beam =  2 , with  ≈ FWHM/2.35 for the Gaussian assumption (Hogg 2001).Following Hogg (2001) and utilizing the final differential counts (Section 5.1), the calculated confusion limit is approximately 1.49 mJy, indicating that the source density reaches one per 30 beams when the flux is below 1.49 mJy.Consequently, at our 3.5 threshold, the estimated confusion noise is ∼ 0.43 mJy beam −1 (Dole et al. 2003;Simpson et al. 2019).
We followed Helou & Beichman (1990) and utilized the equation: where Ω bm represented the 850 μm beam area of SCUBA-2 with a value of 242 arcsec 2 (Geach et al. 2017), and S lim was the confusion limit flux of 1.49 mJy given by the rule of thumb above.Based on our final number counts, the estimated confusion noise was ∼ 0.50 mJy beam −1 , consistent with former.If we considered the equation ( 6) with S lim set to infinity, the derived value of 0.95 mJy was close to the confusion noise of 0.86 mJy proposed by Geach et al. (2017), but this value was influenced by bright sources (Cowie et al. 2017;Geach et al. 2017).Cowie et al. ( 2017) highlighted that the signal-to-noise ratio of the signal residual map exhibits excess over the Gaussian instrument noise due to confusion noise (see Cowie et al. 2017, fig. 3(a)).For a test, we used the equation from Cowie et al. (2017) and Simpson et al. (2019) as follows: where  total was the total measured noise, encompassing both confusion noise  c and statistical noise  stat from the signal map (Cowie et al. 2017).Here,  total and  stat represented the dispersion of the signal subtracted map excluded the 3.5 detected sources and the jackknife true noise map, respectively (Simpson et al. 2019).The jackknifing process removed the bright sources from the astronomical signal map, as well as the confusion noise caused by faint sources below the detection limit (Austermann et al. 2010), thus  stat came from Gaussian instrument noise indeed.Due to the dominance of instrument noise throughout the SCUBA-2 850 μm map and its nonuniform nature, increasing with radius due to the observation pattern (Simpson et al. 2019), we calculated the value only in the deepest region with  rms ⩽ 1 mJy beam −1 .The resulting value of ∼ 0.28 mJy beam −1 was in basic agreement with previous estimates, the confusion noise would be ∼ 0.30 mJy beam −1 when the entire map was considered in the calculation process.
Because of the description above, we estimated the confusion noise to be 0.43 mJy beam −1 .The total uncertainty of the deboosted flux density for each source is obtained as the square root of the sum of squares of confusion noise, instrument noise, and deboosting error (Shim et al. 2020;Hyun et al. 2023).

Source catalogue
We conducted a cross-match between the source catalogue with the S2CLS survey (Geach et al. 2017), and checked both the observed and deboosted fluxes of the matched sources.The sources in our catalogue with a significance greater than 4 were used for matching, with a matching radius of 7", resulting in a total of 107 matching sources.The discrepancy in the observed flux was less than 4.8 per cent, while the deboosted flux density was approximately 2.3 per cent.
At the deepest part of the field of view, where the two observations overlap most, the differences were 0.4 and 2.9 per cent, respectively.Employing the same matching range and radius, we also compared the flux of the SCUBA-2 SSA22 field with the deeper ALMA observations (ADF22; Umehata et al. 2018).We compared the total flux contribution of the corresponding ALMA SMGs to the deboosted flux of the SCUBA-2 source and found that they were in agreement.Given that ALMA SMGs are situated in the deepest part of our observed field of view, we believe that the contribution from faint sources below the detection threshold may be smaller.It is important to note that there was an error in the flux conversion due to using a conversion factor from 1.1 mm to 850 μm for the ALMA sources.Additionally, the absolute flux calibration uncertainty is 15 per cent for S2CLS survey (Geach et al. 2017) and 10 per cent for ADF22 (Umehata et al. 2018).In conclusion, we find that both the observed flux and deboosted flux are consistent with the submillimetre flux studied previously.
The final source catalogue consists of 390 submillimetre sources with a significance greater than 3.5, among which 92 sources reached a 5 level.The Table 1 presents the first 10 most significant detections, including information about source coordinates, measured and deboosted flux densities, signal-to-noise ratios, completeness, false detection rates, and positional offsets.The total uncertainty of deboosted flux combines the deboosting uncertainty, instrument noise, and confusion noise.For visual convenience, the values of FDR shown in the table are presented in logarithmic form.Caution should be exercised when using sources with high false detection risk and location uncertainty for follow-up observation and identification.
We plotted the distribution of deboosted sources in Fig. 9.The 'deep region' corresponds to the area with  rms ⩽ 1 mJy in the science map, which is the overlapped region between new and archival observations.The dark grey bars in the plot represent the number distribution of sources in this region, which contains 171 sources.Additionally, the light grey bars stacked at the top display the distribution of 219 sources in the remaining area.
Among the sources, five have flux densities greater than 10 mJy, three of which fall between 10 and 11 mJy.The other two are the brightest sources in the SSA22 850 μm map, with flux densities of 15.03 and 13.65 mJy, respectively.The former source is listed in the S2CLS catalogue presented by Geach et al. (2017) and is located in a protocluster core with a redshift of ∼ 3.09 (Tamura et al. 2010;Umehata et al. 2015) in the centre of the field.The latter source is found at the edge of the map and requires further investigation through multiwavelength studies to determine its nature, such as whether it is a local emitter, a lensing galaxy, a bright starburst galaxy, or a multiplicity.

Number counts
The surface number density of submillimetre galaxies serves as an effective test for galaxy formation models and provides valuable constraints on cosmological simulations.The submillimetre survey covers a larger area, then the map will have a smaller cosmic variance and the number counts will be closer to unbiased.A JCMT's ongoing large program, the SCUBA-2 Large eXtragalactic Survey (S2LXS, Garratt et al. 2023), is making significant progress in submillimetre survey area.
To construct the number counts for the SSA22 field, we followed a three-step process.First, we deboosted the flux densities of sources using the two-dimensional diagram of the simulation result (Section 3.2.2) to obtain the intrinsic flux density.Second, we calculated the surface density of each source by dividing one by the effective area and bin width, measured in units of deg −2 mJy −1 (Shim et al. 2020).And we applied correction to each source based on its deboosted flux density, which involved weighting each source according to its authenticity and completeness: where   represents the effective area, and d  denotes the width of the flux bin in which the deboosted flux density of each source is located.F  and C  stand for the false detection rate and completeness correction, respectively (Shim et al. 2020;Zhang et al. 2022).To obtain the raw number counts, we summed up the weighted surface number densities of all sources within each bin.
Finally, we adopted a technique from Geach et al. (2017) to refine the raw number counts.Due to the presence of instrument noise, different injected flux densities could be measured as the same flux density, as explained in Section 4.1.Hence, we performed 1000 random samplings of the flux density for each source from its corresponding 'intrinsic' (injected) flux density distribution generated by the Monte Carlo simulation (see Section 4.2).For each sampled flux density, the completeness correction values were computed based on the flux value and the local noise of the source.Conversely, the differential FDR corrections remained the same for all sampling fluxes of a particular source, as they depended solely on the observed S/N.The surface number density within each flux density bin was calculated from the 1000 sampling processes.The mean of each flux bin was then considered the final differential number counts, with the standard deviation providing the uncertainty.
We also constructed cumulative number counts based on this 'smoothing' process.This cumulative representation offers a valuable tool for quantifying the overall number density distribution of the 850 μm sources and facilitates comparison with other submillimetre surveys or models.The number counts, ranging from 2 to 16 mJy, were summarized in Table 2. Here, S represents the bin centre, and S ′ represents the left edge of the flux bin.The cumulative counts described the number densities of sources with flux greater than S ′ , with their uncertainties derived from the 1000 sampling processes.
The real number counts are illustrated in Fig. 10, with black solid circles representing the counts for the entire effective area.The bestfitting line with the Schechter form for the differential number counts is indicated by the black solid line, featuring the following parameters: N 0 = 7785 ± 1095 deg −2 , S 0 = 2.64 ± 0.2 mJy beam −1 , and  = 1.65 ± 0.1.As anticipated, our results closely align with the average number counts from Geach et al. (2017).This finding suggests that the number counts of the SSA22 field remain similar to those of blank fields, even in the presence of dense regions.While the Schechter model offers more physical insight (Geach et al. 2013), we also provide the best-fitting results using the double-power law form: Both fitting results are consistent with each other, despite their different functional forms.The parameters for both fittings are listed at the end of Table 2.
We observed an upturn in our number counts above 10 mJy, with the two brightest sources mentioned earlier (Section 4.6.2) being the primary reason for the count points in the tail being higher than the prediction of the fitting curve.According to Geach et al. (2017), this upturn could be influenced by local sources or gravitational lensing galaxies.In the deep region ( rms ⩽ 1 mJy; area ∼ 0.09 deg 2 ), represented by the light grey spots, we observed an apparent boost between 8 and 10 mJy (also see Fig. 9), aside from the excess caused by the known brightest source.

Comparison with other surveys
In Fig. 10, we compared our number counts with several published submillimetre surveys.Our differential number counts largely overlap with the S2CLS survey (Geach et al. 2017), which covers an area of approximately 5 deg 2 and has a sensitivity of about 1.2 mJy.This indicates that it is reasonable to adopt their Schechter counts model for our Monte Carlo simulation.The S2LXS survey (Garratt et al. Table 2.The number counts for the SSA22 field and the best-fitting parameters of differential counts.The flux density S refers to the centre of the flux bin, while S ′ denotes the left edge of the flux bin.The column labelled d /d represents the differential number counts, and  ( >  ′ ) denotes the cumulative number counts.The uncertainties of the number counts arise from the standard deviations obtained from 1000 times of sampling.The last two rows display the best-fitting parameters for the differential counts.The first row presents the parameters for the Schechter form fitting, and the second row provides the parameters for the double power law form.2023), which provides number counts of the rare bright source population in a uniformly scanned sky area of about 9 deg 2 , also shows consistent results with ours within the flux density range we have presented.Moreover, the number counts of the SSA22 field are in close agreement with the source number density of the S2COSMOS (Simpson et al. 2019) and NEPSC2 (Shim et al. 2020), both of which map a sky zone of about 2 deg 2 .The S2COSMOS survey reaches a depth of approximately 0.5 mJy over a limited portion of the large survey field, while the EGS deep field map (Zavala et al. 2017) has a central instrumental noise of 0.2 mJy.And the sensitivity in the deepest region of the JWST-TDF survey (Hyun et al. 2023) reaches 0.8 mJy , similar to ours, and their survey estimations are essentially consistent with our counts, except for a slight excess at brighter end.
On the other hand, we also present the differential number counts of the lensing field A370 (Chen et al. 2013a), which show values two or three times higher than our observations at flux densities ≳ 2 mJy.Chen et al. (2013b) and Hsu et al. (2016) combined cluster lensing fields with blank fields, benefiting from the remarkable amplification of the gravitational lensing effect, extending significantly into the extreme faint flux regime of ∼ 0.1 mJy.Our statistical distribution and the extrapolation of the Schechter model at the faint end agree well with their consequences.For flux densities > 5 mJy, the lensing effect becomes evident, showing enhancements of ten to tens of times when the flux exceeds 10 mJy.This phenomenon is consistent with the conclusions reported by Negrello et al. (2010), Vieira et al. (2010), andAretxaga et al. (2011), where a clear upturn is observed in the counts at the brighter flux range due to gravitational lensing amplification from foreground structures or cosmic variance.
Meanwhile, we presented and compared the ALMA 870 μm, 1.1  (Casey et al. 2013), S2CLS (Geach et al. 2017), and S2COSMOS (Simpson et al. 2019).Through extensive comparison with results from other surveys, the number counts in the SSA22 field exhibit significant concordance with previously published counts.In the cumulative counts figure, we also present predictions of 850 μm derived from various galaxy formation and evolution models, depicted as dot-dashed curves.These models encompass the semi-empirical model (Hayward et al. 2013a), the semi-analytical model GALFORM (Cowley et al. 2015), the empirical model SIDES (Béthermin et al. 2017), and the cosmological hydrodynamical model SIMBA (Lovell et al. 2021).Except for the model proposed by Hayward et al. (2013a), the remaining models account for the source blending effect introduced by the relatively large beam of 850 μm submillimetre observations.The model by Béthermin et al. (2017) constrains the SFR to values less than 1000 M ⊙ yr −1 and achieves simulation outcomes that align well with observational data.We present both the source extraction number counts of their model and the intrinsic counts denoted by a blue dotted line.
mm, and 1.2 mm interferometric counts.To transform the counts to the 850 μm band for visual comparison, we used the equation ( 7) from Dunne & Eales (2001) which relate the flux ratio at different (sub)millimetre bands to thermal dust emission and dust emissivity.For the dust emissivity , we adopted a value of 1.8 (Planck Collaboration et al. 2011;Dudzevičiūtė et al. 2020).Since the submillimetre regime is on the Rayleigh-Jeans side of the Planck function (Dunne & Eales 2001), the second term in the equation ( 10) is simply ∝  2 .Therefore, we used the scale conversion factor of  3.8 to convert the interferometric results to the 850 μm band for comparison.The expected S 850 /S 1100 ratio was approximately 2.66, similar to conversion factors obtained from using the free spectral index parameter (Austermann et al. 2010) and assuming a simple empirical relation between the stellar formation rate and the gas-phase metal mass (Hayward et al. 2013a).The multiplicative factor of ∼ 1.09 converted 870 μm counts to corresponding 850 μm flux bins, close to the value derived from the SED template assuming submillimetre galaxies located at z = 2 (Béthermin et al. 2017).The earliest report from Karim et al. (2013) found that interferometric counts were generally consistent with results from lowresolution single-dish surveys, but falling below 850 μm single-dish counts for flux regime less than 7 mJy and greater than 10 mJy.This suggested that most single-dish submillimetre sources were not caused by a mixing of unresolved SMGs, while bright sources were significantly affected by source blending.Stach et al. (2018) pointed out that interferometric observations were not affected by blending, and the counts dropped by 28 per cent compared to single-dish surveys.They and Simpson et al. (2020) concluded that about half or more of single-dish sources with 850 μm flux greater than 9-12 mJy were multiple sources.They also believed that about a third of these populations were physically related, which was consistent with the summary from galaxy formation models (see Section 5.2.2;Hayward et al. 2013a;Lovell et al. 2021).Umehata et al. (2017) obtained the ALMA number counts of the SSA22 central region of 6 arcmin 2 (ADF22).We plotted their counts excluding the sources of the protocluster in Fig. 10 left panel.Based on this, Umehata et al. (2018) extended observations to obtain a continuous ALMA interferometer map of 20 arcmin 2 , whose showed several times more counts compared to the blank field.Hatsukade et al. (2018) presented the 26 arcmin 2 ALMA blank field survey in the GOODS-S field, and their estimates were consistent with other single-dish surveys or interferometric counts.
Overall, our number counts were found to be in good agreement with previously published counts observed from submillimetre single-dish and interferometer surveys, as well as the lensing fields.

Comparison with models
The replication of high-redshift submillimetre galaxy counts has constituted a pivotal focus within galaxy formation models (Lovell et al. 2021).In this section, we undertake a comparison between our cumulative counts and model predictions.Hayward et al. (2013a) employed a combination of a semiempirical model, three-dimensional hydrodynamical simulations, and dust radiative transfer to predict SMGs counts.They also considered the contribution of individual subpopulations of SMGs, such as merger-induced starbursts, galaxy pairs, and isolated disc galaxies.Their simulated single-dish counts matched total cumulative counts at S 850 < 5 mJy and were in good agreement with our deep region counts ranging from 3 < S 850 < 12 mJy (refer to Fig. 10).Despite not accounting for the lensing effect, their projections aligned with the bulk of observed counts in the high flux regime (> 10 mJy).Hayward et al. (2013b) further considered and reevaluated the source blending due to a large beam size, finding that over half of blended sources had at least one spatially independent blending component, with a typical redshift separation of about 1. Lacey et al. (2016) introduced an enhanced iteration of GAL-FORM, a semi-analytical framework for galaxy formation and evolution.This updated model integrated a mildly top-heavy initial mass function (IMF), AGN feedback, a dust radiation model, and observational constraints to attain a more physically motivated representation.In the model, starbursts predominantly arose from disc instabilities.Utilizing the GALFORM model and employing the light-cone technique, Cowley et al. (2015) simulated the 15 arcsec beam effect of JCMT, instrument noise, and matched filters to yield the 'submillimetre map' at 850 μm.They subsequently employed the 'top-down' approach for source extraction.The counts derived from the simulated map, depicted in Fig. 10, slightly exceeded the observed counts, manifesting an excess by a factor of ∼ 1.3-1.5 within the flux range of < 8 mJy.Model counts closely matched our counts at higher flux levels.
The empirical SIDES model (Béthermin et al. 2017) used updated version of the two star-formation modes (i.e.main sequence star-forming and starburst galaxies; see Béthermin et al. 2012), combined with observational constraints, realizing a 2 deg 2 extragalactic submillimetre survey simulation.The model imposed an upper limit of 1000 M ⊙ yr −1 on the star formation rate (SFR).We plotted the simulated source extraction and intrinsic number counts with the SFR limit in Fig. 10.Apart from extremely bright flux densities (> 15 mJy), the source extraction counts remain consistent with observed counts across the flux density range.Between S 850 = 5-10 mJy, the model overpredicts cumulative counts by around 50-60 per cent.The intrinsic counts curve is lower than our observed counts by ∼ 30-50 per cent for S 850 < 10 mJy, aligning with interferometric observation results (Simpson et al. 2015;Stach et al. 2018).At the brighter end (> 10 mJy), intrinsic counts are in agreement with single-dish observations.Lovell et al. (2021) harnessed the cosmological hydrodynamic simulation SIMBA and dust radiation transfer to generate 850 μm emission.To consider the source blending effect, they projected the radiation from light-cone onto a 0.5 deg 2 sky plane, which was then convolved with the SCUBA-2 PSF to generate a simulated singledish 'observed map'.Their predicted cumulative counts lie a factor of about 1.3-2 below observed counts at S 850 = 2-5 mJy, and 2-3 times lower than our measurements at S 850 = 5-10 mJy.

SMGs overdensity ?
In this section, we aim to analyse the variation in our field compared to the blank field and quantify the excess of submillimetre sources in the SSA22 field.Before proceeding, we will first estimate the cosmic variance of the SSA22 field.Utilizing the configuration proposed by Moster et al. (2011), we adopted a mean redshift of z = 3.09 and a redshift bin width of Δz = 1.The cosmic variances for typical submillimetre galaxies (with logM ★ /M ⊙ ∼ 10.5) in the total map and the deep region are approximately 19 per cent and 26 per cent, respectively.Consequently, even within the central region, the cosmic variance remains insignificant.When considering solely the extremely narrow redshift slice at redshift 3.09 (Δz = 0.1), cosmological variance becomes somewhat significant.However, due to the broad range of redshifts (z ∼ 1 -6) observable in the submillimetre waveband, and our extensive coverage, the overall detected cosmic volume is vast, diminishing the impact of cosmic variance (Somerville et al. 2004;Driver & Robotham 2010).Based on these results, it appears that the large-scale structure of SSA22 has a limited impact on the counts, suggesting that the fluctuations in the counts are likely primarily due to Poisson randomness.Simpson et al. (2019) suggests that cosmic variance does not strongly affect the 850 μm population on scales of ∼ 0.5-3 deg 2 .Therefore, it is reasonable to regard the average number counts of S2CLS survey (∼ 5 deg 2 ) as a representative of blank field, which we use to estimate the excess of SSA22 850 μm counts.To conduct the comparison, we divided the astronomical map into two regions: the 'deep region' with  rms ⩽ 1 mJy, located at the centre of the map and containing the protocluster core, and the 'shallow region' (1 mJy <  rms < 2 mJy) covering the rest of the area.These two regions occupy ∼ 0.09 and 0.25 deg 2 , respectively.Following the approach used by Geach et al. (2017), we used the equation Δ = (N -N S2CLS ) / N S2CLS to assess the elevation in source density, where N and N S2CLS represented the cumulative number counts of the SSA22 field and the S2CLS survey, respectively.The value of Δ represents the cosmic variance of number counts of each region compared to the blank field.
In Fig. 11, the dotted lines represent 50 per cent of the S2CLS average surface number density.The counts in the shallow region are well consistent with the 50 per cent ranges of the S2CLS source density, indicating that the shallow region is comparable to a blank field.The total number counts in the SSA22 field are also good consistent with the blank fields in the flux range of < 10 mJy.In the deep region, for fluxes greater than 10 mJy, the counts seem to exceed those of the blank field by a considerable margin.However, given the limited number of extremely bright sources in this range, the impact of randomness is relatively significant.Additionally, as discussed in Section 5.1, this upturn may be influenced by local or lensing galaxies.Furthermore, the impact of cosmic variance in the deep region is relatively limited.Therefore, we believe that the apparent excess in this flux regime may not necessarily indicate a significant influence of large-scale structure on the counts and this variation is more likely due to Poisson noise (Geach et al. 2017).In the remaining part of the deep region counts, there is only a slight excess at around 9 mJy, we conducted a simple comparison of the source number density with the cumulative curves presented in Geach et al. (2017) and Simpson et al. (2019).We calculated the source number density N SSA22 / N blank for a specific flux cut, where N SSA22 represents the source number density in our field, and N blank represents that from the blank field.At S 850 ⩾ 4 mJy, the ratios for the total and deep regions are 0.97 ± 0.06 and 1.08 ± 0.10, respectively.At S 850 ⩾ 9 mJy, these ratios become 1.03 ± 0.20 and 1.77 ± 0.48.field compared to the blank field.The red inverted triangles, blue squares, and green circles represent the deviation of cumulative number counts in the deep region, shallow region, and total region, respectively, relative to the mean density of the S2CLS survey.The dotted lines correspond to 50 per cent of the S2CLS average source density.Compared to the blank field counts, the deep region exhibits an obvious excess, while the shallow region shows good agreement with the S2CLS survey.The total number counts lie somewhere in between, and it is reasonable to attribute this discrepancy to the concentration of sources in the deep region at the field centre.
Following the approach of Arrigoni Battaia et al. (2018) and Zhang et al. (2022), we quantified the excess of the 850 μm number counts by comparing them to the blank fields represented by S2CLS (Geach et al. 2017) and S2COSMOS (Simpson et al. 2019).We refitted our differential counts and derived the best-fitting parameter N 0 using their parameters S 0 and  of the Schechter counts model.When fixing S 0 = 2.5 and  = 1.5 (Geach et al. 2017), we found N 0, refit / N 0, S2CLS ratios of approximately 1.33 and 1.04 for the deep and shallow counts, respectively.Fixing S 0 = 3.0 and  = 1.6 (Simpson et al. 2019), the N 0, refit / N 0, S2COSMOS ratios were approximately 1.38 and 1.08 for the deep and shallow counts.For direct comparison with known overdensity fields, we selected smaller central areas with  rms ≲ 0.84 and 0.95 mJy (∼ 0.036 deg 2 and 0.074 deg 2 ) in the SSA22 field, which roughly corresponded to the sizes of individual fields presented by Arrigoni Battaia et al. (2018) and Zhang et al. (2022).In these smaller regions, the density excesses in the counts compared to blank fields were about 48 per cent and 25 per cent, respectively, much lower than the estimations of 4 times and 2 times the number density of the blank field mentioned in the literature.In brief, the counts in the deep region are 1.35 times those of the blank field, while the total counts are only 1.16 times.
Overall, within our submillimetre waveband detection volume, cosmic variance is not substantial.The 850 μm number counts, derived from a total detection coverage of approximately 0.3 deg 2 , conform to the blank field.In the deep region, the number counts for flux density < 8 mJy are in agreement with the blank field.Due to the lower number of sources at the brighter end of the flux density and the insignificant cosmic variance, we believe that the excess at the brighter end is mainly influenced by Poisson noise.The quantified excess of the deep region counts over the blank field is 1.35 times, which is not statistically significant.In the SSA22 field, the detection of SCUBA-2 sources at 850 μm does not exhibit signs of overdensity.

SUMMARY
In this work, we present the deepest 850 μm mosaic map of the SSA22 field to date.The map is a combination of JCMT/SCUBA-2 new observations and archival data spanning from 2012 September to 2015 June, totaling 91 hours of observations.The sensitivity of the map is limited to less than 2 mJy beam −1 , resulting in a total effective area of about 0.34 deg 2 .The deep region, with an  rms ⩽ 1 mJy, covers approximately 0.09 deg 2 , while the central deepest coverage reaches an rms of ∼ 0.79 mJy beam −1 , which is about 1.5 times deeper than the previous S2CLS survey.We estimate the confusion noise to be 0.43 mJy beam −1 .
We process the raw data and adopt an iterative 'top-down' method to extract the potential submillimetre sources.The extraction threshold is set to 3.5, meaning that 17 per cent of sources are false detections.The final source catalogue contains 390 sources, with 92 sources reaching 5.This SSA22 850 μm source catalogue is publicly available.We perform a Monte Carlo simulation 5000 times using the jackknife true noise map to estimate the effects of flux boosting, positional uncertainty, completeness, and false detection rate.Our flux boosting and positional uncertainty are larger than S2CLS due to differences in mapping uniformity, survey area, and cosmic variance.The mean value of boosting factors for all sources is about 1.4, and the average positional uncertainty is about 2.5 arcsec.The average completeness at 50 and 80 per cent levels is 3.7 and 5.1 mJy, respectively.At the 4 level, the cumulative false detection rate is about 5 per cent, which drops to approximately 1 per cent when the source significance reaches 4.5.
We construct the SSA22 850 μm number counts, showing good agreement with some large-area SCUBA-2 surveys, such as S2CLS, S2COSMOS, and NEPSC2 (Geach et al. 2017;Simpson et al. 2019;Shim et al. 2020), within the deboosted flux range of 3-15 mJy.Additionally, the counts are broadly consistent with published 850 μm measurements from single-dish and interferometer surveys, and lensing field, as well as predictions from galaxy formation models.We observe a clear upturn in the counts in the brighter flux regime.
Compared to the blank fields, we measure that the total effective area counts are 1.16 times higher, and the counts in the deep region show an improvement of about 1.35 times.In the flux density range of < 8 mJy, the deep region closely resembles the blank field.Due to the limited number of sources at the brighter end and the insignificant impact of cosmic variance from the large-scale structure, we attribute the excess of deep region counts in the brighter flux regime mainly to Poisson noise.In general, the total number counts align with the blank field, and the excess of the deep region counts relative to the blank field is not significant.Furthermore, we find that the number counts do not indicate overdensity.
Given the abundant multiwavelength data sets available for the SSA22 field, including X-ray, optical, infrared, and radio bands, we plan to conduct research on the physical properties of the submillimetre sources in the SSA22 field, and to explore the evolution process of high-redshift galaxies and the structure of the early universe.
the "Light of West China" Program (No. xbzg-zdsys-202212).The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institute; the National Astronomical Research Institute of Thailand; Center for Astronomical Mega-Science (as well as the National Key R&D Program of China with No. 2017YFA0402700).Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities and organizations in the United Kingdom and Canada.Additional funds for the construction of SCUBA-2 were provided by the Canada Foundation for Innovation.Y.Z.acknowledges financial supports from the Agencia Estatal de Investigación del Ministerio de Ciencia e Innovación (AEIMCINN) under grant (La evolución de los cíumulos de galaxias desde el amanecer hasta el mediodía cósmico) with reference (PID2019105776GB-I00/DOI:10.13039/501100011033) and the China Scholarship Council (202206340048).

Figure 1 .
Figure1.The presented flux density map covers an area of approximately 0.34 deg 2 with  rms ⩽ 2 mJy.The contours illustrate the variation in noise levels at 0.9, 1.0, 1.2, 1.4, 1.6, and 1.8 mJy, smoothly increasing from the centre outwards.The thick white dashed line delineates the region with  rms ⩽ 1 mJy, referred to as the 'deep region' in subsequent chapters.The white circles represent the 92 sources detected above 5, with a circle radius equal to the FWHM of the 850 μm empirical PSF, which measures 14 ′′ .

Figure 2 .
Figure2.The pixel distribution of the signal-to-noise ratio map of the SSA22 field (black solid line).The light grey filled area corresponds to the pixel distribution of the jackknife S/N map, while the grey dashed line represents the Gaussian distribution.Their good alignment indicates that the random instrument noise in the 850 μm map follows a Gaussian distribution.On the other hand, the pixel distribution of the science map extends beyond the Gaussian noise on both the left and right sides.The right side reflects the submillimetre signal, while the left side results from the matched-filter process.

Figure 3 .
Figure3.The profile of the reconstructed PSF model for the 850 μm mosaic map.The grey dots represent the normalized median profile of all sources with significance greater than 5.The black solid line depicts the best fitting curve of the empirical PSF, with the thick grey vertical line flagging the half of FWHM of ∼ 6. ′′ 95.The green dotted curve represents the instrumental PSF profile fromMairs et al. (2021).The circular depression around the summit of the source arises from the matched-filter, which has been updated(Mairs et al. 2021), resulting in a slimmer empirical PSF compared to previous literature.

Figure 4 .
Figure 4.An example of the injected (intrinsic) flux density probability distribution p( true ) (grey shaded bars) corresponding to the two-dimensional bin ( obs , ), obtained through the empirical recovery method from 5000 simulations.The black solid line represents the observed flux distribution, assuming a Gaussian distribution centred on the observed flux, with the standard deviation being the local instrument noise.The observed flux increases from left to right, and the rms decreases from top to bottom, with their values shown in the upper left corner of each stamp.As the observed flux increases and the local instrument noise decreases, the effect of flux boosting becomes less significant.

Figure 5 .
Figure 5.The graphical statistical outcomes from 5000 Monte Carlo simulations.Panel (a) illustrates the number distribution of injected sources, written in log 10 (  ).This distribution provides insights into the mapping sensitivity in relation to the area covered.The presence of horizontal spurs with rms values of approximately 0.8 and 1.15 mJy suggests slightly larger covered areas.Panel (b) presents the input flux density distribution as a function of the output flux density and noise, which is used for deboosting each source with a significance of ⩾ 3.5 in the science map.Panel (c) represents the completeness of the source extraction, which is related to the input flux density and the instrument noise.It is calculated by dividing the observed number of artificial sources by the total number of sources injected into the jackknife map.Panel (d) displays the positional error between the extracted position and the actual location of the 850 μm source, and is a function of the input flux density and the local noise.

Figure 6 .
Figure 6.(Left image)The median flux boosting is shown as a power law function of measured S/N, representing the ratio between recovered and intrinsic flux density of the same source.The '+' markers symbolise the median boosting result from the simulation, and the black solid line is the best-fitting curve.The horizontal black dotted line marks the value of 1.05.The boosting curve shows a significant drop within the range of S/N < 6, followed by a gentle decline to about 1.04.Overall, the boosting effect in our map falls in between other studies.(Right image) The positional uncertainty is defined as the distance between the extracted position and the input position of the same source.The 'x' markers represent the median offset result, and the best-fitting of positional error is depicted in a power-law form (black solid curve).The horizontal black dotted line indicates the value of 0.7.The positional error decreases steadily from 3.3 to 0.73 arcsec with the increase of observed S/N.

Figure 7 .
Figure7.The average completeness as a function of input flux density, is obtained through comparing the recovered sources number to the number of injected artificial sources.The triangle markers symbolize the average completeness at 3.5, the corresponding values of 50 and 80 per cent is 3.7 and 5.1 mJy, respectively.In this work, we only analyse the sources of 3.5 and above, the completeness of 4 displayed in graph is for comparison and shows the change in completeness under different measurement significance.

Figure 8 .
Figure8.The false detection rate is defined as the ratio of the number of 'sources' detected in the random noise of the jackknife map to the number of detected sources on the science map.The red and blue solid lines represent the best-fitting lines for the differential FDR, corresponding to two linear functions in logarithmic space.At the detection significance level of 3.5, the differential and cumulative FDRs are 43 and 17 per cent, respectively.Notably, the differential FDR sharply decreases to zero as the source detection significance reaches 5.2.

Figure 9 .
Figure 9.The number distribution of the 850 μm deboosted flux density for the 390 detected SCUBA-2 sources great than 3.5.The dark grey bars represent the flux density distribution of 171 sources in the deeper region, where the rms is ⩽ 1 mJy in the astronomical map.On the other hand, the light grey bars display the flux density distribution of 219 sources in the shallower region.

Figure 11 .
Figure11.The cosmic variance of 850 μm number counts in the SSA22 field compared to the blank field.The red inverted triangles, blue squares, and green circles represent the deviation of cumulative number counts in the deep region, shallow region, and total region, respectively, relative to the mean density of the S2CLS survey.The dotted lines correspond to 50 per cent of the S2CLS average source density.Compared to the blank field counts, the deep region exhibits an obvious excess, while the shallow region shows good agreement with the S2CLS survey.The total number counts lie somewhere in between, and it is reasonable to attribute this discrepancy to the concentration of sources in the deep region at the field centre.

Table 1 .
The 850 μm catalogue contains 390 submillimetre sources greater than 3.5.Among these, 92 detections achieve a significance of 5.Below are the top 10 most remarkable sources from the catalogue, sorted by significance.The column header S obs 850 ±  inst is the measured flux density and the local instrument noise.The column S deb 850 ±  tot represents the deboosted flux density and the total uncertainty, which is consisted of the deboosting uncertainty, instrument noise and confusion noise.The columns C, log 10 F,   give the source completeness, differential FDR in logarithmic form and position offset, respectively.The full catalogue is available in the online journal.
The left panel displays the differential number counts, while the right panel illustrates the cumulative counts.The black dots correspond to the number counts encompassing the entire effective area, accompanied by a solid black line representing the best-fitting Schechter curve.The grey dots depict findings from the central dense region with  rms ⩽ 1 mJy.Hollow polygons portray SCUBA-2 850 μm number counts from single-dish surveys, while circles represent counts from fields influenced by the lensing effect.Branch markers indicate interferometric counts of ALMA at wavelengths of 870 μm, 1.1 mm, and 1.2 mm, scaled to the 850 μm band using a factor of  3.8 .Colorful dashed lines denote best-fitting counts from representative surveys conducted in the COSMOS field