CHORUS. I. Cosmic HydrOgen Reionization Unveiled with Subaru: Overview

To determine the dominant sources for cosmic reionization, the evolution history of the global ionizing fraction, and the topology of the ionized regions, we have conducted a deep imaging survey using four narrow-band (NB) and one intermediate-band (IB) filters on the Subaru/Hyper Suprime-Cam (HSC), called Cosmic HydrOgen Reionization Unveiled with Subaru (CHORUS). The central wavelengths and full-widths-at-half-maximum of the CHORUS filters are, respectively, 386.2 nm and 5.5 nm for NB387, 526.0 nm and 7.9 nm for NB527, 717.1 nm and 11.1 nm for NB718, 946.2 nm and 33.0 nm for IB945, and 971.2 nm and 11.2 nm for NB973. This combination, including NB921 (921.5 nm and 13.5 nm) from the Subaru Strategic Program with HSC (HSC SSP), are carefully designed, as if they were playing a chorus, to observe multiple spectral features simultaneously, such as Lyman continuum, Ly$\alpha$, C~{\sc iv}, and He~{\sc ii} for $z=2$--$7$. The observing field is the same as that of the deepest footprint of the HSC SSP in the COSMOS field and its effective area is about 1.6 deg$^2$. Here, we present an overview of the CHORUS project, which includes descriptions of the filter design philosophy, observations and data reduction, multiband photometric catalogs, assessments of the imaging quality, measurements of the number counts, and example use cases of the data. All the imaging data, photometric catalogs, masked pixel images, data of limiting magnitudes and point spread functions, results of completeness simulations, and source number counts are publicly available through the HSC SSP database.


Introduction
Understanding cosmic reionization is one of the most important objectives in observational cosmology. Measurements of the Thomson scattering optical depth of free electrons in the intergalactic medium (IGM) in the cosmic microwave background with Wilkinson Microwave Anisotropy Probe and Planck satellites place the epoch of hydrogen reionization at z ∼ 10 (Komatsu et al. 2011;Planck collaboration 2018). Detection of the Gunn-Peterson trough in the spectra of quasars at z > 6 (e.g., Fan et al. 2006;Becker et al. 2007) and the decrement in the luminosity functions (LFs) of Lyα emitters (LAEs) at z > 6 (e.g., Kashikawa et al. 2006;Ouchi et al. 2010;Hu et al. 2010;Kashikawa et al. 2011;Nakamura et al. 2011;Konno et al. 2014;Matthee et al. 2015;Santos et al. 2016;Itoh et al. 2018;Hu et al. 2019;Taylor et al. 2020) suggest a rapid increase in the hydrogen neutral fraction, xHI, beyond z = 6. The xHI measurements obtained so far still has a large scatter due to the small statistics and the systematic uncertainty in each measurement (e.g., Robertson et al. 2015;Greig & Mesinger 2017;Inoue et al. 2018;Finkelstein et al. 2019;Naidu et al. 2020). New xHI measurements near z ∼ 7 based on large samples of galaxies over large cosmic volumes are required to reach a concordance in the history of reionization.
Galaxies, such as LAEs and Lyman break galaxies (LBGs), at z > 6 can be the dominant ionizing sources if their Lyman continuum (LyC) emissivity or escape fraction fesc, is sufficiently high. The required fesc is ∼ 20%, assuming a standard population synthesis model (e.g., Inoue et al. 2006;Robertson et al. 2015;Finkelstein et al. 2019;Naidu et al. 2020). Measuring the LyC emissivity of these galaxies, which is still quite uncertain, is the most critical step in the determination of the ionization sources for reionization. Bright quasars do not contribute much to the LyC emissivity for reionization owing to their rarity (e.g., Bianchi et al. 2001), whereas faint Active Galactic Nuclei (AGNs) might contribute to or even dominate (Giallongo et al. 2015;Madau & Haardt 2015). However, other observations suggest that faint AGNs are not sufficiently abundant (Kashikawa et al. 2015;Onoue et al. 2017;Matsuoka et al. 2018;Parsa et al. 2018). Measuring the LyC emissivity of AGNs is also required to resolve their role in reionization (Micheva et al. 2017;Grazian et al. 2018). The metal-free stellar population, the so-called Population III (Pop III) stars (e.g., Carlberg 1981), can be an important ionizing source, particularly in the early phase of the reionization process (e.g., Sokasian et al. 2004). Because their characteristics, such as the initial mass function (IMF), are unknown observationally owing to the lack of information (Nagao et al. 2005;Nagao et al. 2008;Kashikawa et al. 2012;Vanzella et al. 2020), identifying Pop III clusters and revealing their nature are highly relevant to understanding their role in reionization.
Depending on the dominant ionizing source, the topology of reionization is expected to be different: "inside-out" (Iliev et al. 2006) by galaxies or "outside-in" (Miralda-Escudé et al. 2000) by X-ray sources, such as AGNs. Observationally, the spatial inhomogeneity of xHI at z ∼ 6 has been reported to be based on quasar spectra (Becker et al. 2015) and LAE LFs (Ouchi et al. 2010;Nakamura et al. 2011). Obtaining the overdensities of LAEs at z > ∼ 6.5 may support the concept that these LAEs are located in the ionized bubbles produced by the galaxy overdensities (e.g., Castellano et al. 2016;Bagley et al. 2017;Higuchi et al. 2019;Harikane et al. 2019;Tilvi et al. 2020). Resolving the ionization topology is one of the main scientific goals for future H I 21 cm tomography with the Square Kilometer Array (SKA) (e.g., Hasegawa et al. 2016). However, with a sufficiently large survey of LAEs and LBGs, we may visualize the xHI map before the SKA era. Such topology observations will be important consistency checks of the ionizing source observations, by comparing the observed topology with the prediction based on the identified dominant ionizing sources.
The Hyper Suprime-Cam (HSC) Komiyama et al. 2018;Kawanomoto et al. 2018;Furusawa et al. 2018) has the widest field-of-view (FoV), 1.75 deg 2 , on an 8mclass telescope. This unique capability with a well-considered set of narrowband (NB) filters, for the first time, allows dealing with the central questions of reionization by an unprecedented wide and deep NB survey. Since 2007, we have studied the specifications of HSC NBs carefully (see section 2) and developed them under strict scientific reviews with the financial support of approximately 100 M yen from the Japan Society for the Promotion of Science. Four NBs (NB387, NB816, NB921, and NB101) among them are included in the on-going HSC Subaru Strategic Program (SSP) (Aihara et al. 2018a) for deep surveys of LAEs at z ≃ 2.2, 5.7, 6.6, and 7.3 . The latter three high-z surveys provide new estimates of the cosmic average of xHI at z ≃ 6.6 and 7.3 from comparison with LAE LFs Shibuya et al. 2018a;Shibuya et al. 2018b;Konno et al. 2018;Inoue et al. 2018). However, the HSC SSP cannot address the LyC emissivity and the xHI topology directly.
Thus, we initiated a Subaru intensive program, Cosmic HydrOgen Reionization Unveiled with Subaru (CHORUS), with NB387, NB527, NB718, IB945, and NB973 (Table 1; see  Fig. 1). Supplying the data of the NB921 and broadband (BB) filters from the HSC SSP, we aim to answer the questions of the ionization source and the topology, in addition to the history, by measuring the LyC emissivity values of galaxies and AGNs, abundance of Pop III stars, and LF of faint AGNs, and visualizing the xHI map directly see sections 6.1-6.5 for more information). These measurements are realized only in combination with our NBs. This program will become a legacy survey at least for a decade because no other deep and wide-area survey, such as the Large Synoptic Survey Telescope, has a set of NBs like this program.
This paper presents an overview of the CHORUS project and its first public data release: CHORUS PDR1. In section 2, we describe the CHORUS filter set in detail. Section 3 is a summary of the observations, data reduction, and photometric catalogs of the CHORUS data. Section 4 presents a summary of imaging quality, such as the size of the point spread function and the limiting magnitude of each NB. In section 5, the number counts in the five CHORUS NBs are given as a set of fundamental measurements of the NB imaging data. Section 6 gives a brief instruction of the usage of the CHORUS PDR1 data and some example science cases. A summary of this paper is found in the final section.
We assume a standard set of cosmological parameters of H0 = 70 km s −1 Mpc −1 , ΩM = 0.3 and ΩΛ = 0.7 and use the AB magnitude system (Oke 1990) throughout the study.

Filter set
The CHORUS filter set consists of four NBs (NB387, NB527, NB718, and NB973) and one IB (IB945) in addition to NB921 from HSC SSP. Table 1 is a summary of the characteristics of the filters and the targeted spectral features of high-z galaxies and the redshift ranges. Fig. 1 shows the total efficiency of the NBs as well as HSC BBs. 1 The CHORUS filter set not only covers a wide redshift range of 2 < ∼ z < ∼ 7 of Lyα but also enables us to observe multiple spectral features at specific redshifts at the same time. As found in Fig. 1, the LyC and He II 1640 line for z ≃ 3.3 LAEs selected with NB527 can be observed with NB387 and NB718, respectively, and similarly, the LyC, C IV line, and He II line for z ≃ 4.9 LAEs selected with NB718 can be observed with NB527, NB921, and NB973, respectively. 1 In the UD COSMOS field of the HSC SSP, there are two filters each for the rand i-bands, namely r2/r and i2/i . The newer r2 and i2 filters were used about 60% in the data according to the HSC PDR2 web page (https://hscrelease.mtk.nao.ac.jp/doc/index.php/data-2/). These two filter images were all combined to make final coadd images of the rand i-bands (Aihara et al. 2019    For z ≃ 2.2 LAEs selected with NB387, the He II line can be observed with NB527 (see Table 1). The last CHORUS filter, IB945, is a filter to select LBGs within a redshift range similar to that of LAEs selected with NB921 (z ≃ 6.6) and NB973 (z ≃ 7.0) (see section 6.4 for its science case). It is also useful to trace the line-free continuum between the C IV and He II lines at z ≃ 4.9 (see Fig. 1). In the following, we further describe the philosophy of the HSC SSP and CHORUS filter specifications and the characteristics of the CHORUS filters measured in the laboratory.
The NB filter set for the HSC imaging survey was designed very carefully ). The filter wavelengths in the redder wavelength range were determined to pass through OH sky windows at 816, 921, 973, and 1010 nm. The three filters, except for the one at 973 nm, were used for the NB imaging observations in the HSC SSP survey (Aihara et al. 2018a) called SILVERRUSH . The remaining NB973 filter became a CHORUS filter. The filter wavelengths in the blue wavelength range were determined to ensure that multiple spectral features of the LAEs selected by an NB filter were captured by the other NB filters and strong rest-frame optical emission lines of the LAEs passed through the OH sky windows in the near-infrared range maximally. The selected wavelengths were 387, 527, and 718 nm, which became the CHORUS filters. 2 This combination allows us to observe LyC, Lyα, and He II at 2 < ∼ z < ∼ 5 as described above. The full-width-at-half-maximum (FWHM) of each filter was also determined carefully ). The FWHM of NB921 was first determined, to cover the full width of the OH sky window. The FWHMs of the other NB filters, except for NB973 and IB945, were subsequently determined to ensure the same efficiency for the Lyα line detection as NB921.
Because the observed equivalent width of a line is EW obs = (1 + z)EW0, where EW0 is the intrinsic equivalent width and z is the source redshift, the FWHMs were determined by the relation, F W HM = F W HM921 × (1 + z)/(1 + 6.58). Note that 6.58 is the redshift of the LAEs selected with NB921. NB973 has a narrower FWHM than the relation to enhance the efficiency of the Lyα line detection and compensate for the shallower depth in the band. IB945 has an FWHM to fill the wavelength gap between NB921 and NB973.
The HSC filter is 600 mm in diameter, making it a challenging task to realize a uniform multi-layer interferometric coating 2 Although NB387 was also used for the 'Deep' layer observations in the HSC SSP (Aihara et al. 2018a;Ouchi et al. 2018), there was no data taken in the footprint of the 'UltraDeep' layer in the COSMOS field, where the CHORUS observations were performed.
to produce narrow bandpass filters with sharp cut-on and cutoff wavelengths. Therefore, there was a small deviation from the specifications during the production process. We examined the performance of the NB filters produced by the manufacturers, Optical Coatings Japan (NB387, NB527, NB718, NB921, and NB973) and Materion (IB945), using the method described in Kawanomoto et al. (2018). Slight non uniformities in the central wavelength (CW), FWHM, and transmission peak (TP) were found along the radial direction. The CW tended to increase with the radius, except for NB387 where the CW variation was a complex way. The FWHM also tended to increase in the outer part, except for NB387 where it decreased with the radius. The overall CW variation was less than 0.3-0.5% (the worst case was found in NB527), whereas the FWHM variation is ranging from a few-10% (NB387). The TP also tended to increase with the radius, except for NB387 again, where the TP decreased with the radius. The overall TP variation was less than 1-10% (NB973). A higher degree of uniformity was found in the azimuthal direction. Table 1 presents the CWs, FWHMs, and TPs for the area-weighted mean transmission functions of the CHORUS filters. The values for the transmission functions at the center positions of the filters are also provided for reference. The machine-readable filter transmission curves are available at the Subaru/HSC web site. 3

Observations, data reduction and photometric catalogs
We were awarded 13 nights in the classical mode of the Subaru Telescope for four consecutive semesters from S16B to S18A (i.e., 2 years) under a Subaru open-use intensive program, S16B-001I (PI: A. K. Inoue). We were also awarded 11.5 hours in the queue mode in semester S18B under an open-use normal program, S18B-004 (PI: A. K. Inoue), to supplement the time loss due to the weather during the original four semesters. Under the two programs, we performed deep imaging observations using our NB filters on the HSC on the dates listed in Table 2. The on-source exposure time as well as the effective exposure time used for the final coadd images are also listed in Table 2. Some of the observation dates were not under a photometric condition (i.e., low transparency), and we discarded a fraction of the observed data owing to bad quality. The observing field was the COSMOS field (Scoville et al. 2007), which has the same footprint as the HSC SSP Ultra Deep (UD) layer (Aihara et al. 2018a). Any sky region in the HSC SSP survey is defined by a set of two numbers of tract and patch (Aihara et al. 2018b). Most of the HSC SSP UD COSMOS is covered by tract number 9813. There are 9×9 sub-regionspatchesin a tract. Each patch has a sky area of 12 ′ × 12 ′ and a tract has that of 1.8 • × 1.8 • . The numbers of patches are expressed 3 https://subarutelescope.org/Observing/Instruments/HSC/sensitivity.html as 000 (south-east corner), 001, 002, ..., 008 (north-east corner), 100, 101, ..., 808 (north-west corner). The patch configuration in the UD COSMOS tract 9813 is found in Fig. 2. The data reduction was conducted using the HSC SSP pipeline, hscPipe, version 6.7, which is the same as that used for the SSP PDR2 data processing. We used the obtained imaging data satisfying the following conditions: exposure time = 1200 sec, seeing < 1.5 ′′ , and transparency > 0.4 for NB387, exposure time = 1200 sec (but 3 frames of 468-1016 sec), seeing < 1.2 ′′ , and transparency > 0.75 for NB527, exposure time = 1200 sec, seeing < 1.2 ′′ , and transparency > 0.76 for NB718, exposure time = 1200 sec, seeing < 1.2 ′′ , and transparency > 0.9 for NB973, and exposure time = 720 sec, seeing < 1.2 ′′ , and transparency > 0.7 for IB945. hscPipe generates several multiband photometric catalogs, which we call "CHORUS PDR1 catalog." The source detection is based on the signal-tonoise ratio (S/N) within the PSF against the background noise fluctuation (Bosch et al. 2018) and the detection criterion is 5σ, the same as the SSP PDR papers (Aihara et al. 2018b;Aihara et al. 2019). The readers should refer to pipeline paper (Bosch et al. 2018) and the SSP PDR2 paper (Aihara et al. 2019) for the detailed procedures of the data reduction, source detection, and photometry. In the following, we describe some specific points for the CHORUS data release.
The photometric zero-points for the CHORUS NBs were determined based on the BB photometry of the Galactic stars in the Pan-STARRS1 (PS1) catalog (Chambers et al. 2016) similar to the method for the HSC SSP (Aihara et al. 2018b;Aihara et al. 2019). The color terms translating PS1 BB magnitudes to HSC NB magnitudes were estimated based on the Pickles stellar spectral templates (Pickles 1998). This seems work well, except for NB387, where a systematic zero-point offset was found in Liang et al. (2020) (and see the last paragraph in this section). The multiband photometry was conducted with the so-called priority order (Bosch et al. 2018) of the filters, which was defined as i, r, z, y, g for the BB filters, followed by the order of NB921, NB816, NB973, NB718, NB527, NB1010, NB387, and IB945 for the NB filters. The priority order of the filters was used to determine the object positions. 4 For the line emitters detected in a single NB image but not detected in any other BB and NB images, photometry on the NB image was conducted at the positions of the detections in the NB image. However, for the line emitters that were also detected in other BB (and NB) images, the photometry positions were forced to be the positions in the most prioritized band among the images in which the objects were detected. In such cases, photometry on the NB image, which includes the line flux, can be underestimated if there are spatial offsets between the prioritized band and NB positions. This would not be a problem for the bright objects that are detected well both in the prioritized band and NB images; however, it might be an issue for the objects that are marginally detected in the prioritized band and whose positions in the band have a large uncertainty. Quantifying this effect would need a large set of pipeline simulations that was reserved for future line emitter analyses.
There are several different photometric magnitudes in the CHORUS PDR1 multiband catalog. In this paper, we do not recommend using cmodel magnitudes, which represent the total flux densities obtained by the brightness profile fitting (Bosch et al. 2018) and were considered representative magnitudes in SSP PDR1 (Aihara et al. 2018b). This is because there is an issue of erroneously bright cmodel magnitudes for some objects in SSP PDR2 (Aihara et al. 2019) and a follow-up study (Hayashi et al. 2020). According to section 6.6.3 of the SSP PDR2 paper (Aihara et al. 2019), this problem is likely to be caused by a failure in the deblending procedure in crowded regions of objects. This problem can be severe in the HSC SSP UD footprint in the COSMOS field, where the CHORUS observations were performed, because the very deep BB depth yields a very high density of the detected objects.
We instead recommend using the undeblended convolvedflux magnitudes for estimates of total magnitudes. 5 These measurements are based on the aperture magnitudes of the final coadd images smoothed with a Gaussian function (see section 3.6 of the SSP PDR1 paper, Aihara et al. 2018b). The Gaussian smoothing size is defined as the resultant FWHM of the point sources in the smoothed images, which is 0.59 ′′ , 0.84 ′′ , 1.1 ′′ , or 1.3 ′′ (Aihara et al. 2019). This process produces a set of the smoothed BB and 5 There is another similar measurement called convolvedflux, for which the source deblender algorithm in hscPipe works.
On the other hand, the source deblending process was turned off for undeblended convolvedflux. We have found that the two magnitudes of undeblended convolvedflux and convolvedflux are identical in most cases.
NB images with the selected same PSF FWHM. Subsequently, the flux densities of each object in all the smoothed (i.e., PSF-matched) images are measured in a circular aperture of diameter size 1.1 ′′ , 1.5 ′′ , 2.0 ′′ , or Kron-size, and aperture corrections are also applied assuming a point-source. 6 This process is performed for all the detected objects as a function of hscPipe (Aihara et al. 2018b;Aihara et al. 2019). These PSF-matched magnitudes yield better colors in crowded fields (Aihara et al. 2018b), which is the most important feature to reliably select LAEs in two-color diagrams. This procedure is very similar to the method performed in previous studies with SExtractor (Bertin & Arnouts 1996). In this study, we adopt 1.5 ′′ (2.0 ′′ )-diameter apertures on 0.84 ′′ (1.1 ′′ )-smoothed images. 7 There are two significant known issues in the SSP PDR2 photometric catalog (Aihara et al. 2019), and they are prevalent in the CHORUS PDR1 catalog also produced by hscPipe 6.7. The first one is that the uncertainties of the PSF-matched magnitudes are underestimated (section 5.8.11 of the SSP PDR1 paper, Aihara et al. 2018b). This occurs because covariances are introduced by the Gaussian smoothing process. This problem persists in the SSP PDR2 release (Aihara et al. 2019). Therefore, we do not recommend using these cataloged values of uncertainties for the PSF-matched magnitudes. Instead, we estimated the uncertainties by measuring the rms values of the aperture flux densities in the actual images. The limiting magnitudes provided in the next section are based on these measurements. 8 The second issue is that there is a systematic offset 6 The amount of the aperture correction can be found as undeblended convolvedflux X Y apcorr in the catalog database, where X and Y are the numbers indicating the target convolution size and the aperture size, respectively. 7 If the native PSF of a band is larger than the target smoothing FWHM, any Gaussian smoothing is not applied to the band and simply the native images are used for the aperture photometry (see section 3.6 of Aihara et al. 2018b). Such a case happened for NB387 in the 0.84 ′′ target FWHM smoothing. 8 The uncertainty in a magnitude is given by δm = (2.5/ln10)/(S/N ). The signal-to-noise ratio, S/N , of each object can be estimated from the ratio of undeblended apertureflux in the catalog database to the 1σ aperture in the zero-point of NB387 of the SSP PDR2 release (Aihara et al. 2019). This is due to the metallicity effect of the 4,000 A break in stellar spectra on N B387 − g color. hscPipe uses the Pickles spectral templates based on Solar metallicity stars (Pickles 1998) for zero-point estimations. However, the actual stars used for the calibration seem to be dominated by metalpoor halo stars because these stars are faint (g ∼ 20) and the HSC survey fields are located in high Galactic latitudes. This template mismatch causes a systematic offset of the zero-point only for NB387. A complete description and full analyses are found in appendix of Liang et al. (2020). We did not correct the NB387 magnitudes in the CHORUS PDR1 catalog but recommend users to apply a 0.45 mag subtraction correction to the NB387 magnitudes in the catalog, as recommended by the HSC SSP team. 9 4 Imaging data quality

Definition of masked areas
We defined the masked areas in each CHORUS NB image based on the flags from hscPipe and the visual inspection conducted by one of the authors. Fig. 2 shows the masked area of the NB527 image as an example. Pixels in each image that satisfied the following conditions of the hscPipe flags were selected as masked pixels: pixelflags bright object is True (pixels affected by bright objects) or pixelflags saturatedcenter is True (pixels affected by count saturation). In addition to these pixels, we also masked the pixels affected by the halos of flux density calculated from the limiting magnitudes in the next section. 9 https://hsc-release.mtk.nao.ac.jp/doc/index.php/known-problems-2/#hsclink-10 bright stars out. Generally, these halos are more prominent in NB images than in BB images probably owing to the multilayer interferometric thin film coating of NBs. We supplementary set larger circular (sometimes elliptical) masks to cover these bright stars' halos if the pixels with the above flags do not cover sufficiently large areas. One of the authors (S.Y.) selected these pixels by eyes. The edge of the FoV was also defined by eyes by the same author to avoid lower S/N regions. The effective survey area of each NB was calculated from the total number of non masked pixels in the image and is listed in Table 2. The exact positions of the masked pixels are available as image files in a FITS format.

Sizes of point spread functions
We measured a representative FWHM of the PSF in each patch of the CHORUS NB images using SExtractor 10 ver 2.19.5 (Bertin & Arnouts 1996). We adopted a median of the FWHM IMAGE measured by SExtractor for stellar-like objects, which were not in the masked areas, as the representative FWHM of the PSF. Stellar-like objects were selected depending on the measurements of SExtractor as those satisfying the following three criteria: (1) MAG AUTO in a range of [18.0, 22.0] for NB387 and NB527, [18.0,21.0] for NB718 and NB973, and [19.0,21.0] for IB945, (2) ELONGATION < 1.2 (i.e., nearly circular shape), and (3) FLAG < 4 (i.e., not located in a bad position, such as saturated or truncated pixels; see Bertin & Arnouts 1996 for the details). If there were insufficient number of stellar-like objects in a patch because of an extremely small non masked area, we did not measure the FWHM of the PSF. The measured FWHMs and their spatial variations are shown in Fig. 3, and the values of the area-weighted average and the central patch, 404, are listed in Table 2. The FWHM values range from 1 ′′ to 0.6 ′′ , depending on the band. The spatial variation in the FWHM values in each NB is as small as < 0.02 ′′ , demonstrating the excellent stability of the image quality across the entire FoV of the HSC.

Limiting magnitudes
We measured the limiting magnitudes in each patch of the CHORUS NB images by conducting aperture photometry at random sky positions with Python packages, Astropy 11 /photutils 12 . A large number of random positions (10,000 as default, and 5,000 or 1,000 at the edge of the FoV depending on the available area) were prepared in each patch by avoiding the masked pixels and the pixels of the objects detected by SExtractor in advance. If the non masked area in a patch was less than 1 arcmin 2 , we did not measure the lim-  iting magnitudes. For convenience, we measured the limiting magnitudes with two aperture sizes (1.5 ′′ and 2.0 ′′ in diameter) in each patch of each CHORUS NB image. We estimated the background brightness locally by adopting a median count in an annulus of 1.5 ′′ (2.0 ′′ ) width and 2.5 ′′ (3.0 ′′ ) inner diameter set around the 1.5 ′′ (2.0 ′′ ) aperture and subtracted it from the aperture count as a correction for the background contribution. The standard deviation, σ, was obtained from the histogram of the background-subtracted aperture counts by fitting a Gaussian function. Note that we did not apply any aperture correction to the counts. The 5σ values of the area-weighted average and those in the central patch, 404, for the two aperture sizes are listed in Table 2. The spatial variation in the 5σ values for the 1.5 ′′ -diameter aperture case are shown in Fig. 4. We can see a radial dependence of the limiting magnitudes; it is ∼ 1 mag shallower at the outer edge compared with that at the center.
The patch, 403, shows a shallower depth than the radial trend, which was also reported in Hayashi et al. (2020). This is probably because a half of the four channels in a CCD chip (SDO-ID=0 20/DET-ID=33) close to the patch center was unavailable since November-December 2016 (see the CCD information page of HSC on Subaru Telescope website 13 ) and the dithering amount was too small to compensate the low-sensitivity area.

Completeness
We estimated the detection completeness based on a simulation of embedding and recovering numerous mock galaxies in actual CHORUS NB images using hscPipe, Astropy, and GALSIM (Rowe et al. 2015) 14 . These mock galaxies were assumed to have an intrinsic brightness profile described by a circular Sersic profile with index n = 1.5 and a half-light radius of 1 pix of the HSC images (0.17 ′′ ). The adopted radius corresponds to 1.4-0.9 physical kpc for 2 < z < 7 and very consistent with those of typical LBGs/LAEs (Shibuya et al. 2015;Kawamata et al. 2018). When embedding the mock galaxies into each patch image, we convoluted the brightness profile with the PSF of the image. Because the intrinsic sizes of the mock galaxies were very small compared to the PSF size, the embedded profiles were similar to the PSF but slightly extended. The locations to embed the mock galaxies were chosen randomly, not avoiding the actual objects, to consider the cases in which the extended and bright galaxies overlap and hide faint ones. We considered nine intrinsic magnitudes of the mock galaxies in a range of ±2.0 magnitudes with a 0.5 magnitude step around the 5σ limiting magnitude. The number of the input mock galaxies was 500 for each intrinsic magnitude in each patch of each band. We ran hscPipe to detect the mock galaxies and measure the magnitudes. The criteria for successful recovery were detec-13 https://www.naoj.org/Observing/Instruments/HSC/ccd.html 14 https://github.com/GalSim-developers/GalSim tion within a radius of 0.5 ′′ from the embedded location and a magnitude difference of less than 0.5 magnitude. We adopted the undeblended convolvedflux magnitude (0.84 ′′ FWHM convolution and φ1.5 ′′ aperture case, except for NB387, where 1.1 ′′ FWHM convolution and φ2.0 ′′ aperture case is adopted) as the output total magnitude for Fig. 5; however, the results with the convolvedflux magnitude were the essentially same. Note that the undeblended convolvedflux and convolvedflux magnitudes are aperture-corrected by a PSF model and equivalent to the total magnitude for the point sources. The detection completeness is defined as the ratio of the number of successfully recovered galaxies to that of embedded mock galaxies. The mock galaxies placed in the masked areas were excluded from the analysis in this final step. To moderate the computation time, we conducted this experiment in only 4 patches of 202, 303, 404, and 606, instead of all the 81 patches. Fig. 5 shows the results of the simulations of the detection completeness. We find that the completeness is > ∼ 60% at the 5σ limiting magnitude (φ2.0 ′′ not aperture-corrected) in each patch, which is reasonable. Even at magnitudes brighter than the limiting magnitudes, the completeness is only 80-95% and does not reach 100% because of the hiding effect by bright objects. For NB527, the completeness is relatively lower than that of other bands. This seems to be caused by the source confusion due to the very deep depth with a moderate PSF size. Indeed, the surface number density of the detected objects in NB527 up to a 5σ limit of 26.46 AB (φ2.0 ′′ ) corresponds to an areal fraction of ∼ 13% even if we assume a circular area with a radius equal to the PSF FWHM size (r = 0.82 ′′ ) per object. For NB973, we also present the average completeness over all the patches estimated similarly but independently by ; it is slightly lower than ours but is broadly consistent. Actually, Itoh et al. (2018) used hscPipe version 4.0.5, and their 5σ-limiting magnitude in the 1.5 ′′ -diameter aperture of 25.0 was 0.4 magnitude shallower than that of this study.

Number counts
For the verification of the CHORUS PDR1 catalog, we calculated the source number counts in the CHORUS NBs and compared the results with those reported in the literature. We adopted the undeblended convolvedflux magnitudes (0.84 ′′ FWHM convolution and φ1.5 ′′ aperture case, except for NB387, where 1.1 ′′ FWHM convolution and φ2.0 ′′ aperture case is adopted) as the total magnitude of the objects and excluded the objects and area in the masked regions. The total magnitudes for each object were corrected for the Galactic extinction using the values listed in the CHORUS PDR1 catalog which are estimated from Schlegel et al. (1998)  instance, ∼ 0.1 mag even in NB387 in which the extinction effect is the severest. We also applied a catalog flag, forced.merge peak XXXX is True, where XXXX = n387, n527, n718, i945, or n973, to ensure that we count only the sources detected in the NB that we are considering. Otherwise, the number count included the sources that were not actually detected in the NB because hscPipe lists all the sources detected in any single band in the catalogs. Because the BB images of the HSC SSP are deeper than the CHORUS NB images, numerous faint BB sources, which are not detected in any of the CHORUS NBs, are also listed in the CHORUS PDR1 catalog. Fig. 6 shows the number count measurements. We only present the data points for those magnitudes at which the completeness (see Fig. 5) is larger than 10%. We also display the literature data of both the BBs and NBs whose wavelengths are close to those of each NB. There are excellent agreements between all the CHORUS number counts and the literature data. For example, the number counts of NB718, NB973, and IB945 present excellent agreement with the data of Subaru/Suprime-Cam (i from Kashikawa et al. 2004 for NB718, and z from Kashikawa et al. 2004 andNB921 from Ouchi et al. 2010 for NB973 andIB945). The NB527 number count also agrees with those of Subaru/Suprim-Cam NB503 from Ouchi et al. (2008) and V from Kashikawa et al. (2004); however, there are slight differences from those in bluer bands, like B from Kashikawa et al. (2004) and NB497 from Yamada et al. (2012). Because the slopes of the number counts reported in the literature are steeper at shorter wavelengths, the CHORUS results seem to follow this trend. For NB387 after the zero-point correction (see sec. 3), the number count data agree with those of Suprime-Cam/NB387 from Nakajima et al. (2012), Suprime-Cam/B from Kashikawa et al. (2004), and CFHT/MegaCam/u from Sawicki et al. (2019) at magnitudes fainter than 22. At brighter magnitudes, the NB387 data of ours and Suprime-Cam lie between the two BB literature data, which is reasonable given the trend mentioned above, because the wavelength of NB387 is indeed between those of the Band u-bands. The zero-point uncorrected NB387 number count is less than those reported in the literature data. These results suggest the necessity and validity of the zero-point correction.

How to use the catalog and example use cases of the CHORUS data
The CHORUS PDR1 catalog is available at the Catalog Archive Server (CAS) 15 of the HSC-SSP. An example SQL query to obtain a list of the coordinates, patch, the Galactic extinction in NB387, φ2.0 ′′ aperture NB387 magnitudes, NB387 total magnitudes (undeblended convolvedflux with 1.1 ′′ FWHM convolution and φ2.0 ′′ aperture), and the corresponding aperture The CHORUS PDR1 (forced) catalog consists of 7 files of s18a chorus.forced, s18a chorus.forced2, ..., s18a chorus.forced7. The FITS files of the masked areas (Fig. 2), spatial maps of PSF sizes (Fig. 3) and limiting magnitudes (Fig. 4), numerical data of detection completeness simulations (Fig. 5), and numerical data of measured number counts (Fig. 6) are also available through the website 16 .
In the following, we describe four specific studies ongoing with the CHORUS data and other ancillary studies including results published so far.

Lyman continuum survey of galaxies and AGNs
In the CHORUS project, we are performing a LyC survey for galaxies and AGNs at z ≃ 3.3 and 4.9 using NB387 and NB527, respectively (see Fig. 1). The samples are mainly the LAEs selected by NB527 (z ≃ 3.3) and NB718 (z ≃ 4.9). An average LyC-to-UV flux density ratio, or fesc, of the LAEs will be estimated by stacking analysis. Combining the LyC-to-UV ratios corrected for the IGM transmission and the observed UV luminosity densities, we can estimate the LyC luminosity densities and examine whether galaxies dominate the LyC emissivity at z ∼ 3-5. We can also detect the LyC from individual LAEs and AGNs. The number of such LyC galaxies at z > 3 is still limited (e.g., Iwata et al. 2009;Vanzella et al. 2012;Steidel et al. 2018;Iwata et al. 2019). This program will increase that number significantly. Follow-up spectroscopy of these reliable LyC galaxies allows an fesc calibration from the non-ionizing properties, such as a high [O III]/[O II] ratio de Barros et al. 2016;Nakajima et al. 2020). Once we establish the relation between the LyC emissivity (or fesc) values and the properties of galaxies observable at longer wavelengths at z < ∼ 5, we can estimate fesc of the galaxies in the epoch of reionization (z > 6). In this epoch, direct LyC observations are impossible owing to the severe absorption caused by the neutral hydrogen atoms in the IGM (Inoue & Iwata 2008;Inoue et al. 2014).

Surveys of faint AGNs as dual-emitters
Using NB718 and NB921, we can select a population of faint AGNs at z ≃ 4.9 as Lyα-C IV dual emitters if we appropriately treat the effects of broad lines whose line widths are broader than the NB widths (Fig. 1). This selection method is unique because these faint AGNs with a magnitude range similar to that discussed by Giallongo et al. (2015) are difficult to be isolated from the LBGs by any combination of BB colors. Moreover, our NB527 data allows the direct measurement of the LyC of these faint AGNs. The faint-AGN population reported by Giallongo et al. (2015) can reproduce the LyC emissivity required to keep the Universe ionized at z ∼ 5-6 if fesc = 1 (Madau & Haardt 2015). We examine if this is the case not only by measuring the faint-AGN LF but also by measuring the fesc of this population. Simultaneously, we measure fesc of galaxies at the same redshift and study their contribution directly. Subsequently, we identify the population that dominates the cosmic LyC emissivity at z ∼ 3-5 and infer the source of the cosmic reionization at z > 6.

Surveys of Pop III-dominated galaxies
Pop III stars are expected to have high LyC emissivity and may play a significant role in the early phase of reionization (e.g., Sokasian et al. 2004). However, observational evidence for this stellar population is still inconclusive (e.g., Kashikawa et al. 2012;Rydberg et al. 2015;Vanzella et al. 2020), and its characteristics, e.g., the initial mass function and the LyC emissivity, are unknown. Theoretically, Pop III star clusters can be formed even at z ∼ 2-3, depending on the efficiency of the IGM metalenrichment (e.g., Tornatore et al. 2007;Johnson 2010;de Souza et al. 2011). Indeed, pristine (Z < 10 −4 ) gas clouds, from which Pop III stars may form, are found at z ≃ 3 (e.g., Fumagalli et al. 2011). Identifying the epoch of the termination of a metal-free star formation is also a very important problem. The He II emission line can be an indicator of Pop III stars (e.g., Tumlinson et al. 2001;Schaerer 2003). Using our well-tuned set of NBs, we can search for Lyα-He II dual emitters at z ≃ 2.2, 3.3, and 4.9 (see Table 1 and Fig. 1). Such a dual emitter search was performed by Nagao et al. (2008) with the Suprime-Cam for z ∼ 4. We update the upper limit on the He II line luminosity density or the Pop III star formation rate density by Nagao et al. (2008) with a 7 times larger area and up to ∼ 1 mag deeper imaging data at a wider redshift range of z = 2-5. This will enable us to constrain the metal-enrichment efficiency in the IGM stronger than before.

Mapping spatial distributions of ionized bubbles
Because Lyα photons are sensitive to xHI, the distribution of the LAEs depends on both xHI and the large-scale structure (LSS). However, LBGs are selected by the continuum and are not very sensitive to xHI but simply trace the LSS. If we normalize the number density of the LAEs for each location by that of the LBGs, we can isolate the spatial distribution of xHI; we can visualize the xHI map directly from the LAE-to-LBG density ratio. An essential point is to select LBGs with a similar redshift to the NB-selected LAEs whose redshift variation is as small as ∆z ≃ 0.1. Although such an LBG selection is not feasible by BB filters, we can reach ∆z ≃ 0.3 using an intermediate filter, IB945 for redshifts around z = 6.6 and 7.0. From the spatial variation of the observed number density ratio, we can visualize an xHI map far ahead of the future 21 cm experiments, such as the SKA. Such an xHI mapping opens a new pathway to examine the ionization topology, for example, allowing the measurements of the ionized bubble size directly.

Ancillary sciences
Very wide (1.6 deg 2 ) and very deep ( > ∼ 0.1L * ) LAE samples at z ≃ 2.2, 3.3, 4.9, 6.8, and 7.0 are obtained by this program. An LAE LF at z = 7.0 has already been published as the CHORUS paper II ) that reports a very similar LF shape to those at z = 5.7 and 6.6 with a small overall decrement compared to the lower-z ones. This indicates that the reionization process does not alter the shape of the LAE LF in contrast to a previous claim of a bright-end hump (Zheng et al. 2017). A comprehensive LAE search in the CHORUS data in conjunction with the HSC SSP data is ongoing, adopting a machine learning technique (Ono, Y. et al. in preparation). A study of intensity mapping of Lyα emission is also ongoing (Kikuchihara, S. et al. in preparation). More LAE related studies including environmental ones and optical emission line properties in the near-infrared band are also ongoing. We are also searching for spatially extended LAEs -Lyα blobs -in the imaging data. A result of a Lyα blob survey at 4.9 ≤ z ≤ 7.0 has been published recently as the CHORUS paper III (Zhang et al. 2020) that reports a Lyα blob at z ≃ 7.0, the highest-z blob found so far. [O II], [O III], Hβ, and Hα emitters at various redshifts are also studied with the CHORUS data in conjunction with the HSC SSP data (Hayashi et al. 2020) and Spitzer data . More line emitter studies including AGN surveys will be performed in future. The CHORUS PDR1 data will be useful for various studies performed all over the world. We strongly encourage the readers' own studies using the dataset.

Summary
In this study, we conducted a deep imaging survey using a set of narrow-band (NB) and intermediate-band (IB) filters, named as Cosmic HydrOgen Reionization Unveiled with Subaru (CHORUS), with the Hyper Suprime-Cam (HSC) equipped on the Subaru Telescope. The filter set used in the survey is shown in Figure 1. The wavelengths, NB387, NB527, NB718, NB973, and IB945, summarized in Table 1, were chosen in a very wellorganized manner to maximize the scientific outcomes when the filters were used in combination with another NB, NB921, from the Subaru Strategic Program for HSC (HSC SSP). The observations were performed on 18 nights during the period from January 2017 to December 2018 in an intensive program of the Subaru Telescope. The observing field was the COSMOS field (Scoville et al. 2007). The data reduction, source detection, and multiband photometry were conducted using hscPipe, the official pipeline software for the HSC developed by the HSC SSP team (Bosch et al. 2018). The hscPipe version used was 6.7, the same as that used for HSC SSP public data release (PDR) 2 (Aihara et al. 2019). Because there is a known problem for the cmodel magnitudes in this version, we recommended using the undeblended convolvedflux magnitudes instead. The imaging data qualities were examined thoroughly in a spatially dependent manner, as shown in Figures 3 and 4. For each NB/IB, typical values of the size of point spread functions, 5σ limiting magnitude, and survey area excluding the carefully defined masked regions (see Figure 2 as an example) are summarized in Table 2. An extensive set of mock observation simulations was performed to estimate the source detection completeness, as shown in Figure 5. We checked the source number counts and compared them with previous studies in Figure 6. An excellent agreement with the literature data was obtained in each NB/IB. All the images and photometric catalogs are publicly available through the HSC SSP database. This data release is called CHORUS PDR 1.