The Quasar Luminosity Function at Redshift 4 with Hyper Suprime-Cam Wide Survey

TL;DR

This study measures the quasar luminosity function at redshift 4 using Hyper Suprime-Cam data, revealing a flat faint-end slope and no significant evolution from lower redshifts, impacting our understanding of quasar populations and their contribution to cosmic UV light.

Contribution

It provides the first wide-area measurement of the z=4 quasar luminosity function with a large sample, showing a flat faint-end slope and no evolution from z~2, and compares UV and X-ray selected AGN populations.

Findings

The luminosity function is well described by a double power-law with a flat faint-end slope.

No clear evidence of redshift evolution in the knee and faint-end slope from z~2 to 4.

UV-selected quasars underrepresent obscured AGNs compared to X-ray luminosity functions.

Abstract

We present the luminosity function of z=4 quasars based on the Hyper Suprime-Cam Subaru Strategic Program Wide layer imaging data in the g, r, i, z, and y bands covering 339.8 deg^2. From stellar objects, 1666 z~4 quasar candidates are selected by the g-dropout selection down to i=24.0 mag. Their photometric redshifts cover the redshift range between 3.6 and 4.3 with an average of 3.9. In combination with the quasar sample from the Sloan Digital Sky Survey in the same redshift range, the quasar luminosity function covering the wide luminosity range of M1450=-22 to -29 mag is constructed. It is well described by a double power-law model with a knee at M1450=-25.36+-0.13 mag and a flat faint-end slope with a power-law index of -1.30+-0.05. The knee and faint-end slope show no clear evidence of redshift evolution from those at z~2. The flat slope implies that the UV luminosity density of the quasar population is dominated by the quasars around the knee, and does not support the steeper faint-end slope at higher redshifts reported at z>5. If we convert the M1450 luminosity function to the hard X-ray 2-10keV luminosity function using the relation between UV and X-ray luminosity of quasars and its scatter, the number density of UV-selected quasars matches well with that of the X-ray-selected AGNs above the knee of the luminosity function. Below the knee, the UV-selected quasars show a deficiency compared to the hard X-ray luminosity function. The deficiency can be explained by the lack of obscured AGNs among the UV-selected quasars.

Tables

Table 1. Table 1: Sub-regions defined in this paper

Sub-regions	Field name	Central coordinate	Galactic coordinate	Area	Eff.Area††{\dagger}††{\dagger}footnotemark: $†$	Eff.Area‡‡{\ddagger}‡‡{\ddagger}footnotemark: $‡$	$N_{\mathrm{org}}$	$N_{\mathrm{sample}}$
		(RA,Dec)	($l$, $b$)	(deg2)	(deg2)	(deg2)
WideA	XMM-LSS	(\timeform34,\timeform-4)	(\timeform168,\timeform-59)	62.8	41.6	35.1	594	372
WideB	GAMA09H	(\timeform135,\timeform0)	(\timeform229,\timeform28)	81.1	45.2	38.1	747	341
WideC	WIDE12H	(\timeform180,\timeform0)	(\timeform276,\timeform60)	107.9	68.5	58.6	984	527
	GAMA15H	(\timeform218,\timeform0)	(\timeform349,\timeform54)
WideD	HECTOMAP	(\timeform243,\timeform43)	(\timeform68,\timeform47)	88.0	46.4	40.2	902	428
	VVDS	(\timeform340,\timeform0)	(\timeform68,\timeform-48)
	AEGIS	(\timeform215,\timeform52)	(\timeform95,\timeform60)
Total				339.8	201.7	172.0	3,227	1,668∗*∗*footnotemark: $\ast$

Table 2. Table 2: Database selection criteria

Flag	condition
detect_is_primary	True
flags_pixel_edge	not True
flags_pixel_saturated_center	not True
flags_pixel_cr_center	not True
flags_pixel_bad	not True

Table 3. Table 3: Number count of the z = 4 𝑧 4 z=4 quasars

$i$	$N_{\mathrm{obs}}$	$n$	$n_{\mathrm{corr}}$
(mag)		(deg-2 mag-1)	(deg-2 mag-1)
$20.00-20.50$	34	$5.50\times {10}^{-1}$	$5.48\times {10}^{-1}$
$20.50-21.00$	78	$1.30\times {10}^0$	$1.30\times {10}^0$
$21.00-21.50$	121	$1.97\times {10}^0$	$1.97\times {10}^0$
$21.50-22.00$	162	$2.70\times {10}^0$	$2.69\times {10}^0$
$22.00-22.50$	170	$3.35\times {10}^0$	$3.31\times {10}^0$
$22.50-23.00$	189	$4.17\times {10}^0$	$3.88\times {10}^0$
$23.00-23.50$	258	$6.65\times {10}^0$	$4.62\times {10}^0$
$23.50-24.00$	656	$2.54\times {10}^1$	$1.03\times {10}^1$

Table 4. Table 4: Binned z = 4 𝑧 4 z=4 quasar luminosity function

$M_{1450}$	$N_{\mathrm{corr}}$	$\log \left(\Phi \right)$††{\dagger}††{\dagger}footnotemark: $†$	$\sigma$∗*∗*footnotemark: $\ast$
(mag)
HSC S16A-Wide2
$-21.875$	47.7	$-6.253$	80.733
$-22.125$	129.0	$-6.088$	71.892
$-22.375$	103.7	$-6.219$	59.367
$-22.625$	92.8	$-6.298$	52.208
$-22.875$	78.5	$-6.382$	46.855
$-23.125$	95.8	$-6.297$	51.538
$-23.375$	81.4	$-6.369$	47.366
$-23.625$	90.4	$-6.341$	47.960
$-23.875$	84.7	$-6.405$	42.783
$-24.125$	74.0	$-6.493$	37.343
$-24.375$	76.0	$-6.487$	37.385
$-24.625$	63.0	$-6.577$	33.384
$-24.875$	43.0	$-6.745$	27.422
$-25.125$	42.0	$-6.756$	27.088
$-25.375$	30.0	$-6.899$	23.032
$-25.625$	14.0	$-7.189$	17.287
$-25.875$	8.0	$-7.077$	29.581
SDSS DR7 $z=3.6-4.2$
$-26.125$	260.0	$-7.387$	2.545
$-26.375$	287.0	$-7.552$	1.657
$-26.625$	234.0	$-7.642$	1.490
$-26.875$	144.0	$-7.853$	1.169
$-27.125$	111.0	$-7.966$	1.026
$-27.375$	66.0	$-8.192$	0.791
$-27.625$	30.0	$-8.534$	0.534
$-27.875$	17.0	$-8.781$	0.402
$-28.125$	9.0	$-9.057$	0.292
$-28.375$	2.0	$-9.710$	0.138
$-28.625$	3.0	$-9.534$	0.169
$-28.875$	2.0	$-9.710$	0.138

Table 5. Table 5: Parameters of the best fit z = 4 𝑧 4 z=4 quasar luminosity function with a double power-law model

Method	$\log \left({\Phi }^{\ast }\right)$††{\dagger}††{\dagger}footnotemark: $†$	$\alpha$	$\beta$	$M_{1450}^{\ast }$	Note
		[Faint End]	[Bright End]	(mag)
Maximum Likelihood	$2.66\pm 0.05$	$-1.30\pm 0.05$	$-3.11\pm 0.07$	$-25.36\pm 0.13$	Fitting to the sample
${\chi }^2$ minimization	$2.41\pm 0.56$	$-1.32\pm 0.08$	$-3.18\pm 0.11$	$-25.44\pm 0.19$	Fitting to the binned luminosity function

Equations

$$ishape\_hsm\_moment\_11$$

$$ishape\_hsm\_moment\_22$$

$$0.65(g-r)-0.30$$

$$3.50(g-r)-2.90$$

$$(g-r)$$

$$-2.25(i-z)+0.400$$

$$(i-z)$$

$$g_{\mathrm{HSC}}$$

$$r_{\mathrm{HSC}}$$

$$i_{\mathrm{HSC}}$$

$$z_{\mathrm{HSC}}$$

$$y_{\mathrm{HSC}}$$

$$\frac{\mathrm{d}\Phi }{\mathrm{d}{M}_{1450}}(M_{1450})=\frac{{\sum }_i{V}_a(M_{1450,i}{)}^{-1}}{\Delta M_{1450}}$$

$$V_a(M_{1450})=\int d_A(z{)}^2(1+z{)}^{3}c\frac{\mathrm{d}\tau }{\mathrm{d}z}(z)A(M_{1450},z)\mathrm{d}z,$$

$$\varphi (M_{1450},z)$$

$$=$$

$$\mathcal{L}=$$

$$N(M_{1450},z)=\frac{\mathrm{d}{\Phi }^{\mathrm{model}}}{\mathrm{d}{M}_{1450}}{d}_A(z{)}^2(1+z{)}^{3}c\frac{\mathrm{d}\tau }{\mathrm{d}z}(z)A(M_{1450},z).$$

$$\log (l_{2keV})=0.721\log (l_{2500\mathring{\mathrm{A}}})+\mathrm{4.531.}$$

Full text

1 Introduction

After the discovery that every massive galaxy harbors a super massive black hole (SMBH) in its center, the issue of how these SMBHs formed and evolved over cosmic history became one of the major unanswered questions in observational cosmology. The cosmological evolution of the AGN luminosity function, which reflects the growth history of SMBHs through accretion, has been intensively investigated using large AGN samples (e.g, [Richards et al. (2006)]; [Croom et al. (2009)]; [McGreer et al. (2013)]; [Ross et al. (2013)]; [Ueda et al. (2014)]). There is an overall trend that the number density of more luminous AGNs peaks at higher redshift, so-called "down-sizing" and "antihierarchical" growth, and the number density of the most luminous AGNs, i.e. quasars, shows rapid decline from the peak redshift, $z\sim 3$ to the local universe.

On the other hand, the number density of less-luminous AGNs beyond its peak redshift would show a hint of milder decline (Ueda et al., 2014). Such evolution would be the first evidence of a different evolutionary trend in the early universe that less-luminous AGNs are more numerous at first, then luminous AGNs grow later, i.e. "up-sizing" and "hierarchical" growth of SMBHs. In order to disclose the discussion of the evolutionary trend seen in the early universe, it is important to determine the shape of the AGN luminosity function, especially around its knee, at each redshift, and examine the redshift evolution of the shape.

The faint-end slope of the AGN luminosity function in the early universe is also important to evaluate the contribution of AGNs to the UV ionizing photon budget. Recently utilizing X-ray-selected AGNs, Giallongo et al. (2015) derived $z=4,5,6$ AGN luminosity functions with a steep slope and high number density in the faint-end. If we take the high number density of less-luminous AGNs at face value, the AGN emissivity of UV ionizing photons would be as high as the value required to keep the intergalactic medium highly ionized, and can significantly contribute to the cosmic reionization.

The strategic survey program of the Subaru telescope with the Hyper Suprime-Cam (HSC-SSP; Aihara et al. (2017a)), which started in Mar. 2014 and is assigned 300 nights for 5 years, provides us a unique opportunity to determine the shape of the quasar luminosity function in the early universe with unprecedented accuracy. An overview of the camera and details of the dewar system are given in Miyazaki et al. (2017) and Komiyama et al. (2017), respectively. The Wide-layer component of the survey covers 1,400 deg2 in the equatorial region down to 26.8, 26.4, 26.4, 25.5, and 24.7 mags in the $g$,$r$,$i$,$z$,$y$ bands, respectively, with 5$\sigma$ for point sources in the 5 year survey. The latest internal release of the data, S16A-Wide2, covers 339.8 deg2, including the edge regions where the final depth have not been achieved (Aihara et al., 2017b). The depths of the Wide-layer component are about 1 mag deeper than previous wide-field surveys covering similar area (e.g. Canada-France-Hawaii Telescope Legacy Survey wide fields covering 150 deg2). Furthermore, the image quality is better than other wide-field surveys; the median seeing size in the $i$ band is \timeform0.61" (Aihara et al., 2017b). Thanks to the depth and image quality of the data, we can select candidates of $z=4$ quasars with stellarity and $g$-dropout selections reliably down to $i=24.0$ mag, which corresponds to $M_{\rm 1450}=-22$ mag, i.e. well below the knee of the quasar luminosity function. In combination with the spectroscopically identified $z=4$ quasars from the Sloan Digital Sky Survey (SDSS) data release 7 (DR7) (Schneider et al., 2010), we can construct the $z=4$ quasar luminosity function covering a wide-luminosity range.

In this paper, we select candidates of $z\sim 4$ quasars applying a $g$-band dropout selection to objects with stellar morphology. The selection criteria and definition of the sample are described in section 2. We evaluate the effective survey area by constructing $z\sim 4$ quasar photometric models based on a SED library of quasars and a noise model of the HSC Wide-layer dataset. The statistical contamination rate of the sample is estimated with photometric data of Galactic stars and galaxies. The effective survey area and the contamination rate are discussed in section 3. Because most of the selected candidates do not have spectroscopic information, we derive their photometric redshifts with a Bayesian method using the library of quasar photometric models. In order to construct the $z=4$ quasar luminosity function covering a wide luminosity range, we combine our results with a sample of luminous quasars in the same redshift range utilizing the SDSS DR7 quasar catalog. The properties of the sample and the derivation of the $z=4$ quasar luminosity function is described in section 4. In section 5, we discuss the evolution of the quasar luminosity function at $z>2$, and we compare the $z=4$ quasar luminosity function with the luminosity function of X-ray selected AGNs at $z=4$. Throughout the paper, we use cosmological parameters of $H_{0}=70$ km s*-1* Mpc*-1*, $\Omega_{m}=0.3$, and $\Omega_{\lambda}=0.7$. All magnitudes are described in the AB magnitude system.

2 Sample definition

2.1 HSC-SSP Wide-layer catalog

Candidates of $z\sim 4$ quasars have been selected from the Wide-layer component of the HSC-SSP database. The S16A-Wide2 internal release of the survey covers an area of 339.8 deg2, where the $g$, $r$, $i$, $z$, and, $y$ bands data are available. The area includes edge regions where the integration time has not reached the target value and the final depth has not been achieved, therefore the area is larger than the full-color and full-depth area. The data are reduced with hscPipe-4.0.2 (Bosch et al., 2017).

In the S16A-Wide2 release, the Wide-layer survey data is separated into 7 continuous fields. Rough central coordinates of the fields are summarized in table 1. We divide the fields into 4 sub-regions of WideA-D. The total area of each sub-region is listed in table 1. The total area only considers area covered by all 5 bands. We remove some patches where the color sequence of stars brighter than $i_{\rm PSF}<22$ mag shows significant offset from that expected from the Gunn-Stryker stellar spectro-photometric library (Gunn & Stryker, 1983); patches which have offset in the patch_qa table larger than 0.075 mag either in the $g-r$ vs. $r-i$, $r-i$ vs. $i-z$, or $i-z$ vs. $z-y$ color-color plane (see section 5.8.4 in Aihara et al. (2017b)). We also remove tract 8284, because the photometric data in the region are unreliable.

We first select objects which are not flagged with saturation, bad pixel, and cosmic-ray hit and located at the edge of the detector. In the HSC pipeline, crowded objects are deblended by a deblending process. We consider the objects after the deblending process. The selection conditions are summarized in table 2.

We use PSF magnitudes for stellar objects. The PSF magnitude is determined by fitting a model PSF at each position to an image of an object. For discussions on extended objects, we refer their CMODEL magnitudes. CMODEL magnitudes are determined by at first fitting exponential and de Vaucouleurs profiles separately, then fitting them together (Bosch et al., 2017). Both profiles are convolved with the model PSF at the position of each object. The profiles of stellar objects are fitted with exponential and/or de Vaucouleurs profiles with radius of 0, CMODEL magnitudes should be consistent with the PSF magnitudes. However, PSF magnitudes are more stable, because the fitting is not affected by the uncertainty in the profile model. All the magnitudes are corrected for Galactic extinction based on dust maps by Schlegel et al. (1998).

2.2 Selecting stellar objects

We consider AGNs whose morphology is dominated by their nuclear stellar component as quasars. Stellar objects are selected by comparing the second order adaptive moments of an object with those of the PSF at the position of the object. We employ the adaptive moments measured in the $i$ band, because $i$ band images of the Wide-layer survey are selectively taken under good seeing conditions for weak-gravitational lensing study. The median seeing size of the $i$ band images is \timeform0.61". The adaptive moment evaluated with the algorithm described in Hirata & Seljak (2003) is available in the HSC-SSP database. We use the conditions,

[TABLE]

to select stellar objects. We remove objects whose adaptive moments are not measured correctly and listed as "nan". We do not refer CLASSIFICATION_EXTENDEDNESS, which is provided in the database, because the classification criterion is more optimized to select extended objects, and contamination of extended objects among stellar object is not negligible compared to the selection criteria applied in this paper.

The incompleteness and contamination of the above criteria are evaluated with the simulated images of the Wide-layer dataset in the Cosmic Evolution Survey (COSMOS) region, where $i$-band imaging data with Advanced Camera for Surveys (ACS) on Hubble Space Telescope (HST) is available. The simulated images are constructed with selected images from ultra-deep HSC-SSP survey in the COSMOS region. Three stacked images simulating good, median, and bad seeing conditions during the Wide-layer observations are provided in the database as the COSMOS wide-depth stacks, they have FWHMs of \timeform0.5", \timeform0.7", and \timeform1.0". It should be noted that the $i$ band images are selectively taken under good seeing conditions, therefore most of the data are taken under the good and median seeing conditions, and rarely taken under the bad seeing condition (see figure 3 of Aihara et al. (2017b)). The ACS image is deeper than the Wide-layer survey, and stellarity classification is available in the object catalog (Leauthaud et al., 2007). The classification is more robust than that based on the HSC images thanks to the sharper PSF with HST. The incompleteness is defined by the fraction of ACS stellar objects mis-classified as extended objects with the above criteria on the HSC catalog. The contamination is defined by the fraction of the ACS extended objects that are classified as stellar with the above criteria among the entire HSC stellar objects.

Figure 1 shows the resulting incompleteness and contamination of the classification. In the median condition, the completeness is more than 80% down to $i=23$ mag and decreases to 40% at $i=24$ mag, while the fraction of contaminating objects is less than 5% down to $i=23$ mag, and increases to 30% at $i=24$ mag. In the good and bad seeing conditions, the completeness and contamination vary accordingly. As we explained above, most of the $i$ band data are taken under the condition simulated as the good and median conditions, therefore hereafter we refer the completeness and contamination evaluated with the median condition. We need to note that the current survey area includes regions with shallower depth, and in such regions, the above incompleteness and contamination rate may not be applicable. Because the contamination rate increases rapidly at $i>24$ mag, hereafter we only consider objects brighter than $i=24$ mag.

2.3 Color selection criteria for $z\sim 4$ quasars

Candidates of $z\sim 4$ quasars are selected from stellar objects on the $g-r$ vs. $r-z$ color-color diagram. Figure 2 summarizes the distribution of the S16A-Wide2 stellar objects with known spectroscopic information on the color-color diagram. The selection criteria are determined aimed at including as many quasars between $3.5<z<4.0$ as possible, while minimizing the contamination by other objects. The determined selection criteria are shown with solid lines. They are

[TABLE]

We use the third selection criterion to limit the redshift range of the sample to $z\sim 4.5$.

In order to constrain the quasar luminosity function at $z\sim 4$ without conducting further spectroscopic follow-up observation, we set relatively tight color selection criteria to minimize contamination of red Galactic stars. For example, in Ikeda et al. (2011), they use wider selection criteria to select $z\sim 4$ quasar candidates for spectroscopic follow-up observations in the COSMOS region, and they found a significant number of contamination by Galactic stars. The HSC colors of their Galactic stars are shown with open red circles in the left panel of the figure. We set the tighter selection criteria not to include these stars by loosing a non-negligible fraction of $z\sim 4$ quasars with known spectroscopic redshift. The fraction of quasars missed in the color selection is accounted for in the statistical evaluation of the survey effective area discussed in the next section.

In addition to the criteria on the $g-r$ vs $r-z$ plane, we also apply criteria on the $i-z$ vs $z-y$ plane with

[TABLE]

These criteria are necessary to further remove contamination by red Galactic stars, and to remove some outliers with unreliable photometry. The distribution of the spectroscopically identified stellar objects that meet the $g-r$ vs. $r-z$ color selection criteria in the $i-z$ vs. $z-y$ plane is shown in figure 3. The color selection in the $i-z$ vs. $z-y$ plane removes reddened quasars. In addition to the above color selection criteria, in order not to be affected by objects with low signal-to-noise ratio, we only consider objects with magnitude error less than 0.1 mag in both of the $r$ and $i$ bands.

In the left panel of figure 4, we plot the color and magnitude distributions of the 534 quasars at $z=3.0-4.5$ from the spectroscopic database of the twelfth data release (DR12) of the SDSS (Alam et al., 2015) in the S16A-Wide2 coverage. They are within the coverage and neither masked nor flagged by the masking process described in sections 2.1 and 2.4, except for the HSC $i<22$ mask. Red open and blue filled circles represent quasars which are recovered and missed in the above selection, respectively, with the photometric data of the HSC. Because the bright end of the sample can be affected by saturation or non-linearity, figure 4 is plotted with the colors in the SDSS photometry converted to the HSC system with the conversion discussed in section 3.3. Quasars with red $g-r$ color should be selected with the color selection, but they are missed above around $i<19.8$ mag. Although in this plot, we remove objects which have saturated pixel at the central $3\times 3$ pixels by the flagging discussed in section 2.1, still some bright $i<20.0$ mag stellar objects seem to be affected by saturation or non-linearity effects. Therefore, hereafter we limit the sample fainter than $>20.0$ mag in all of the $r$, $i$, $z$, and $y$ bands, for statistical discussions. There are 379 objects fainter than $i>20.0$ mag among the 534 objects.

In the right panel of the figure, we plot $g-r$ color as a function of spectroscopic redshift of the 379 SDSS quasars. The $g-r$ color of the quasars at $z>3.5$ are red, and most of them are selected by the above color selection criteria; 92 quasars are above $z=3.5$ and 61 are selected by the above $z=4$ quasar selections, i.e. 66% of the SDSS quasars are selected. Among the 31 missed quasars, 1, 26, and 4 objects are missed by the first stellarity, the second $g-r$ vs. $r-z$ color-color, and the third $i-z$ vs. $z-y$ color-color selections, respectively.

After applying the stellarity and color criteria to the $20.0<i<24.0$ objects which meet the flag conditions in table 2, we select 3,227 candidates of $z\sim 4$ quasars. The number of candidates in each sub-region is summarized in the column $N_{\rm org}$ of table 1. However, if we check the 5-band images of the selected candidates, we find that junk objects contaminate the sample. Therefore, we apply further masking based on the mask image and the bright star catalogs as described below.

2.4 Additional masking and masking junk objects

In addition to the flags listed in the database as shown in table 2, we further check similar flags in the mask images provided in the stacked image database. The primary aim of this further masking process is to remove the junk objects. Additionally, we try to easily mimic the flag conditions shown in the table 2 when we construct mock random objects uniformly distributed in the survey region, which is necessary in the evaluation of the effective survey area (see section 3.4). We mask objects which have either of the MP_EDGE, MP_BAD, MP_SAT, MP_NODATA, or MP_NOT_DEBLENDED flag within the radius of \timeform5" and MP_CR in the central $3\times 3$ pixels. These masking is usually stricter than the flag in table 2, therefore by considering these masking in the process of the random mock objects, we expect all of the flagging and masking conditions are simulated.

The above flags and masking are not sufficient to remove objects with bad photometry in the catalog. For example, satellite tracks remaining in the stacked images can be cataloged as objects detected in one band. Such tracks themselves are not selected as a $z\sim 4$ quasar candidate, but affect the photometry of real objects close to them. Ghost images and faint halo around bright stars can also be cataloged as objects meeting the above flag conditions, and also affect the photometry around them.

At first, we remove objects around bright stars and galaxies. In the HSC catalog, objects around bright stars are flagged based on Naval Observatory Merged Astrometric Dataset (Zacharias et al., 2005). We remove objects that have MP_BRIGHT_OBJECT flag within the radius of \timeform5" in the mask image. The process should mimic flags_pixel_bright_object_any flag in the database. In addition to the flag, we further consider bright objects cataloged in Guide Star Catalog (GSC) version 2.3.2. We remove objects within $r_{\rm mask}[^{\prime\prime}]=150.0$ from stars brighter than 10 mag, and $r_{\rm mask}$[′′] $=20.0+17.1\times(13.5-{\rm mag})$ from stars down to 15 mag. Among the multiple magnitudes available for one object in the GSC, we refer the brightest magnitude of the object in the above calculation. The size of the additional masks are empirically determined by checking distribution of "junk" objects around bright stars.

We also remove objects close to satellite tracks and faint halos around bright stars, because photometries of such objects are often unreliable. In the current database, such satellite tracks and faint halos are detected as a widely extended object, and the pixels associated with them are flagged as MP_DETECTED in the mask image same as for other real astronomical objects. We check $10^{\prime\prime}\times 10^{\prime\prime}$ region around each cataloged object and if more than 60% of the pixels around a object are flagged as MP_DETECTED in either of the 5 bands, the object is removed. Additionally, we examine the connected pixels around the object with this flag and if the total number of the connected pixels is more than 30% of the $10^{\prime\prime}\times 10^{\prime\prime}$ pixels, we also mask the object. Furthermore, if the number of detected pixels in the $i$ band is 2.0$\times$ larger than that in either of the $r$ or $z$ band, or those in the $g$, $r$, $z$, $y$ bands are 2.0$\times$ larger than the $i$ band, we mask the object. The size and fraction of the mask are determined by checking objects around remaining satellite tracks and faint halos.

The deblending process of the HSC pipeline does not provide reliable photometry for faint objects in the outskirts of relatively bright galaxies ($i<21$). The above masks around very bright objects are not sufficient to remove such faint objects around bright galaxies. Therefore, we remove objects around stars and galaxies brighter than $i<22.0$ mag. The radius of the mask, $r_{\rm mask}$, is determined with the adaptive moment of the object with $r_{\rm mask}=1.8\times\sqrt{\rm i\_hsm\_moment\_11}$. In order to recover real candidates with $i<22.0$ mag, we only consider the mask for objects fainter than $i>22.0$ mag.

In the selection process, we apply all of the above additional masking processes after selecting 3,227 candidates of $z\sim 4$ quasars with the stellarity and color selections described above. Once we apply the above masking processes, there are 1,668 candidates left. The number in each sub-region is summarized in table 1. We check the selected candidates by eye, and confirm almost all of the junk objects and objects in the outskirt of bright objects are removed. We do not apply a junk object removal with the eye-ball check. The same masking processes are also applied to random mock objects described in section 3.4 to mimic the object detection process.

3 $z\sim 4$ quasar number counts

3.1 Modeling photometric uncertainty at each position

In order to derive the number counts of the $z\sim 4$ quasar sample, we evaluate the effective survey area of the S16A-Wide2 dataset as a function of magnitude and redshift. The effective survey area of the S16A-Wide2 dataset needs to be evaluated with considering the variation of the depth of the images due to the variation of seeing and atmospheric transparency conditions and the number of available images at each location. At first, we establish a method to evaluate uncertainty at each position based on the variance images of the stacked images. Then, we construct libraries of quasar photometric models, randomly locate the quasar models within the survey region, and add random photometric error evaluated with the variance value at the random position. In this process, we assume that the photometric errors are dominated by the background noise. Finally, we apply the same magnitude and color selection to the random quasar models and evaluate effective survey area by the ratio of recovered $z\sim 4$ quasar models to total input models at each redshift and magnitude.

The photometric uncertainties of real objects correlate well with the variance values of the object positions. The correlation shows the dependence on the size of the object. Therefore, we construct equations to calculate photometric uncertainties based on the model PSF size and the variance at each position. In this uncertainty model, it is assumed that photometric uncertainty is not dominated by the object Poisson noise. In reality the objects at the detection limits are faint enough to meet the assumption.

3.2 Constructing a library of quasar photometric models with templates

The quasar photometric models are constructed based on a library of spectral energy distributions (SEDs) of $z\sim 4$ quasars, which describes the distribution of the SEDs of population of broad-line quasars. We make a library of quasar SEDs considering the scatter of the power-law slope, the broad-line equivalent width (EW), and strength of absorption by the inter-galactic medium (IGM). We assume quasar SEDs depend on luminosity, but do not depend on redshift.

The library of the quasar SEDs is constructed in the same way as described in Niida et al. (2016). The underlying power-law component, $f_{\nu}\propto\nu^{\alpha}$ is modeled with three components covering different wavelength ranges. Below 1100 Å, we use $\alpha$ of $-1.76$ following Telfer et al. (2002). We apply $\alpha$ of $-0.46$ and $-1.58$ between 1100 Å and 5011 Å and above 5011 Å, respectively, following Vanden Berk et al. (2001). For the two power-law components below 5011 Å, we assume a scatter around the average $\alpha$ with standard deviation of 0.30 following Francis (1996).

The strength of emission lines are modeled in relative to the C\emissiontypeIV emission line following the emission line ratios tabulated in Vanden Berk et al. (2001). The strength of the C\emissiontypeIV emission line is modeled by considering the dependence of the EW on quasar luminosity, the so-called Baldwin effect. We do not consider the luminosity dependence of the emission line ratios. The average and standard deviation of the C\emissiontypeIV emission line EW are measured as a function of luminosity from quasar spectra of the Baryon Oscillation Spectroscopic Survey (BOSS) in the SDSS III. The measurement covers the luminosity range of $M_{\rm 1450}=-21.5$ mag to $-29.5$ mag with 1 magnitude bin (Niida et al., 2016). In addition to the isolated emission lines, the Fe\emissiontypeII multiplet emission features and Balmer continuum are added by using the template in Kawara et al. (1996).

Finally we apply the effect of the absorption by the IGM. We use the updated number density of Ly$\alpha$ absorption system in Inoue et al. (2014). We also consider the scatter of the number density by a Monte Carlo method described in Inoue & Iwata (2008).

We construct 1,000 quasar SEDs in each magnitude bin from $M_{\rm 1450}=-21.5$ mag to $-29.5$ mag. Therefore in total 8,000 quasar SEDs are constructed. Each spectrum is redshifted, and converted to observed flux. In figure 5, the distributions of the model quasars in the redshift ranges $z=2.5-3.0$ and $z=3.5-4.3$ are plotted with contours in the left and right panels, respectively. The contours enclose 30, 60, 90% of the models in the redshift range from top to bottom. We only consider model quasars with apparent magnitude between $i=20.0$ to $22.0$ mag. The distributions are compared to those of the SDSS DR12 quasars which are in the same magnitude range covered by the HSC S16A-Wide2.

The 60% enclosing contour of the SED library reproduces that of the SDSS $z=2.5-3.0$ quasars. However, the 90% enclosing contour of the SDSS $z=2.5-3.0$ quasars shows extended distribution toward red $r-z$ color. Dust reddening to the quasar spectrum can explain the extended distribution. The direction of the extension is consistent with the reddening vector with Small-Magellanic Cloud-like dust extinction; the red arrow in the figure show the reddening vector with $E(B-V)=0.04$ mag. Richards et al. (2003) report 6% of SDSS quasars show red color, which is explained with $E(B-V)$ larger than 0.04 mag. After correcting for the bias against reddened quasars, they estimate total contribution of such dust-reddened broad-line quasars in the total quasar population as 15%.

In the current SED library, we do not include such population of quasars with dust reddening, primarily because current quasar samples at higher redshifts do not cover such a population (e.g. Richards et al. (2006), Ross et al. (2013); McGreer et al. (2013)). In the higher redshift range, $z=3.5-4.3$, the $r-z$ color distribution of the SDSS quasars is well reproduced by the distribution of the quasar SED library. That means current sample of $z=3.5-4.3$ SDSS quasars does not include the reddened quasar population seen at $z<3$ or the fraction of such reddened quasar is really decreasing. In the current photometric sample, it is hard to tackle such reddened quasar population due to heavy contamination by Galactic stars, and we need future wide-field multi-wavelength surveys to pick up the population.

On the other hand, $g-r$ color distribution of the SDSS quasars is narrower than that of the model quasars. The $g-r$ color distribution in this redshift range is mostly determined by the strength and scatter of the IGM absorption. Currently, the number of $z=3.5-4.3$ SDSS quasars covered by HSC photometry is rather limited to disclose the range of the color distribution, therefore we use the color distribution of the model quasars as a baseline sample of the full quasar population in this paper. Hereafter, we refer to the library of the quasar models as the library of quasar photometric models with templates.

It should be noted that the quasar SED library is constructed over a wide luminosity range using the spectra of less-luminous quasars at lower redshifts, assuming that the spectral shapes of the quasars do depend on luminosity, but do not depend on redshift. Therefore, the library extends to fainter luminosity range than the SDSS sample at $z=4$ and is therefore applicable in this HSC study.

3.3 Quasar photometric models with a SDSS photometric sample

The other library of quasar photometric models is constructed by converting quasar PSF photometry from the SDSS DR12 quasar catalog into the HSC photometric system. Because the HSC survey area is not wide enough to cover large numbers of SDSS quasars at high-redshifts, and quasars brighter than $i<20.0$ mag are affected by saturation in the HSC photometry, we construct a library of model quasars by converting SDSS photometry to HSC photometry instead of directly using the HSC photometry of the SDSS quasars.

We derive the conversion by applying filter response curves of the HSC (Kawanomoto et al., 2017) and SDSS 111http://classic.sdss.org/dr7/instruments/imager system to the spectrophotometric library of Galactic stars (Gunn & Stryker, 1983). The conversion in each band is determined by a linear dependence on one color, and they are

[TABLE]

Although these equations are determined for the SEDs of Galactic stars, they are effective to quasars. If we compare the HSC photometry of SDSS quasars with the SDSS photometry converted with the above equations, converted SDSS colors of SDSS quasars are consistent with their colors measured in the HSC dataset within a scatter of $\sim 0.2$ mag without any systematic offset. It needs to be noted that the above equations are applicable only to smooth SEDs and may not work with SEDs with a break or strong emission line.

In figure 6, we compare the distribution of the model quasars with SDSS photometry at $z=2.5-3.0$ and $z=3.5-4.3$ with that of the SDSS quasars observed by the HSC system in the magnitude range $i=20-22$ mag. The distribution of the converted photometry is broadly consistent with the HSC colors of SDSS quasars directly measured in the HSC images. For quasars at $z=2.5-3.0$, the $g-r$ colors of the quasars with red $r-z$ colors show systematically bluer distribution. Hereafter, we refer to the library as the library of quasar photometric models with SDSS photometry.

3.4 Estimating survey area as a function of redshifts and magnitudes

The effective survey area is determined as a function of magnitude and redshift for the quasar sample. We use the above 2 libraries of quasar models to derive the effective survey area correcting the fraction of missed quasars among the "complete" sample. In order to derive the survey area at a certain magnitude limit, firstly random positions are selected within the survey region, then we apply the same masking processes described in section 2.4. We pick-up 6,000 random positions per deg2. The concept of the random positions is the same as what is available in the database as the random objects (Aihara et al., 2017b). Secondly we randomly assign one quasar model from one of the libraries. For the quasar models with template, we randomly draw a quasar model in the considered magnitude range, and for the quasar models with SDSS photometry, we randomly select one quasar, and normalize the SED to match the considered magnitude. Therefore, in the calculation with the quasar models with SDSS photometry, we ignore the luminosity dependence of the quasar SEDs. Then, we calculate uncertainties associated with the photometry in each band and we apply random fluctuations to the photometric model. Finally, we apply the magnitude ($20.0<i<24.0$) and color selection criteria to examine the fraction of recovered objects at the magnitude.

The resulting effective survey areas as a function of redshift at each magnitude are shown in figure 7. The left and right panels show the effective areas estimated with the quasar models with templates and those with SDSS photometry, respectively. They are broadly consistent with each other, though the area derived with the models with templates shows a broader redshift distribution than that with the models with SDSS photometry as expected from the color distributions shown in figure 5. Both areas show a break at $i=22.0$ mag, caused by the mask around objects with $i<22.0$ mag. The peak of the survey area is 133.0 deg2 for $i=21.4$ mag quasars at $z=3.65$. In table 1, effective areas after removing masked regions are summarized. The effective areas are calculated applying the same masking process described in section 2.4. For the 6-th and 7-th columns, the areas without and with the masks with objects brighter than $i<22$ mag, respectively, are shown. The total effective survey area for $i<22$ mag objects is 172.0 deg2. It should be noted that the effective area for $z\sim 4$ quasars includes not only the effect of the masked regions and shallow areas, but also the selection efficiency of the quasars at each redshift. The peak of the survey area for $z\sim 4$ quasar is 77% of the total effective survey area for $i<22$ mag objects. The fraction is broadly consistent with the fraction of $z>3.5$ SDSS quasars meeting the $z\sim 4$ quasar selection criteria (66%) as discussed in section 2.3.

The differential number counts of the $z\sim 4$ quasars are shown in figure 8 and summarized in table 3. The second and third columns show the raw number and surface number density of the quasars, respectively. We apply the average effective survey area between $z=3.6-3.8$, where the color selection efficiency is maximized. In the number counts, we correct for the incompleteness of the stellarity selection assuming the median condition in figure 1. We do not correct for the contamination by the extended objects, because we consider the contamination in the next sub-section. The number counts show steady increase toward the faint-end, and excess at magnitudes fainter than $i>23.0$ mag. We also examine the number counts in each sub-regions. The results are shown with open squares with error bars estimated with the Poisson statistics. The number counts are consistent with those derived in the entire region, and the effect of the cosmic variance seems to be small, thanks to the large survey area. The scatter at the brightest end is large due to the limited number of the sample.

3.5 Contamination by extended objects

The expected number of contamination by extended objects that meet the $z\sim 4$ quasar color selection criteria is evaluated. In addition to the Lyman Break Galaxies at $z\sim 4$, foreground $z<1$ galaxies can be contaminants to the $z\sim 4$ quasar sample, because they can have similar colors to $z\sim 4$ quasars due to their 4000 Å break. It should be noted that even with the ground-based HSC images, most of the $z\sim 4$ LBGs are extended, thanks to the good image quality in the $i$ band, thus only compact LBGs are contaminating the stellar object selection.

In figure 9, the distribution of the measured adaptive moment ratios in the S16A-Wide2 dataset of the 306 $z>3$ and $i>23$ galaxies that are spectroscopically identified in deep surveys. Broad-line AGNs are removed in this plot. The red line indicate the selection criteria of the stellar objects. Thanks to the good image quality of the $i$-band image of the Wide layer dataset, they are significantly extended compared to the stellar objects. The fraction of the galaxies that are classified as stellar is 6%. The fraction is larger than that observed among all galaxies with $i=24$ mag (2%), in this evaluation we use the distribution of all galaxies in each magnitude range, because it is possible that the spectroscopically identified sample can be biased toward compact galaxies, and we do not know the true size distribution of $z>3$ galaxies.

At first, we construct a catalog of extended objects that meet the color criteria from the HSC-SSP S16A-Wide2 database. Then, we apply the masking processes described above. The number counts of the object without flag is calculated with the same survey area used for the $z\sim 4$ quasar number count. The expected number counts of contamination by extended objects is calculated by multiplying the fraction of ACS extended objects classified as stellar in HSC criteria as a function of magnitude. The fraction is determined in the same way as described in section 2.2, but this time the fraction of HSC stellar objects among the ACS extended objects is calculated. We use the results with the median condition. In this calculation, we assume that the contaminating extended galaxies have the same size distribution as the entire galaxies in the same magnitude range. The resulting number of contamination by extended objects is shown with red open circles in figure 8. The estimated result suggests that the contamination by extended objects can contribute to the number counts below $i=23.5$ mag, and one third of the $z\sim 4$ quasar candidates in the $i=23.5-24.0$ mag bin can be contamination of extended objects, which are non-AGN galaxies. The estimation is consistent with the rapid increase of the contamination of extended objects shown in figure 1.

3.6 Contamination by Galactic stars

The number of contamination by Galactic stars is evaluated by multiplying the observed number counts of Galactic stars in each sub-region with the fraction of Galactic stars meeting the $z\sim 4$ quasar color selection criteria evaluated as a function of $i$ band magnitude. Both of the number counts and the fraction are estimated with a photometric library of Galactic stars in the HSC system.

A photometric library of Galactic stars is constructed as follows. At first, we collect a list of spectroscopically identified Galactic stars in the S16A-Wide2 database which are flagged and masked in the same way as described in sections 2.1 and 2.4. Then, we remove some outliers, which are outside of the stellar sequence seen in the $g-r$ vs. $r-z$ color-color diagram. The distribution of resulting stars on the color-color diagram is shown in figure 10 with red dots. This list of stars cover a wide range of stellar types. However, they can have different mixture of types from that at the HSC Wide-layer depth because most of the spectroscopic identifications are from SDSS DR12 database, whose depth is shallower. Therefore, we match the list of flagged and masked stellar objects in the S16A-Wide2 database with the above list of the spectroscopically identified Galactic stars on the $g-r$ vs. $r-z$ color-color diagram based on the distance on the diagram, $r_{\rm col}$, less than 0.1 mag. The distribution of resulting stars are shown in the figure with gray scale. It follows the distribution of the spectroscopically identified stars, but more heavily populated with late-type stars compared to that of the spectroscopically identified stars. Open green circles indicate the colors of stars in the spectro-photometric library of Gunn & Stryker (1983). They show systematic offset from the observed sequence of the late type stars redder than $r-z>1.0$. Difference in the stellar metallicity is thought to be the cause of the systematic offset (Fukugita et al., 2011).

Then, the number counts of Galactic stars are derived in the same way as for the $z\sim 4$ quasars. At first, we apply the masking process described in sections 2.1 and 2.4 to the list of stellar objects in the HSC database. Then the area of survey region is determined by random model objects constructed from a photometric library of Galactic stars. We assume that the photometric error in the photometric library is negligible, because we only consider stars which have similar colors to the bright spectroscopically identified stars with negligible photometric errors. Then, we add random photometric error associated at each location in the same way described in section 3.4 and pick up random object above the detection limit, i.e. $i<24.0$ mag and photometric errors in $i$ and $z$ bands are smaller than 0.1 mag. The resulting number counts are shown in figure 11. The raw number counts in each sub-region are shown in thin blue line. In order to derive the intrinsic number counts of Galactic stars, we correct for the incompleteness and contamination of the star-galaxy separation applying the fraction shown in figure 1. The resulting number count after the corrections are shown in thick red lines. The number count shows steady increase toward $i=22-24$ mag, but then shows a plateau or decline toward the faint-end.

The expected contamination rate to the quasar number counts in each sub-region is shown in figure 8 with the green solid line. Because photometric uncertainty increases toward fainter objects, the expected number of contamination of Galactic stars increases rapidly, although the number counts themselves show plateau at $i=22-24$ mag. The expected number of contamination varies in the 4 sub-regions depending on the distance from Galactic plane. The contamination by Galactic stars can contribute one third of the number counts in the $i=23.5-24.0$ and $i=23.0-23.5$ mag bins.

3.7 Number counts after correcting for the contamination

Based on the estimated number counts from the contamination, we evaluate the contamination rate, which is the fraction of contaminating objects among the candidates of $z=4$ quasars. For the contamination by Galactic stars, the expected number depends on the Galactic coordinate of a sub-region, we average the number after weighting the survey area of each sub-region. The resulting contamination rate in each $i$ band magnitude bin is shown in figure 12 with crosses. In the magnitude range fainter than $i>23.5$ mag, the contamination rate is estimated to be higher than 50%. In order to correct for the contamination statistically, we fit the contamination rate with the error function. The best-fit function is derived as $[p_{\rm cont}={\rm erfc}(-1.15(i-23.59))/2]$, and plotted as solid line in figure 12.

We correct for the expected contamination by assigning weight for each object based on the probability of non-contamination, i.e. the contamination rate subtracted from 1. If the contamination rate is 0.8 at the magnitude of an object, we assign weight of 0.2 for the object, and the weight is summed when we calculate the number count and the luminosity function. The total sum of the weight of 1,666 candidates is 1155.2, i.e. one third of the sample can be contamination mostly contributing $i>23.5$ mag. The number counts after correcting for the contamination are shown in figure 8 with blue filled circles and summarized in the fourth column of table 3. Once we correct for the contamination statistically, the number density shows monotonic increase toward the faint-end. The excess seen in the magnitude range $i>23$ mag can be explained with the contamination by Galactic stars and compact galaxies that meet the $g$-dropout selection.

4 Results

4.1 Redshifts and absolute magnitudes

Among the 1,668 $z=4$ quasar candidates, 76 of them have spectroscopic redshift information in the literature. Most of the redshift identifications come from the SDSS quasar surveys. A few of them are from deeper surveys (e.g. Akiyama et al. (2015)). Their redshift distribution is shown with a red histogram in figure 13. They distribute redshift range between 3.4 and 4.2. Two SDSS quasars that have redshift $z\sim 1$ are contaminating the sample. The two quasars have $g-r$ colors of 0.37 and 0.50 in the SDSS database, while those in the HSC are 0.81 and 0.95. Because they are located within \timeform34’ each other, their HSC photometry could be affected by an unknown photometric problem in the specific region.

For the remaining candidates, no spectroscopic redshift information is available. We derive their photometric redshifts using the library of quasar models with templates described in section 3.2. For each candidate, we calculate $\chi^{2}$ with all of the templates in the redshift range between $2.5<z<6.0$ with step of $dz=0.1$. As a Bayesian photometric redshift estimation, we take the average of redshift with weighting by $\exp(-\chi^{2}/2)$. We assume uniform prior in the above redshift range. The resulting photometric redshift is plotted against the spectroscopic redshift in figure 14. The derived photometric redshifts correlate well with the spectroscopic redshifts. If we remove 2 outliers whose redshifts are $z\sim 1$, the average and $\sigma$ of the difference between $z_{\rm phot}$ and $z_{\rm spec}$ are are $-0.007$ and $0.143$, respectively. The systematic offset between the $z_{\rm phot}$ and $z_{\rm spec}$ is negligible.

The distribution of the $z_{\rm phot}$ and $z_{\rm spec}$ is shown in figure 13 with green histogram. Most of the selected quasars distribute between $z=3.5$ and $4.3$. The range is consistent with the selection function evaluated in figure 7. The average of the redshift of the sample is 3.9. The redshift distribution of quasars in $21.0<i<23.5$, which are less affected by contamination, shows a similar distribution as the entire sample.

We derive $M_{\rm 1450}$ either with $z_{\rm phot}$ or $z_{\rm spec}$. $M_{\rm 1450}$ is derived by interpolating the broad-band photometry around rest-frame 1450 Å. The uncertainty of $M_{\rm 1450}$ associated with $z_{\rm phot}$ is evaluated by the difference of $M_{\rm 1450}$ with $z_{\rm spec}$ and $z_{\rm phot}$ for the 74 objects with $z_{\rm spec}$. The resulting $\sigma$ of the difference is 0.082 mag. The estimated $M_{\rm 1450}$ are plotted against redshift in figure 15 with filled red and open blue circles for objects with $z_{\rm spec}$ and $z_{\rm phot}$. Most of the spectroscopically identified objects are from the SDSS catalog, and distributed brighter than $-24$ mag.

4.2 Binned quasar luminosity function at $z\sim 4$

The luminosity function of quasars at $z\sim 4$ is calculated, using the photometric redshifts, $M_{\rm 1450}$ and the survey area as a function of redshift and $M_{\rm 1450}$. We evaluate the survey area again with the library of quasar models with templates. The resulting effective survey area as a function of redshift and $M_{\rm 1450}$ is shown with gray scale in figure 16. We also plot the effective survey area with contour in figure 15. The distribution of the sample in the redshift-$M_{\rm 1450}$ plane is consistent with that expected from the effective survey area. It should be noted that the effective survey area includes the efficiency of the quasar color selection, and the the quasar models does not contain reddened quasar SEDs (see section 3.2).

We estimate the binned luminosity function with $\sum 1/V_{a}$ estimator (e.g. Miyaji et al. (2000)) with

[TABLE]

where the sum is taken through objects in a $M_{\rm 1450}$ bin with width of $\Delta M_{\rm 1450}$. $V_{a}$ represents the available volume for the $i$-th object with $M_{{\rm 1450},i}$ in the bin, and is calculated as

[TABLE]

where $d_{A}(z)$ as the angular diameter distance and $\frac{{\rm d}\tau}{{\rm d}z}(z)$ as the look-back time per unit $z$. $A(M_{\rm 1450},z)$ represents the survey area shown in figure 16.

The resulting luminosity function is shown in figure 17 with red filled squares. Thanks to the wide and deep coverage of the HSC Wide-layer dataset, the faint-end of the quasar luminosity function at $z\sim 4$ is constrained with unprecedented accuracy, and the break of the luminosity function is clearly detected at around $M_{\rm 1450}\sim-25$ mag. The binned luminosity function is shown in table 4. The uncertainty, $\sigma$, in each bin is determined with the Poisson statistics.

4.3 $z=4$ quasar sample from SDSS DR7

In order to construct the luminosity function covering a wide luminosity range based on quasars selected with consistent selection criteria, we combine $z=3.6$ to $4.2$ quasars from the spectroscopically identified quasar catalog from SDSS DR7 (Schneider et al., 2010). We use the DR7 catalog, because the SDSS legacy survey with a uniform color selection for high-redshift quasars is fully covered in the catalog. In the SDSS legacy survey, candidates of high-redshift quasars are selected with stellar morphology and multi-color selection for spectroscopic follow-up (Richards et al., 2002). A uniform target selection is applied for spectroscopy, except for the early phase of the SDSS survey. In order to construct a statistical sample, we only consider objects which are spectroscopically observed as a science primary object based on the uniform target selection (SCIENCEPRIMARY$=1$ and TARGET selection flag in table 2 of Schneider et al. (2010)). The effective survey area of the SDSS DR7 sample is determined to be 6,248 deg2 by Shen & Kelly (2012).

The selection function of the uniform target selection is evaluated in Richards et al. (2006) as a function of redshift and magnitude. In the redshift range between $z=3.6$ and $4.2$, the SDSS target selection efficiency is estimated to be larger than 90%, mostly $\sim$100% (figure 6 of Richards et al. (2006)). The target selection for high-redshift quasar is complete down to $i<20.2$ mag, therefore we select 1,260 quasars brighter than $M_{\rm i}(z=2)<(z-3.0)-26.85$ as a statistically complete sample in the above effective area (figure 17 of Richards et al. (2006)). We assume that the sample is complete in the redshift range above the absolute magnitude limit. Utilizing the absolute magnitude limit as a function of redshift, we calculate effective survey volume for each object with their $M_{\rm i}(z=2)$. UV absolute magnitudes, $M_{\rm 1450}$, of the quasars in the sample are derived from the broad-band SDSS photometry of the quasars in the same way for the HSC $z=4$ quasars.

The resulting luminosity function is shown in figure 17 with blue filled circles and summarized in table 4. The second column indicates the number of the quasars in each absolute magnitude bin after correcting for the contamination rate. The luminosity function smoothly connects to the luminosity function determined with the HSC $z=4$ quasar candidates. Because the magnitude range of the HSC $z=4$ quasar candidates is fainter than $i>20.0$ mag and the SDSS DR7 sample is $i<20.2$ mag, there is no overlap in the bin of the luminosity function.

4.4 Double power-law model

The quasar luminosity function is generally well described by a double power-law form with

[TABLE]

where $M_{*}$ is the absolute magnitude of the knee and $\phi^{*}$ is the number density at that luminosity. $\alpha$ and $\beta$ are power-law slopes of the faint- and bright-ends, respectively. The double power-law model is fitted to the $z=4$ quasar samples from the HSC and SDSS DR7 with the maximum likelihood method (Marshall et al., 1983).

We use

[TABLE]

which is modified version of the maximum likelihood estimator for luminosity function (e.g. Miyaji et al. (2000)) $N(M_{\rm 1450},z)$ is the expected number of object with a model per unit absolute magnitude and redshift inverval,

[TABLE]

The normalization of the best-fit model is not constrained with the likelihood estimator, we determine the normalization such that the expected number of object from the model matches $N_{\rm obj}$. Uncertainties are evaluated by fixing a parameter to a value around its best-fit value and minimizing $\mathcal{L}$ with the other parameters. One sigma uncertainty of the parameter is determined by the range whose minimum $\mathcal{L}$ is larger by less than 1 from the best-fit $\mathcal{L}$. The uncertainty associated with the normalization is determined by the Poisson statistics. The best-fit parameters and associated uncertainties are summarized in the first line of table 5. The best-fit model is shown with red solid line in figure 17. The red solid line matches well with the binned luminosity function.

We also fit the binned luminosity function with a double power-law model through the $\chi^{2}$ minimization. The resulting best fit parameters are summarized in the second line of table 5 and the best-fit function is shown in figure 17 as green dashed line. The best-fit parameters are consistent with each other.

5 Discussion

5.1 Comparison with previous results on the $z=4$ quasar luminosity function

The luminosity function derived in this work is compared with previous work in figure 18. In the bright-end, the binned luminosity function of $z=4$ SDSS DR7 quasar is fully consistent with those at $z=3.75$ with the DR3 sample plotted with blue crosses (Richards et al., 2006) and the DR7 sample plotted with green open triangles (Shen & Kelly, 2012). We convert $M_{i}(z=2)$ used in those papers to $M_{\rm 1450}$ with $M_{\rm 1450}=M_{i}(z=2)+1.486$ (Richards et al., 2006). Thanks to the larger effective area of the DR7 sample, the luminosity function from the DR7 sample extends to higher luminosities than the DR3 sample.

Around the knee of the luminosity function, Palanque-Delabrouille et al. (2013) evaluate the luminosity function based on a quasar sample selected by a variability selection. The luminosity function is described with $M_{g}(z=2)$ and we convert $M_{\rm 1450}=M_{g}(z=2)+1.272$ based on figure 12 of Palanque-Delabrouille et al. (2013). The binned luminosity function covering $3.5<z<4$ is plotted with blue open circles in the figure. The binned luminosity function is fully consistent with the luminosity function of this work within the uncertainty.

In Ikeda et al. (2011), Glikman et al. (2011), Masters et al. (2012), and Niida et al. (2016), the quasar luminosity functions are derived down to $M_{\rm 1450}=-21$ mag in deeper but narrower surveys than the SDSS. These functions are determined in the redshift range $3.7<z<4.7$, $3.7<z<5.1$, $3.5<z<5.0$, and $3.7<z<4.7$ respectively. They cover a higher redshift range with the average redshift of the samples of $z=4.0$, therefore we compare with them after correcting for the evolution effect. In this redshift range, the number density of quasars evolves as $(1+z)^{-6.9}$ (Richards et al., 2006), and we correct for the difference by multiplying 1.2 with the averages of the samples. The binned luminosity functions in the results are shown with pink open diamonds, red crosses, and blue open squares for Niida et al. (2016), Glikman et al. (2011), and Masters et al. (2012) samples, respectively, in figure 18. The binned luminosity functions derived in the COSMOS regions (Niida et al. (2016); Masters et al. (2012)) matches well that of this work, except for the faintest luminosity bin at $M_{\rm 1450}=-21$ mag. They select $z=4$ quasars with stellar morphology and dropout selection method in the similar way to this study. On the other hand, the luminosity function by Glikman et al. (2011) based on NOAO Deep Wide-Field Survey (NDWFS) dataset shows significantly higher number density than the other results. They use similar selection method with dropout color and stellar morphology, but they use loose selection criteria for stellarity based on ground-based imaging data with FWHM of \timeform1", and it would be possible that their sample is significantly contaminated by extended non-AGN objects below $R>22.5$ mag (figure 4 in Glikman et al. (2010)). It should be noted that the number density of LBGs rapidly increases below $i>23$ mag, and they can severely contaminate with the loose stellarity selection as discussed in section 3.5.

Although the binned luminosity functions in the COSMOS regions are consistent with this work, the best-fit luminosity function in Masters et al. (2012) has much steeper slope than that from this work. The best-fit luminosity function is shown with blue dashed line after correcting for the redshift difference by multiplying 1.2. The steeper slope is mostly caused by the excess number counts at $M_{\rm 1450}=-21$ mag, and cannot reproduce the luminosity function brighter than the knee. In section 5.3, we further discuss the difference in the best-fit luminosity functions at $z=4$.

By integrating the best-fit luminosity function, we estimate the UV luminosity density at 1450 Å of the quasars at $z=4$ to be $\epsilon_{\rm 1450}=3.2\times 10^{24}$ erg s*-1* Hz*-1* Mpc*-3*, which is similar to that with the best-fit luminosity function by Masters et al. (2012) ($\epsilon_{\rm 1450}$ = $3.1\times 10^{24}$ erg s*-1* Hz*-1* Mpc*-3* after correcting for the evolution factor of 1.2). Giallongo et al. (2015) derive the UV luminosity function of X-ray selected AGNs and argue that the AGN emissivity of UV ionizing photons could be as high as the value required to keep the intergalactic medium highly ionized. Their UV luminosity function is shown with red pentagons in figure 18. We correct for the redshift evolution by multiplying 1.6 to their number density assuming the middle point of their redshift coverage ($z=4.0-4.5$) and the number density evolution with $(1+z)^{-6.9}$. Their number density is more than one order of magnitude higher than the extrapolation of our best-fit luminosity function at $z=4$. Their estimated UV luminosity density is $18.3\times 10^{24}$ erg s*-1* Hz*-1* Mpc*-3* after correcting for the evolution factor. Our estimate is about 6 times smaller than their value, and suggests that optically-selected stellar blue quasars are not the main contributor to the cosmic reionization.

5.2 Evolution of the quasar luminosity function in the early universe

The evolutionary trend in the luminosity function of quasars above the peak of their number density at $z\sim 2-3$ is a fundamental observable to understand the early growth of SMBHs. In the left panel of figure 19, the binned quasar luminosity function at $z=4$ is compared with those at $z\sim 2.3$ (green open circles; Ross et al. (2013)), $z\sim 2.7$ (pink open pentagons; Ross et al. (2013)), $z\sim 3.2$ (blue open squares; Ross et al. (2013); blue open triangles; Masters et al. (2012)), and $z\sim 5$ (red crosses; McGreer et al. (2013)). In the right panel of figure 19, the ratios of the binned quasar luminosity functions to the best-fit double power-law luminosity function at $z=4$ are shown. In figure 20, we also plot the best-fit double power-law model parameters as a function of redshift. If multiple fitting results with different evolutionary scenarios, for example pure-luminosity and pure-density evolution scenarios, are provided in a paper, we pick up the model parameters that reproduce the binned luminosity function better.

The slopes of the bright-end of the luminosity functions distribute around $\beta=3$ above $z=3$, except for the Luminosity Evolution and Density Evolution model in Ross et al. (2013) and Yang et al. (2016) with steeper bright-end slopes. The value in Ross et al. (2013) is closer to those observed at $z<2$ (e.g. Croom et al. (2009)) and could be affected by quasars in the redshift range around $z=2$. If we compare their binned luminosity functions at $z=2.3$ and $z=2.7$ with that at $z=4$, the slopes above the knee seem to be similar to each other except for the most luminous bins. The systematically flatter bright-end slope above $z=3$ than that at $z<2$ is consistent with the trend reported in Richards et al. (2006), Bongiorno et al. (2007), and Richards et al. (2009). The constant bright-end slope above $z=3$ would suggest that we do not need a luminosity dependent evolution model above the knee to explain the increase of the number density at a fixed luminosity (Steffen et al., 2006).

The evolution of the faint-end shows a different trend. The luminosity functions at $z=2.3$ and $2.7$ have a similar shape to the $z=4$ luminosity function and the best-fit faint-end slope at $z=2.2-3.5$ is $\alpha_{\rm z=2.2-3.5}=-1.29^{+0.15}_{-0.03}$ (Ross et al., 2013), which is consistent with the best-fit value at $z=4$ within the 1$\sigma$ uncertainty. Flat faint-end slope is also supported by Bongiorno et al. (2007), Fontanot et al. (2007), and Niida et al. (2016). On the other hand, the best-fit model at $z=3.2$ has $\alpha_{\rm z=3.2}=-1.73\pm 0.11$ (Masters et al., 2012), which is much steeper than that at $z=4$. Such difference can be seen in the shape of the luminosity functions, the binned $z=3.2$ luminosity function does not show clear knee in it. The faint-end slopes at $z=5$ and $6$ could be steeper than that at $z=4$ (McGreer et al. (2013); Kashikawa et al. (2015)). However, the current quasar luminosity function at $z=5.0$ could not be deep enough to clearly detect the knee of the luminosity function, although the best-fit value for the faint-end would imply a steeper slope, $\alpha_{\rm z=5}=-2.03^{+0.15}_{-0.14}$ (McGreer et al., 2013). Furthermore, the constraint on the faint-end of the quasar luminosity function at $z=6.0$ is based on a few quasars (Kashikawa et al., 2015).

The overall evolutionary trend from $z=2$ to $4$ does not support the higher break luminosities and steeper faint-end slopes with redshifts suggested above $z=5$ (McGreer et al. (2013); Kashikawa et al. (2015)). In the right panel of figure 19, luminosity functions at $z=2.3$, $2.7$, and $5$ show relatively flat distribution of the ratio, which would imply pure density evolution of the quasar luminosity function. Such trend is consistent with the evolutionary trend seen in the X-ray luminosity function at $z>3$ (Vito et al., 2014). The shape of the luminosity function at $z=3.2$ shows a different evolutionary trend. The number ratio increases toward the high-luminosity and low-luminosity ends and the ratio is small in the mid-luminosity range around the knee.

5.3 Comparison with the X-ray AGN luminosity function

Comparing the luminosity function of the optically-selected quasars with that of X-ray-selected AGNs in the same redshift range, we can infer the fraction of blue stellar quasars among the X-ray AGN population, which includes heavily-obscured AGNs. Recent analysis of samples of X-ray selected AGNs indicates that the fraction of heavily-obscured AGNs is higher at $z>3$ than that in the local universe (Hiroi et al. (2012); Vito et al. \yearciteVito2013, \yearciteVito2014), and it is possible that the optical selections are missing such heavily-obscured AGNs. In this discussion, we assume both of the AGN populations have the same overall SED, but the optical selection can be affected by the dust extinction and host galaxy contamination.

We convert the $z=4$ quasar $M_{\rm 1450}$ luminosity function to the hard X-ray 2–10 keV luminosity function based on the relation between the monochromatic luminosity at 2500 Å, $l_{\rm 2500{\AA}}$, and 2 keV, $l_{\rm 2keV}$ for non-obscured broad-line AGNs (e.g. Vignali et al. (2003); Strateva et al. (2005); Steffen et al. (2006); Young et al. (2010)). $l_{\rm 2500{\AA}}$ is calculated from $M_{\rm 1450}$ assuming a power-law spectrum with index of $\alpha=-0.46$. Then, we apply the $l_{\rm 2500{\AA}}$ - $l_{\rm 2keV}$ relation in Steffen et al. (2006), which is derived from 333 optically-selected AGNs, as

[TABLE]

Finally, we convert $l_{\rm 2keV}$ to $L_{\rm 2-10~{}keV}$ assuming the typical X-ray spectrum of broad-line quasars, a power-law spectrum with a photon index of $\Gamma=1.8$.

The resulting hard X-ray luminosity function derived from the optical quasar luminosity function is shown in the left panel of figure 21 as red thin dashed line. The thick solid line indicates the hard X-ray luminosity function of Compton-thin AGNs at $z=3.9$ based on the best-fit luminosity-dependent density evolution (LDDE) model by Ueda et al. (2014). The luminosity function of X-ray-selected AGNs outnumber in the entire luminosity range compared to that of the optically-selected quasars, and the discrepancy is broadly consistent with recent comparisons (Marchesi et al., 2016). However, we need to note that the simple one-to-one conversion from UV to X-ray luminosity (Masters et al. (2012); Marchesi et al. (2016)) may not be appropriate to calculate the converted luminosity function. As discussed in Steffen et al. (2006), at a fixed $l_{\rm 2500{\AA}}$, the $l_{\rm 2keV}$ values show significant scatter around the mean value, and such scatter can broaden the luminosity function towards the bright-end due to the steep slope of the luminosity function. Therefore, we convert the UV luminosity function by assuming that the UV luminosity is the primary parameter of the quasar luminosity and that the X-ray luminosity have a scatter around the above relation. Using the rms scatter of $\log(l_{\rm 2keV})$ measured in each $l_{\rm 2500{\AA}}$ bin (table 5 of Steffen et al. (2006)), we derive the rms scatter as a function of $\log(l_{\rm 2500{\AA}})$.

The resulting hard X-ray luminosity function including the scatter is shown by the red thick dashed line in the left panel of figure 21. We consider the luminosity range above $M_{\rm 1450}<-16.0$ mag, which corresponds to $\log L_{\rm 2-10~{}keV}({\rm erg}\ {\rm s}^{-1})=41.91$ for the above conversions. Due to the minimum luminosity, the faint-end of the converted luminosity function shows decline toward fainter luminosity. The resulting X-ray luminosity function is consistent with the best-fit LDDE model of the X-ray luminosity function at $z=3.9$ above the knee of the luminosity function. The consistency suggests that the current selection of the $z\sim 4$ quasars does not miss a significant fraction of quasars that are affected by heavy-obscuration.

Below the knee of the luminosity function, the number density of the converted luminosity function shows a deficiency compared to the hard X-ray luminosity function. The deficiency can caused by the contribution of obscured and/or less-luminous AGNs which are missed in our selection of $z\sim 4$ quasars based on stellar morphology and blue UV color. Obscured AGNs and less-luminous AGNs, i.e. Seyfert galaxies and low-luminosity AGNs, can have extended morphology dominated by their host galaxies and/or red UV color due to the dust extinction. In the left panel of figure 21, the X-ray AGN luminosity function for the less-absorbed AGNs with $\log N_{\rm H}({\rm cm}^{-2})<22$ is shown with thick dotted line. We assume the column density distribution derived in Ueda et al. (2014). The luminosity function of less-absorbed AGNs matches well with the converted luminosity function of the optically-selected quasars down to $\log L_{\rm 2-10keV}$ (erg  s*-1*)$=43$, supporting the deficiency below the knee is caused by the contribution from the obscured AGNs.

The consistency between the luminosity functions of optically-selected quasars and X-ray-selected AGNs above the knee could show different luminosity dependence of obscured fraction at $z>3$ from that determined with the X-ray-selected AGNs (Vito et al. (2014)). Vito et al. (2014) show that the fraction of obscured AGN with the hydrogen column density, $\log N_{\rm H}$ (cm*-2*), larger than 23, does not strongly depend on luminosity in the luminosity range between $\log L_{\rm 2-10keV}$ (erg s*-1*)$=43$ and 45 with $0.54\pm 0.05$. The upper luminosity corresponds to $M_{\rm 1450}=-26.7$ mag and the luminosity function does not show significant deficiency of optically-selected quasars in the luminosity range. This discrepancy would be explained with the luminosity dependence of the relation between broad-line AGN fraction and hydrogen column density measured in X-ray (Ueda et al., 2003); large hydrogen column density can be observed among luminous blue broad-line quasars (e.g. Akiyama et al. (2000); Merloni et al. (2014)).

We need to note that the above conversion depends on the assumed $l_{\rm 2500{\AA}}$-$l_{\rm 2keV}$ relation. If we apply the relation derived in Young et al. (2010), the converted luminosity functions with their results based on CENSORED and SPECTRA samples show smaller and larger number density than the X-ray luminosity function above the knee, respectively.

We also apply the same conversion to the best-fit model of the $z=4$ luminosity function from Masters et al. (2012) and from Giallongo et al. (2015), and the resulting luminosity functions are shown with blue dot-dashed and green long-dashed lines, respectively, in the right panel of figure 21. Converted luminosity functions with and without the scatter are shown with thick and thin lines, respectively. The resulting luminosity function based on Masters et al. (2012) broadly overlaps with that of the X-ray luminosity function, and they are apparently consistent with each other. However, it should be noted that their quasar sample is selected based on similar criteria to ours, employing blue UV color and stellar morphology. Since we expect a significant contribution by obscured and/or fainter AGN at least below $\log L_{\rm 2-10~{}keV}({\rm erg}\ {\rm s}^{-1})\sim 43$, this would imply an over estimation of their luminosity function at the faint-end. Similarly, if we use the best-fit luminosity function in Giallongo et al. (2015), the converted luminosity function significantly over-predict the number density, especially below the knee of the luminosity function.

6 Summary

We construct a statistical sample of 1,666 $z=4$ quasars based on 339.8 deg2 of $g$, $r$, $i$, $z$, $y$ band imaging data from the S16A-Wide2 release of the HSC-SSP program. The $z=4$ quasar candidates are selected by their stellar morphology and $g$-band dropout selection criteria. Our selection is optimized for selecting stellar objects with as low contamination from extended galaxies as possible. The selection is effective down to $i=24$ mag. Our $g$-dropout color selection criteria is based on the color distributions of SDSS quasars at $z=3.5-4.0$ and Galactic stars, which are one of the most numerous contaminants to the selection.

The number counts of the selected $z=4$ quasars show a monotonic increase in the covered magnitude range, $i=20.0-24.0$, once we correct for the effect of contamination to the selection. The survey selection efficiency and effective survey area are evaluated with libraries of quasar photometric models and a noise model of the HSC stacked images. We evaluate the contamination by Galactic stars and compact galaxies that meet the g-dropout selection criteria. The results indicate that the contamination rate can be more than 50% in the magnitude range fainter than $i>23.5$ mag, and the apparent excess seen in the magnitude range of the raw number counts can be explained with this contamination.

Most of the $z=4$ quasar candidates do not have spectroscopic redshift information, and we estimate photometric redshifts using a Bayesian method based on our library of quasar models with templates . The distribution of the resulting photometric redshifts implies the sample covers the redshift range between 3.6 to 4.3 with mean redshift of 3.9. The distribution of the UV absolute magnitudes, $M_{\rm 1450}$, estimated from the broad-band photometry, covers an absolute magnitude range down to $M_{\rm 1450}\sim-22$ mag, which is more than 2 magnitudes deeper than the SDSS sample. Thanks to the large sample, we determine the faint-end of the $z=4$ quasar luminosity function with unprecedented statistics. In order to extend our luminosity coverage to the bright-end, we include the uniform $z=4$ quasar sample from the SDSS DR7, providing a combined coverage of $M_{\rm 1450}=-22$ to $-29$ mag.

The luminosity function is well described by a double power-law model with a clear break around $M_{\rm 1450}\sim-25$. The bright-end slope, $\beta=-3.11\pm 0.07$ is consistent with those derived at $z=2.2-3.5$ and at $z=5$. The faint-end slope, $\alpha=-1.30\pm 0.05$, is consistent with that at $z=2.2-3.5$ obtained from BOSS, but flatter than that derived at $z=3.2$ with deeper coverage from the COSMOS survey. The overall shape of the $z=4$ luminosity function does not support the higher break luminosities and steeper faint-end slopes at higher redshifts suggested at $z=5$ (McGreer et al., 2013).

We convert the $M_{\rm 1450}$ luminosity function to the hard X-ray 2–10 keV luminosity function using the relation between $l_{\rm 2500{\AA}}$ and $l_{\rm 2keV}$ (Steffen et al., 2006). Once we consider the scatter of the relation, the number density of UV selected quasars matches well with that of the X-ray selected AGNs above the knee of the luminosity function. Below the knee of the luminosity function, the UV-selected quasars show a deficiency compared to the hard X-ray luminosity function. The deficiency can be explained by obscured AGNs.

{ack}

We thank Dr. Akio K. Inoue for providing us with the updated Monte Carlo model of the IGM absorption.

The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University.

This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at http://dm.lsst.org.

The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE) and the Los Alamos National Laboratory.

Based in part on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center, National Astronomical Observatory of Japan.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.