The Obscured Fraction of Quasars at Cosmic Noon

The XMM-SERVS X-ray point-source catalog in the XMM-LSS region (Chen et al. 2018) was chosen for the primary selection of AGNs. The X-ray survey observations were performed by the XMM-Newton satellite over 5.3 deg² of the XMM-LSS survey field. XMM-Newton has three detectors, MOS1, MOS2, and PN. The catalog contains the combined detection of 5242 X-ray point sources using the three detectors in the 0.5–2, 2–10, and 0.5–10 keV bands. Figure 1 shows the 0.5–10 keV flare-filtered exposure time within the survey area, which reaches ∼50 ks per pointing continuously over the survey area. In some regions, deeper X-ray observations are available from other projects as listed in Table 2 in Chen et al. (2018). The survey flux limits over 90% of the total area are 1.7 × 10⁻¹⁵, 1.3 × 10⁻¹⁴, and 6.5 × 10⁻¹⁵ erg cm⁻² s⁻¹ in the 0.5–2, 2–10, and 0.5–10 keV bands, respectively. The flux from each detector in the catalog was derived using an energy conversion factor assuming a power-law continuum with photon index, Γ = 1.7, and galactic absorption column density. The catalog provides the combined source flux calculated using the error-weighted average of the flux estimated by all three detectors. In this paper, the count rates were converted to the expected count rates detected with the PN detector (PN-equivalent count rates) in later discussions.

Zoom In Zoom Out Reset image size

The HSC-SSP is an optical imaging survey performed by the 8.2 m Subaru telescope using the HSC imager (Miyazaki et al. 2018). HSC is a wide-field camera with a field of view of 1.5 deg². The camera is made up of 116 2K × 4K fully depleted back-illuminated CCDs (Kamata et al. 2012) mounted at the prime focus of the telescope. Among the HSC-SSP data sets with three different depths, imaging data of deep and ultradeep layers in the XMM-LSS field are utilized from the S19A and S20A internal releases. Each HSC pointing is shown in Figure 1 in red and blue dashed circles. S19A has better seeing in the HSC r band than that in the S20A data release. The depth of the HSC deep layer is 27.4, 27.1, 26.9, 26.3, and 25.3 mag in the grizy bands (Aihara et al. 2022). The depth and typical seeing size are summarized in Table 1. The depth was calculated from the median value of the 5σ point-source limiting magnitudes of all patches. Image data from HSC were reduced using the HSC pipeline (Jurić et al. 2017; Bosch et al. 2018, 2019; Ivezić et al. 2019). Photometric and astrometric calibration was performed against the first data release of the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS1; Schlafly et al. 2012; Tonry et al. 2012; Magnier et al. 2013; Chambers et al. 2016; Magnier et al. 2020). The HSC images are warped onto predefined grids called "Tracts" and "Patchs". Each tract is approximately 15 ×15 and divided into 9 × 9 patches.

CLAUDS is a deep U-band imaging survey of the HSC deep layer performed by the 3.6 m CFHT with MegaCam (Sawicki et al. 2019). MegaCam (Boulade et al. 2003) is a wide-field camera with 40 2048 × 4612 pixel back-illuminated CCDs, and covers an area of 1.02 deg².

The CLAUDS survey was performed using two u-band filters; u* and U. The two u-band filters are significantly different from each other. The new U-band filter has better transmission than the u*-band filter and has no red leak at 5000 Å, which is observed in the u*-band filter. For the XMM-LSS deep field, the CLAUDS survey was performed entirely using the u*-band filter. The depth of the CLAUDS in the XMM-LSS reaches 26.6 mag (5σ detection in 2'' diameter aperture), while in the ultradeep region, which corresponds to the Subaru XMM-Newton Deep Survey (SXDS; Furusawa et al. 2008) region, reaches 1 mag deeper. The survey properties are summarized in Table 1.

Basic calibration and data reduction of CLAUDS MegaCam data were performed with the Elixir software (Magnier & Cuillandre 2004) at CFHT. Elixir performs the basic data reduction before sending the data to the Canadian Astronomy Data Center to be processed with MegaPipe (Gwyn 2008). MegaPipe performs astrometric calibration against Gaia astrometry while photometry was calibrated against a Sloan Digital Sky Survey (SDSS) u-band photometry, cross checked with synthetic u-band photometry produced using a combination of Pan-STARRS (Magnier et al. 2020) g-band and GALEX NUV photometry.

The VISTA Deep Extragalactic Observations survey (VIDEO; Jarvis et al. 2013) is a deep near-infrared (NIR) imaging survey performed by the 4.1 m Visible and Infrared Survey Telescope for Astronomy (VISTA) at Cerro Paranal with the VISTA InfraRed CAMera (VIRCAM; Dalton et al. 2006). VIRCAM consists of 16 2K × 2K Raytheon VIRGO HgCdTe detectors.

VIDEO data release 5 mosaic images of the XMM-SERVS field are provided in the ESO Phase 3 data archive. ⁹ In the XMM-SERVS field, the mosaic images in each band are separated into three smaller areas designated as XMM1, XMM2, and XMM3. Among them, XMM1 covers the SXDS. The data were reduced at the Cambridge astronomical survey unit (CASU) using the VISTA Data Flow System (VDFS; Irwin et al. 2004). The astrometry and photometry of the survey were calibrated against the Two Micron All Sky Survey (2MASS) point-source catalog (Skrutskie et al. 2006). The final 5σ depths in a 2'' diameter aperture are 24.51, 24.44, 24.12, and 23.77 mag in the Y, J, H, and Ks bands, respectively. The survey properties are summarized in Table 1.

SERVS (Mauduit et al. 2012) is a deep mid-infrared imaging survey performed by the Spitzer space telescope using the Infrared Array Camera (IRAC; Fazio et al. 2004) during the post-cryogenic mission. Only the IRAC channel 1 (3.6 μm) and channel 2 (4.5 μm) are usable due to the high background in the other bands due to the outage of the cryogenic coolant.

SERVS covers five deep multiwavelength extragalactic fields (ELAIS-N1, ELAIS-S1, Lockman Hole, Chandra Deep Field South, and XMM-LSS), in total 18 deg². The mean integration time per pixel is approximately 1200 s, which is close to the confusion limit of the Spitzer IRAC data. The 5σ depths in 3.6 and 4.5 μm bands are 1.9 and 2.2 μJy in 38 diameter aperture. They correspond to 5σ magnitudes of 23.20 and 23.04, respectively. The survey properties are summarized in Table 1.

The mosaic and uncertainty images were retrieved from the NASA/IPAC Infrared Science Archive. ¹⁰ The data were processed at the Spitzer Science Center. The data reduction pipeline performs the standard image reduction and additional detector-specific processing. The images were co-added and reprojected using MOPEX to a pixel scale of 06 pixel⁻¹. Original photometric calibration of the IRAC data was performed using dedicated calibration observations. Crosschecks against the SWIRE survey (Lonsdale et al. 2003) suggest that a correction factor of 1.02 is needed for the 3.6 μm band but not for the 4.5 μm band. We apply the correction during catalog construction.

Similar to other multiwavelength survey fields, there are a large number of spectroscopic redshift measurements in the XMM-LSS region. Spectroscopic redshift measurements within the survey area were compiled from various spectroscopic surveys, including those from the SDSS data releases 9 and 16 (Ahumada et al. 2020; Ahn et al. 2012), VIMOS Public Extragalactic Redshift Survey (VIPERS; Scodeggio et al. 2018), Galaxy and Mass Assembly (Liske et al. 2015), VIMOS VLT Deep Survey (VVDS; Le Fèvre et al. 2013), VANDELS (Garilli et al. 2021), Milliquas (Flesch 2019), MOSFIRE Deep Evolution Field Survey (Kriek et al. 2015), Ultradeep Survey ¹¹ (UDS; McLure et al. 2013; Bradshaw et al. 2013), 3D-Hubble Space Telescope (3D-HST; Brammer et al. 2012; Momcheva et al. 2016), and the SXDS multiwavelength catalog (Akiyama et al. 2015), as well as from individual studies in the SXDS and XMM-LSS regions (Yamada et al. 2005; Geach et al. 2007; Ouchi et al. 2008; Saito et al. 2008; Smail et al. 2008; van Breukelen et al. 2009; Ono et al. 2010; Simpson et al. 2012; Díaz Tello et al. 2013; Melnyk et al. 2013; Yabe et al. 2014; Wang et al. 2016; Menzel et al. 2016; Ono et al. 2018). In total, 294,536 secure spectroscopic redshift records associated with 238,403 unique galaxies and AGN were compiled. The majority of the spectroscopic redshift records are from the objects in the SDSS and VIPERS catalogs, which have an i-band magnitude up to 22.5. However, the faintest magnitude of deep spectroscopic surveys such as VVDS-UDEEP, 3D-HST, and spectroscopic follow-up of X-ray sources in the SXDS from Akiyama et al. (2015) reaches an i-band magnitude of 24.75.

In order to obtain the multiwavelength properties of the optical counterparts of the X-ray sources, this work uses deep imaging data in the 12 photometric filters from the data sets described above. Proper treatment of the point-spread function (PSF) shape and size differences between data sets is needed in order to obtain accurate colors for photometric redshift estimation and SED fitting. This is especially important for the SERVS mid-infrared data set, which suffers from severe blending due to the larger PSF size than the other data sets.

Prior-based PSF-convolved photometry was performed using T-PHOT (Merlin et al. 2015, 2016). T-PHOT uses morphological information of an object in a high-resolution image to measure its flux in images with low spatial resolution. The low-resolution image needs to have the same World Coordinate System (WCS) and the same or integers-times pixel scale as the high-resolution prior image. In this analysis, all data were resampled to the WCS defined in the HSC S19A internal data set. The process of image alignment, background subtraction, and preparation of variance images are described in the following subsections.

3.1.1. Image Alignment and Background Subtraction

3.1.2. PSF Modeling and Convolution Kernel Construction

3.1.3. Pixel–pixel Correlation Correction

It is well known that pixel–pixel correlation occurs during image resampling with SWarp. Resampling smooths the image; thus, the variance in the image becomes artificially smaller than the original image. Therefore, it is necessary to correct the underestimation of the variance.

Fixed-aperture analysis similar to Bielby et al. (2012) was performed on the background variance images to estimate the effects of pixel–pixel correlation. SExtractor was used to produce segmentation images of all CLAUDS and VIDEO images and 1000 fixed apertures of 2'' radius were placed in regions with no source detections. The sky background was measured from the science images in the 1000 apertures while the photometric uncertainties were estimated from the weight images in the same aperture using PhotUtils (Bradley et al. 2020). A 4σ clip was used to remove outlier apertures that may have fallen on the image boundary or bad detector regions. The sky background variance (${\sigma }_{\mathrm{sky}}^{2}$) was estimated by fitting a Gaussian distribution to the distribution of the sky background values and compared to the variance with the median of the variance image in the same aperture (${\sigma }_{\mathrm{wei}}^{2}$). If there is no pixel–pixel correlation, the two values are consistent with each other.

The HSC pipeline has already performed the correction for pixel–pixel correlation and no resampling was done for the reduced image. For the SERVS data set, thanks to the resampling process described in Section 3.1.1, the SERVS variance data preserved correct variance and is consistent with the variance in the sky background. On the other hand, The variances measured in CLAUDS and VIDEO were significantly smaller than the variance measured in the sky background. The correction factor $k={\sigma }_{\mathrm{sky}}^{2}/{\sigma }_{\mathrm{wei}}^{2}$ was estimated in each patch of CLAUDS and VIDEO data.

The correction factors for the VIDEO images are constant over the survey area of the VIDEO survey. A single correction factor determined from the median of the correction factors in all of the VIDEO images was adopted for simplicity and is summarized in Table 2. On the other hand, the correction factor of CLAUDS shows significant variation due to the variation in the survey depth. The correction factor of each patch was applied individually and the median of all correction factors in each tract was used when the patch-level correction factor could not be determined. Figure 2 shows the distribution of the adopted correction factor over the CLAUDS survey area.

Zoom In Zoom Out Reset image size

Table 2. Pixel–pixel Correlation Correction Factor for Variance in the VIDEO Data Set

Band	Correction ($k={\sigma }_{\mathrm{sky}}^{2}/{\sigma }_{\mathrm{wei}}^{2}$)
VIDEO Y	16.78
VIDEO J	11.88
VIDEO H	11.56
VIDEO Ks	8.32

Download table as:  ASCIITypeset image

The primary source catalogs and segmentation images were constructed from HSC S19A r-band images using SExtractor. As summarized in Table 1, the HSC S19A r-band image has the highest resolution among the imaging data sets and is deep enough to detect a large fraction of objects in the other bands. Therefore, the r-band image is used as the high-resolution prior.

The T-PHOT fitting process fails in some regions where bright stars are present in the image. This is likely due to the fact that the SExtractor deblending algorithm divides the star and halo into many individual sources. To solve this problem, sources within the HSC S19A r-band bright star masks (Coupon et al. 2018) were removed. T-PHOT was run twice in each patch to take into account small sub-pixel astrometric offsets. The astrometric offsets were determined in the first run and applied automatically during the second run. The results include catalogs containing the fitting results, model images, residual images, and residual statistics.

Figure 3 shows optical-IR images, model images, and residual images of an optical-infrared counterpart of an X-ray source. Several sources can be seen from the u* to Ks bands but are blended together in the 3.6 and 4.5 μm images. The residual images in the 3.6 and 4.5 μm bands show systematic residuals. This is possibly due to the PSF asymmetry and its variation over the field of view of the IRAC data sets.

Zoom In Zoom Out Reset image size

The T-PHOT catalogs were combined together and duplicate objects in overlapping patch regions were removed from the catalog. The photometric magnitudes and uncertainties in each band were calculated using the zero-points and zero-point uncertainties shown in Table 3. If the signal-to-noise ratio of the measurement is below 2σ, then 2σ upper limits were adopted instead.

Objects that fall into the HSC S20A bright star mask, bad detector regions affected by stray light, or detector defects, sources on the edge of the survey area, and sources with failed fitting results are flagged in the catalog. In addition, sources that are likely local galaxies and were broken up by the detection algorithm were also identified using region files created from the HyperLeda catalog (Makarov et al. 2014). Lastly, sources containing saturated pixels in the HSC images were flagged using the HSC mask images. The galactic reddening value for each object was retrieved from the IRAS reddening map ¹³ in Schlegel et al. (1998). The galactic dust attenuation in each band was calculated using the Galactic dust extinction law (Fitzpatrick 1999).

Because the multiwavelength photometry sample does not cover the entire sample of the X-ray point sources in Chen et al. (2018), we redefined the survey area of the X-ray sample based on the availability of the multiwavelength photometry considering the following conditions:

• 1.  
  The optical-IR counterpart is in the HSC, VIDEO H band, and SERVS-IRAC1 coverage.
• 2.  
  The optical-IR counterpart is not in any bright star masks of HSC nor in the bad regions of the VIDEO H band.

The first condition was imposed to maximize the coverage of multiwavelength photometry, while the second condition was imposed to remove regions where the optical-IR images were affected by image artifacts such as stray light, bright star halos, and edges of the images. Figure 4 shows the distribution of the X-ray sources that meet the above criteria shown as orange symbols. Out of the 5237 XMM-SERVS X-ray sources with an optical-IR counterpart, 3542 X-ray sources are selected as the primary sample for the statistical discussion.

Zoom In Zoom Out Reset image size

The total area within the redefined area was estimated using a Monte Carlo simulation by randomly distributing 100,000 mock data points within the original survey area presented in Chen et al. (2018) and calculating the fraction of data points that satisfy the criteria. The estimated survey area for statistical analysis is 3.52 deg². The area curve in the 0.5–2, 2–10, and 0.5–10 keV bands was calculated by normalizing the maximum area of the area curve presented in Chen et al. (2018) to be 3.52 deg².

In order to check the updated area curve, the logN-logS of the primary sample based on the redefined survey area is compared to that of the entire XMM-SERVS sample. Figure 5 shows the $\mathrm{log}N-\mathrm{log}S$ relation in the 0.5–2, 2–10, and 0.5–10 keV bands for the primary sample (orange) and the original XMM-SERVS (blue). The $\mathrm{log}N-\mathrm{log}S$ based on the redefined survey area is consistent with the $\mathrm{log}N-\mathrm{log}S$ of the original survey area within the Poisson uncertainty. We conclude that the normalized survey area reproduces the survey area of the primary sample well.

Zoom In Zoom Out Reset image size

In order to select the optical counterpart of each X-ray source, first, the identification in Chen et al. (2018) was adopted. The original optical counterparts of 3180 X-ray sources have a corresponding object in the PSF-convolved photometric catalog within a 2'' radius. In Chen et al. (2018), the majority of X-ray sources (2762 sources) were matched with their counterparts in the SERVS catalog. The SERVS-matched counterparts sometimes contain multiple counterparts in the photometric catalog due to the large PSF size of IRAC. For each SERVS-matched source, we examined the number of r-band detected neighbors within a 2'' radius of the counterpart, and 342 sources have multiple optical counterparts suggesting they are blended. The brightest source in the SERVS 3.6 μm band was chosen as the counterpart of the X-ray source among the blended sources. This modification changed the optical counterpart of 23 of the blended sources.

The remaining 362 sources with no corresponding counterpart within the PSF-convolved catalog show no or faint objects in the HSC S19A r band but show a significant detection in the near-IR (NIR) bands. In order to recover the optically faint counterparts, we produced PSF-matched cutouts by matching the PSF to that of the VIDEO H band using the PSF described in Section 3.1.2. Aperture photometry using a 2'' diameter aperture was performed using SExtractor in dual-image mode on the optical and NIR images with the VIDEO H-band image as the detection image. The 2'' diameter aperture contains ∼80% of the PSF flux. Aperture correction factors were calculated from the VIDEO H-band growth curves. The factor was calculated for each patch and applied to the patch individually in order to account for the PSF variation. For the mid-infrared data sets, prior-based PSF-convolved photometry was performed using the VIDEO H-band images as the high-resolution images. Photometry for 282 of the 362 HSC S19A r-band non-detected sources was successfully obtained.

The images of the remaining 80 sources show neither the HSC S19A r-band nor the VIDEO H-band source but some show faint HSC S19A i-band objects close to the detection limit. We do not attempt to recover these sources because the optical identification with such faint sources is uncertain. The photometry of the 282 VIDEO H-band detected sources and 3180 r-band detected sources were combined together into a single catalog. In summary, multiwavelength photometry is obtained for 97.7%, 3462 out of the 3542 primary X-ray sources. A description of the catalog of the primary X-ray AGN sample and multiwavelength photometry in the HSC-Deep XMM-LSS region is provided in the Appendix.

Out of the 3462 primary X-ray sources, 1321 sources have prior spectroscopic redshift measurements. For the remaining X-ray sources, we calculated the photometric redshift using the photometric redshift code LePhare (Arnouts et al. 1999; Ilbert et al. 2006). LePhare estimates the photometric redshift by minimizing the χ² between the observed and model photometry that was derived from template SEDs. For X-ray detected sources, empirical templates of galaxies, local AGN, as well as composite AGN templates constructed from AGN and galaxy templates in Salvato et al. (2011) were used. For X-ray non-detected sources, galaxy models from Ilbert et al. (2009) were used. The model magnitudes were calculated between redshift 0 and 6 with steps of 0.05. We use both SMC and Calzetti extinction laws (Prevot et al. 1984; Calzetti et al. 1994) and fit the reddening as a free parameter with E(B – V) of 0, 0.025, 0.005, 0.075, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.125, 0.15, 0.175, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, and 0.6 mag. No additional emission lines were considered for the AGN models but are added to the galaxy models. LePhare has the capability of estimating systematic zero-point shifts and applying them to the photometry to reduce the deviation from spectroscopic redshifts. The systematic shifts are derived using the spectroscopic redshift sample. The zero-point correction determined from X-ray non-detected sources was adopted to the X-ray detected sources and is shown in Table 4. In both cases, an additional uncertainty (ERR_SCALE) of 0.05 mag was also adopted in order to take into account any additional systematic uncertainty associated with the photometry.

Photometric redshift performance was evaluated based on the median absolute deviation (σ_NMAD: Hoaglin et al. 1983) defined as

$$\begin{eqnarray*}&&{\sigma }_{\mathrm{NMAD}}=1.48\times \mathrm{median}\,\left(\frac{| {z}_{\mathrm{phz}}-{z}_{\mathrm{spec}}| }{1+{z}_{\mathrm{spec}}}\right),\end{eqnarray*}$$

where z_phz and z_spec are the photometric redshift and spectroscopic redshift, respectively. The outlier fraction (f_out) is the fraction of objects whose normalized absolution deviation (z_phot − z_spec)/(1 + z_spec) is larger than 0.15 in the spectroscopic redshift sample. A threshold of 0.15 was adopted following studies in the COSMOS field (Ilbert et al. 2009; Laigle et al. 2016). The photometric redshift performance was evaluated using sources that do not lie within the stellar locus in the g − z against z − 3.6 μm plane

$$\begin{eqnarray*}&&(g-z)-0.5937({\text{}}z-[3.6])\lt -1.7,\end{eqnarray*}$$

where g, z, and [3.6] are the g, z, and 3.6 μm band magnitudes. Galaxy templates were used to evaluate the photometric redshift performance of the X-ray non-detected sources. Similarly, AGN or galaxy templates were used to evaluate the photometric redshift performance of the X-ray detected source. Figure 6 shows the comparison between photo-z and spec-z of X-ray non-detected (left) and X-ray detected sources (right).

Zoom In Zoom Out Reset image size

A photometric redshift scatter (σ_NMAD) of 0.07 and an outlier fraction (f_out) of 26% was achieved for X-ray detected sources with four catastrophic failures where the photometric redshift could not be determined. For comparison, a photometric redshift scatter σ_NMAD of 0.031 and outlier fraction f_out of 3% was achieved for the X-ray non-detected galaxies. The σ_NMAD and f_out for X-ray detected sources are worse than X-ray non-detected sources. Photometric redshifts for AGN can be difficult to accurately determine due to the flat featureless UV-optical continuum of unobscured AGN.

Figure 7 shows the distribution of the spectroscopic and photometric redshift of the primary X-ray sources as a stacked histogram. Above redshift 2, most of the redshift estimates are from the photometric redshift. In summary, 3458 AGN of the primary sample have a spectroscopic redshift or 12-band photometric redshift; thus, we achieved a total of 99.8% redshift completeness thanks to the deep multiwavelength data set.

Zoom In Zoom Out Reset image size

The hydrogen column density ($\mathrm{log}{N}_{{\rm{H}}}$) associated with the nuclear X-ray emission of the primary X-ray source was calculated using the X-ray hardness ratio (HR) and best-redshift estimation assuming an intrinsic AGN X-ray spectrum. The HR is defined as (H − S)/(H + S) where S is the 0.5–2 keV band count rate and H is the 2–10 keV band count rate. In the XMM-SERVS catalog, the flux of each source is derived by combining measurements from multiple detectors with weighting. The PN-equivalent count rate was calculated by dividing the flux by the energy conversion factor (ECF) of the PN used in Chen et al. (2018). As a check, the PN-equivalent count rate was compared with the reported count rate of sources detected with the PN detector in Chen et al. (2018). The two count rates are consistent with each other for the 2–10 keV band but show a systematic offset in the 0.5–2 keV band. A correction factor of 1.12 was applied to the 0.5–2 keV band to make the PN-equivalent count rate consistent with that measured with the PN detector.

A phenomenological AGN model constructed from a linear combination of a cutoff power-law, pexrav reflection component (Magdziarz & Zdziarski 1995), and scattered AGN continuum was assumed as the AGN X-ray spectra.

$$\begin{eqnarray*}\begin{array}{rcl}\mathrm{Model} & := & \mathrm{tbabs}\,\cdot \,(\mathrm{zphabs}\,\cdot \,\mathrm{cabs}\,\cdot \,\mathrm{zcutoffpl}\\ & & +\mathrm{pexrav}+\mathrm{constant}\,\cdot \,\mathrm{zcutoffpl}),\end{array}\end{eqnarray*}$$

where tbabs is the galactic X-ray extinction and zphabs is the absorption associated with the AGN, cabs is the additional Compton scattering, and constant is the fraction of the scattered AGN continuum. The parameters of the phenomenological model were set according to the best-fitted values from Ricci et al. (2017b) based on AGN in the local universe assuming a photon index of Γ = 1.8. For the reflection component, the reflection strength was set to 1.0 and the inclination angle was set to 30°. The fraction of the scattered AGN continuum was set to 1% of the transmitted AGN continuum. The exponential cutoff was set to 381 keV for all components. The elemental abundances were set to the solar abundance in Anders & Grevesse (1989).

The HR model grid was calculated in the redshift range between 0.001 to 6 and $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})$ grid between 20 and 26 using XSPEC (Arnaud 1996) assuming an on-axis response matrix file (RMF) and ancillary matrix files (ARF) of the XMM-Newton PN detector. ¹⁴ Figure 8 shows a comparison between the HR of the primary sample and the HR of the model grid. Spectroscopically identified type 1 broad-line AGN (BL-AGN) are marked with open circles, and their HR shows a mild increasing trend with increasing redshift. Since most of the BL-AGN should have low column density ($\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\lt 22$), the increasing X-ray hardness is unlikely driven by increasing obscuration. The increasing X-ray hardness can be explained by the Compton reflection component at rest frame 10–30 keV shifting into the observed hard 2–10 keV band at high redshift. It should be noted that broad-line classification is only available for a limited sample.

Zoom In Zoom Out Reset image size

The column density was estimated only for AGN above redshift 2 and only between $\mathrm{log}{N}_{{\rm{H}}}\,({\mathrm{cm}}^{-2})\,=\,$20 and 24 where the HR can distinguish the obscuration of CTN-AGN. The redshift dependence of the HR indicates that obscured AGNs with column density less than $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\lt 22$ can hardly be distinguished. The column density can only be constrained for sources detected in both the 0.5–2 and 2–10 keV bands. Only upper or lower limits on the column density can be derived for X-ray sources detected only in a single band. It should be reminded that the lower limit of the column density for sources that are detected only in the 2–10 keV band is mostly larger than $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\gt 23$ above z > 2. Therefore, sources detected only in the 2–10 keV band were assumed to have $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\gt 23$ in the following discussion. Furthermore, some AGN show a smaller HR than that of an AGN model with Γ = 1.8. These sources likely have a softer AGN X-ray spectrum with a larger photon index (Γ > 1.8) and no obscuration. We assign a column density of $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})=20$ to them.

In order to examine the obscured fraction of luminous AGNs above redshift 2, 673 AGN between redshifts 2 and 5 were selected based on the best available redshift (either spec-z or photo-z) as the high-redshift AGN sample. Of the 672 AGN, 251 AGN were detected in both the 0.5–2 and 2–10 keV bands, while 275 and 53 were detected only in the 0.5–2 and 2–10 keV bands, respectively. Finally, 93 AGN were detected only in the 0.5–10 keV band. Within the 672 AGN, 203 (30%) have spectroscopic redshift and more than half of the spectroscopic sample is BL-AGN (71%).

The intrinsic X-ray 2–10 keV luminosity was derived using

$$\begin{eqnarray*}&&{L}_{X}=\mathrm{ext}(z,{N}_{{\rm{H}}})k(z)4\pi {D}_{L}^{2}(z){f}_{X}\end{eqnarray*}$$

by assuming the best available redshift and $\mathrm{log}{N}_{{\rm{H}}}$, where D_L is the luminosity distance, f_X is the observed X-ray flux in the 2–10 keV band, k(z) is the k-correction term, and ext(z, N_H) is the extinction correction in the observed frame. The observed fluxes were calculated from the PN-equivalent count rate in the 2–10 keV band assuming an ECF of 1.26 × 10¹¹ counts⁻¹/erg s⁻¹ cm⁻². ¹⁵ The ECF and extinction correction was calculated using XSPEC by assuming the phenomenological AGN model presented in Section 4.2. The ECF was from the model without absorption while the extinction correction was the ratio between the model without absorption to that with $\mathrm{log}{N}_{{\rm{H}}}\,=\,$ 20–24.

Figures 9 and 10 show the 2–10 keV luminosity and $\mathrm{log}{N}_{{\rm{H}}}$ of the high-redshift AGN. BL-AGN are marked with a black open circle. Most of the AGN detected in both the 0.5–2 and 2–10 keV bands are those with $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\lt 23.5$ and possess quasar-level luminosity ($\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})\gt 44.5$). For AGN detected only in the 0.5–2 keV band, only upper limits can be placed on $\mathrm{log}{N}_{{\rm{H}}}$ and $\mathrm{log}{L}_{X}$ and the upper limits show a large scatter. On the other hand, the lower limits of $\mathrm{log}{N}_{{\rm{H}}}$ and $\mathrm{log}{L}_{X}$ are derived for AGN detected only in the 2–10 keV band. Most of them are located above $\mathrm{log}{N}_{{\rm{H}}}\gt 23$, which implies that they are heavily obscured AGN.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to model the intrinsic number of AGN at z > 2, the functional form of the luminosity and absorption functions were assumed following Ueda et al. (2014) because the size and luminosity coverage of the current sample are not large enough to determine the overall shape of the luminosity function. The hard X-ray AGN luminosity function describes the number density of CTN-AGN with $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})=20\mbox{--}24$. It is expressed as an evolving double power law following the luminosity-dependent density evolution (LDDE) model. The hard X-ray luminosity function of CTN-AGN in the local universe is described as follows:

$$\begin{eqnarray*}&&\displaystyle \frac{d{{\rm{\Phi }}}^{\mathrm{CTN}}({L}_{X},z=0)}{d\mathrm{log}{L}_{X}}=A{\left[{\left(\displaystyle \frac{{L}_{X}}{{L}_{* }}\right)}^{{\gamma }_{1}}+{\left(\displaystyle \frac{{L}_{X}}{{L}_{* }}\right)}^{{\gamma }_{2}}\right]}^{-1},\end{eqnarray*}$$

where A is the normalization of the luminosity function, L_* is the break luminosity, and γ₁ and γ₂ are the slopes of the connected power laws. The luminosity function outside the local universe follows a luminosity and redshift-dependent evolution:

$$\begin{eqnarray*}&&\displaystyle \frac{d{{\rm{\Phi }}}^{\mathrm{CTN}}({L}_{X},z)}{d\mathrm{log}{L}_{X}}=\displaystyle \frac{d{{\rm{\Phi }}}^{\mathrm{CTN}}({L}_{X},z=0)}{d\mathrm{log}{L}_{X}}e({L}_{X},z),\end{eqnarray*}$$

where e(L_X , z) is the evolutionary term, which describes the redshift dependence of the luminosity function following Equation (16) in Ueda et al. (2014).

Ueda et al. (2003) introduced the absorption function f_abs, which is the probability distribution function defined between $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\,=\,$20 and 26. The function is normalized to 1 between $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\,=\,$ 20 and 24 as follows:

$$\begin{eqnarray*}&&{\int }_{20}^{24}{f}_{\mathrm{abs}}({L}_{X},z;{N}_{{\rm{H}}})d\mathrm{log}{N}_{{\rm{H}}}=1.\end{eqnarray*}$$

As shown in Figure 8, the amount of absorption in the column density range between $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\,=$ 20 and 22 cannot be distinguished from the HR used in this analysis. The absorption function was modified by combining the bins between $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\,=\,$20 and 22 together. The definition is as follows:

$$\begin{eqnarray}{f}_{\mathrm{abs}}({L}_{X},z;{N}_{{\rm{H}}})=\left\{\begin{array}{cc}\tfrac{1-\psi ({L}_{X},z)}{2} & [20\leqslant \mathrm{log}{N}_{{\rm{H}}}\lt 22]\\ \tfrac{1}{1+\epsilon }\psi ({L}_{X},z) & [22\leqslant \mathrm{log}{N}_{{\rm{H}}}\lt 23]\\ \tfrac{\epsilon }{1+\epsilon }\psi ({L}_{X},z) & [23\leqslant \mathrm{log}{N}_{{\rm{H}}}\lt 24]\\ \tfrac{{f}_{\mathrm{CTK}}}{2}\psi ({L}_{X},z) & [24\leqslant \mathrm{log}{N}_{{\rm{H}}}\lt 26]\end{array}\right.,\end{eqnarray} \tag{ 1 }$$

where ψ(L_X , z) is the fraction of CTN-AGN with $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\geqslant 22$ and is the ratio of AGN with $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\,=$ 23–24 to those with $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\,=$ 22–23, and f_CTK is the relative number of Compton-thick AGN (CTK-AGN) relative to obscured CTN-AGN, which is fixed to 1. The constraints on f_CTK are discussed in Section 6.3.

The redshift and luminosity dependence of the absorption function can be described with the function ψ(L_X , z). The obscured fraction in the local universe shows a linearly decreasing dependence on the X-ray luminosity:

$$\begin{eqnarray}\begin{array}{rcl}\psi ({L}_{X},z) & = & \min [{\psi }_{\max },\max [{\psi }_{43.75}(z)\\ & & -\,\beta (\mathrm{log}{L}_{X}-43.75),{\psi }_{\min }]],\end{array}\end{eqnarray} \tag{ 2 }$$

where ${\psi }_{\min }$ and ${\psi }_{\max }$ are the minimum and maximum obscured fractions of CTN-AGN, which are determined to be 0.2 and 0.84, respectively (Ueda et al. 2014). The maximum obscured fraction in Ueda et al. (2014) was set to 0.84 since larger values of ψ will return a negative probability for the absorption function in the $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\,=\,$20–21 bin. Our modification of the absorption function in $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})=20-22$ removes this effect; hence, larger values of ${\psi }_{\max }$ are allowed as discussed in the later sections. β is the slope of the decrease, which is set to 0.24 (Ueda et al. 2014), and ψ_43.75(z) is the obscured fraction at $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})=43.75$.

The redshift evolution of the obscured fraction ψ(L_X , z) is described by ψ_43.75(z)

$$\begin{eqnarray}{\psi }_{43.75}(z)=\left\{\begin{array}{cc}{\psi }_{43.75}^{0}{\left(1+z\right)}^{a1} & z\lt 2.0\\ {\psi }_{43.75}^{0}{\left(1+2\right)}^{a1}={\psi }_{43.75}^{2} & z\geqslant 2.0\end{array}\right.,\end{eqnarray} \tag{ 3 }$$

where ${\psi }_{43.75}^{0}$ is the obscured fraction at $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})=43.75$ in the local universe. It was determined to be 0.43 ± 0.03 (Ueda et al. 2014). The redshift dependence parameter a₁ is determined to be 0.48 ± 0.05. We assumed a constant ${\psi }_{43.75}(z=2)={\psi }_{43.75}^{2}$ in the redshift range above redshift 2.

It should be pointed out that ψ_43.75(z) is a parameter used to control the evolution of the obscured fraction ψ(L_X , z) and is not limited to between ${\psi }_{\min }$ and ${\psi }_{\max }$. As defined in Equation (2), ψ_43.75(z) is equal to ψ(43.75, z) only when ${\psi }_{43.75}(z)\leqslant {\psi }_{\max }$, and beyond that would suggest that ψ(43.75, z) has saturated at ${\psi }_{\max }$ with no further evolution with redshift but the obscured fraction of higher luminosity AGN can still be lower. Therefore, the obscured fraction was estimated with ψ(L_X , z) (see Section 5.1).

In order to estimate the obscured fraction at high redshifts, it is important to evaluate the dependence of the effective survey volume as a function of luminosity, redshift, and amount of obscuration. The survey area that is sensitive enough to detect obscured AGN can be smaller than that for unobscured AGN at the same intrinsic luminosity. The survey area function ${\rm{\Omega }}({L}_{X},z,\mathrm{log}{N}_{{\rm{H}}})$ was estimated as a function of AGN 2–10 keV luminosity, redshift, and $\mathrm{log}{N}_{{\rm{H}}}$ using XSPEC (Arnaud 1996). The same phenomenological AGN model, RMF, and ARF calibration files of the PN detector used to estimate the N_H, and the survey area curve defined in Section 3.4 are considered.

Figure 11 shows the survey area at redshift 2.5 as a function of $\mathrm{log}{L}_{X}$ and $\mathrm{log}{N}_{{\rm{H}}}$ in the 2–10 keV band. At a fixed X-ray luminosity, the accessible survey area decreases with increasing column density due to increasing X-ray extinction. Above $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})=24$, the area becomes constant because the scattered and reflected components dominate the spectra.

Zoom In Zoom Out Reset image size

Maximum-likelihood (ML) fitting was performed to estimate the obscured fraction of quasars using the high-redshift AGN sample. The sample for the ML fitting was constructed from the 304 z = 2–5 AGN detected in the 2–10 keV band. For each source, the survey area based on $\mathrm{log}{L}_{X}$, z, and $\mathrm{log}{N}_{{\rm{H}}}$ was calculated in order to evaluate the possibility of detecting that AGN. One AGN was removed since the corresponding survey area is zero. For the AGN detected only in the 2–10 keV band, we assign the column density of $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})=23.5$.

The ML method is a parametric fitting method, which uses the observed parameters of each object without binning. The likelihood function (${ \mathcal L }$) is generally defined as the product of all probability densities (P) in the sample. The probability density is defined as the probability of finding the ith object with $\mathrm{log}{N}_{{\rm{H}}}^{i}$ at $\mathrm{log}{L}_{X}^{i}$ and zⁱ as P_i :

$$\begin{eqnarray*}&&{P}_{i}=\displaystyle \frac{{f}_{\mathrm{abs}}({L}_{X}^{i},{z}^{i};\mathrm{log}{N}_{{\rm{H}}}^{i}){\rm{\Omega }}({L}_{X}^{i},{z}^{i},\mathrm{log}{N}_{{\rm{H}}}^{i})}{{\int }_{20}^{24}{f}_{\mathrm{abs}}({L}_{X}^{i},{z}^{i};\mathrm{log}{N}_{{\rm{H}}}){\rm{\Omega }}({L}_{X}^{i},{z}^{i},\mathrm{log}{N}_{{\rm{H}}})d\mathrm{log}{N}_{{\rm{H}}}},\end{eqnarray*}$$

where f_abs and Ω are the absorption function and survey area function, respectively.

The likelihood function in the logarithmic form ${ \mathcal M }$ is then defined as the sum of the logarithmic probability density of the ith AGN within the fitting sample

$$\begin{eqnarray*}&&{ \mathcal M }({\psi }_{43.75}^{2},\epsilon )=-2\sum _{i}\mathrm{ln}{P}_{i}({\psi }_{43.75}^{2},\epsilon ).\end{eqnarray*}$$

The best-fit parameters are obtained by minimizing the likelihood function over the parameter space of interest. In our case, we set and ${\psi }_{43.75}^{2}$ as free parameters. The 1σ uncertainty of the best-fit parameters was estimated based on the range where the log-likelihood value changes from the minimum by 1.

The ML fitting was performed in two ways, (1) by fitting ${\psi }_{43.75}^{2}$ with fixed = 1.7, which is the parameter determined in the local universe (Ueda et al. 2014) (hereafter we refer to as "1D") and (2) by fitting both ${\psi }_{43.75}^{2}$ and simultaneously (hereafter "2D").

First, the fitting was performed with the maximum obscured fraction (${\psi }_{\max }$) set to be 0.84 following Ueda et al. (2014). The results suggested that the fitting is affected by the choice of ${\psi }_{\max }$. Thus, the maximum obscured fraction ${\psi }_{\max }$ was set to 0.99 instead. The best-fit results from the 1D and 2D cases are shown in Table 5.

Table 5. ML Fitting Results

Parameter	1D	2D
	1.7 (fixed)	1.4±0.3
${\psi }_{43.75}^{2}$	${1.01}_{-0.04}^{+0.03}$	${0.99}_{-0.03}^{+0.04}$

Download table as:  ASCIITypeset image

Figure 12 shows the observed $\mathrm{log}{N}_{{\rm{H}}}$ distribution of the 2–10 keV band detected AGN. The comparison suggests that the fitting results from the 2D case reproduce the observed distribution better than 1D. Figure 13 shows the predicted distribution of redshift, 2–10 keV luminosity, and $\mathrm{log}{N}_{{\rm{H}}}$ based on the best-fitted parameters from the 2D-ML-fit. The predicted distributions reproduce the observed distribution well. Therefore, the results of the 2D-ML fit are adopted for further discussion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It should be pointed out that the parameter ψ_43.75 represents the obscured fraction of moderately luminous quasars with $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})=43.75$; however, the current sample covers the luminosity range of $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})\,\sim \,$44–45. The obscured fraction of luminous quasars was estimated by calculating the intrinsic number of luminous obscured CTN-AGN. The intrinsic number of AGN in each $\mathrm{log}{N}_{{\rm{H}}}$ bin is calculated with the following equation:

$$\begin{eqnarray}\begin{array}{rcl}N & = & \displaystyle \int \int \displaystyle \int {f}_{\mathrm{abs}}(\mathrm{log}{L}_{X},z,{N}_{{\rm{H}}})\displaystyle \frac{d{\rm{\Phi }}(\mathrm{log}{L}_{X},z)}{d\mathrm{log}{L}_{X}}\\ & & {\rm{\Omega }}{\left(1+z\right)}^{3}{d}_{A}{\left(z\right)}^{2}\displaystyle \frac{d\tau }{{dz}}\,d\mathrm{log}{L}_{X}\ {dz}\ d\mathrm{log}{N}_{{\rm{H}}},\end{array}\end{eqnarray} \tag{ 4 }$$

where f_abs, dΦ, Ω, $d\mathrm{log}{L}_{X}$, d_A , and τ are the absorption function, 2–10 keV luminosity function, survey area, angular size distance, and look-back time, respectively. The integration limits are between the redshift, luminosity, and column density of interest. For the expected number of AGN, the survey area Ω is replaced with the survey area function ${\rm{\Omega }}({L}_{X},z,\mathrm{log}{N}_{{\rm{H}}})$.

Based on the best-fit parameters from the 2D-ML fit, the obscured fraction of the luminous quasars with $\mathrm{log}{L}_{X}\ (\mathrm{erg}\,{{\rm{s}}}^{-1})\gt 44.5$ is estimated to be ${0.76}_{-0.03}^{+0.04}$. The lower limit of the luminosity integration was set to $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})=44.5$, as it corresponds approximately to the break in the hard X-ray luminosity function at redshift 2.

Our best-fit obscured fraction is the largest compared to the obscured fraction, which used the same model and in the same luminosity range. Hiroi et al. (2012) estimated the obscured fraction of quasars with $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})\,=\,$44–45 at redshift 3–5 to be ${0.54}_{-0.19}^{+0.17}$, while the best-fit parameters in Ueda et al. (2014) suggest an obscured fraction of 0.50 ± 0.09 above redshift 2.

For comparison with other studies, the obscured fraction of AGN with $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})\,=\,$ 44–45 based on the 2D-ML best-fit parameters is adopted. For CTN-AGN, the obscured fraction is ${0.85}_{-0.03}^{+0.04}$. If we calculate the obscured fraction with $\mathrm{log}{N}_{{\rm{H}}}({\mathrm{cm}}^{-2})\gt 23$, the obscured fraction is ${0.50}_{-0.06}^{+0.07}$.

The ML method applied in Section 5.1 assumes that the slope of the obscured fraction dependence on luminosity does not depend on redshift. Therefore, the rate at that the obscured fraction increases is the same at all luminosity as long as it has not saturated. Here, the slope of the luminosity dependence on the luminous end at high redshift was investigated by applying the ML method to determine the slope parameter β assuming that ${\psi }_{43.75}^{2}=0.73$ following Ueda et al. (2014) and = 1.4 following the 2D-ML-fit results.

The best-fit slope determined using the ML fit is $-{0.09}_{-0.03}^{+0.04}$, which suggests an obscured fraction of luminous CTN-AGN with $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})\,=\,$44–45 is 0.78 ± 0.02. The best-fit slope suggests that the obscured fraction shows almost no luminosity dependence. However, it should be noted that the luminosity coverage is limited to between $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})\,\sim \,$44 and 45, and the obscured fraction at higher luminosity cannot be constrained.

First, systematic uncertainties arise from the photon index and the reflection strength assumed in the phenomenological AGN model, which is used to calculate the column density using the hardness ratio.

If we assume a photon index of 1.7 with the same reflection strength, the obscured fraction for luminous quasars with $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})\,=\,$44–45 is estimated to be ${0.77}_{-0.04}^{+0.03}$, which is approximately 9% lower than when assuming a photon index of 1.8 (${0.85}_{-0.03}^{+0.04}$). Assuming a photon index of 1.8 with a stronger reflection strength of 1.3 reduces the obscured fraction to 0.82 ± 0.03 for luminous quasars, corresponding to a systematic change of approximately 3%. The assumption of the photon index and the reflection strength does not strongly affect the estimate of the column density ratio () since the results are consistent with each other within uncertainties.

Second, the obscured fraction based on a hard X-ray selected AGN sample could suffer from a bias in which the 2–10 keV count rates close to the detection limit are larger than the true count rates due to statistical fluctuations of the photon count rates (Eddington bias). As a result, hard X-ray-selected AGN samples may have harder hardness ratios and as a result, larger column densities overall. This statistical fluctuation may also affect intrinsically unobscured AGN, which makes them show a large column density consistent with obscured AGN due to positive fluctuation in the 2–10 keV band.

Last, the estimates of the column density and luminosity may also be affected by the uncertainty in the photometric redshift. Contamination from low-redshift AGNs (z < 2) can also affect the best-fit results. We estimate the contamination rate from low-redshift AGN (z_spec < 2) based on the fraction of outlier spectroscopically confirmed quasars above redshift 2 (AGN above the red-dashed line with z_phz > 2 as shown in Figure 6) over all AGN above redshift 2 (z_phz > 2) to be 31%.

Using the best-fit parameters from the 2D-ML method, we calculate the obscured fraction of CTN-AGN as a function of the 2–10 keV luminosity by assuming the luminosity-dependent obscured fraction presented in Section 4.4. Figure 14 shows the best-fit obscured fraction based on the 2D-ML best-fit parameters compared with previous measurements in the local universe (bottom panel) (Burlon et al. 2011; Ueda et al. 2014; Georgakakis et al. 2017) and at z = 2–5 (top panel) (Hasinger 2008; Iwasawa et al. 2012; Hiroi et al. 2012; Kalfountzou et al. 2014; Ueda et al. 2014; Liu et al. 2017). Our estimate of the obscured fraction is larger than the obscured fraction in the local universe, which supports the increasing trend in the obscured fraction toward high redshift. At high redshift, our estimate of the obscured fraction at $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})=43.75$ is the largest compared with studies in the same redshift range, except for Liu et al. (2017), in which the obscured fraction in the same redshift range and with $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})=43.5\mbox{--}44.2$ is determined to be 0.91 ± 0.03. The luminosity coverage of their sample is below the luminosity range of our sample but consistent with our estimate of the obscured fraction at $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})=43.75$ if we extrapolate the fraction toward lower luminosity without the maximum limit of 0.84 and compare it with the obscured fraction of AGN with $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\gt 22$.

Zoom In Zoom Out Reset image size

The upper and lower panels of Figure 15 show the obscured fraction of AGN with $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\gt 22$ and >23, respectively. The obscured fraction is larger than the obscured fraction in the same luminosity range as in Ueda et al. (2014) (black solid line). The larger obscured fraction suggests that the evolution of the obscured fraction needs to be stronger at z < 2. In order to reconcile with our estimate of the obscured fraction, the evolution parameter (a1) of the obscured fraction must be ${0.75}_{-0.08}^{+0.10}$. However, this will systematically increase the obscured fraction at z < 2. Alternatively, the obscured fraction may still evolve above redshift 2 and the obscured fraction of high-luminosity AGN saturates at a higher redshift than that of low-luminosity AGN.

Zoom In Zoom Out Reset image size

Vito et al. (2014) estimated the obscured fraction of $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\gt 23$ and $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})\geqslant 43$ AGN at redshift 3–5 to be 0.54 ± 0.05. Our sample resides at a lower redshift than that in Vito et al. (2014). However, our model assumes no evolution of the obscured fraction above redshift two. The obscured fraction of AGN based on the same definition and redshift and luminosity range was estimated to be 0.57 ± 0.05 based on the best-fit parameters of the 2D-ML fit. This is consistent with that in Vito et al. (2014) within 1σ uncertainty.

Liu et al. (2017) also examined the obscured fraction using AGN samples from the C-COSMOS legacy (Civano et al. 2011, 2016) combined with the Chandra Deep Field South 7M catalog (CDF-S; Luo et al. 2017). The obscured fraction of AGN with $\mathrm{log}{N}_{{\rm{H}}}({\mathrm{cm}}^{-2})\ =$ 22–23 among those with $\mathrm{log}{N}_{{\rm{H}}}({\mathrm{cm}}^{-2})\lt 23$. For AGN with $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})\ =\,$44.1–44.9 at redshift 2–3, the fraction is 0.67 ± 0.07. The obscured fraction estimated based on the same definition between $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})\ =\,$44 and 45 using the 2D-ML best-fit parameters is ${0.70}_{-0.01}^{+0.05}$, which is consistent with the above value within 1σ uncertainty.

Recently, Gilli et al. (2022) estimated the obscured fraction of AGNs as a function of redshift by constructing an ISM model based on ALMA data. Our results are consistent with those of Gilli et al. (2022) for AGN with $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\gt 23$ with $\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})\sim 44.$ but higher than the 1σ uncertainties for the obscured fraction of AGN with $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\gt 22$ at the same luminosity. The larger obscured fraction in our study in the case of $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\gt 22$ may be due to the difficulty in separating unobscured AGN from mildly obscured AGN through the hardness ratios.

The increasing trend in the obscured fraction has strong implications for the physical structure and evolution in the nuclear environment of the AGN. The larger obscured fraction at high redshift implies that there is a larger amount of obscuring material within the nuclear, circumnuclear, or host galaxy scale compared to AGN in the local universe.

The trend may be a direct result of the evolution in the host galaxy-scale properties of ISM. Observations of massive galaxies at high redshift show that they are more compact (van der Wel et al. 2014) and have higher gas fractions (Tacconi et al. 2010; Carilli & Walter 2013) compared to those in the local universe with the same stellar mass. This may suggest that the gas density is higher on all spatial scales of the host galaxy compared to those in the local universe. As a result, the higher occurrence of obscuration in the host galaxy or circumnuclear region may explain the increasing trend of the obscured fraction of CTN-AGN (Buchner & Bauer 2017; Buchner et al. 2017; Circosta et al. 2019; Fabian et al. 2008; Gilli et al. 2022) and may explain the Compton-thick fraction (D'Amato et al. 2020). Due to the denser ISM and metal abundance of gas in the circumnuclear region and host galaxy at high redshift, AGN feedback can be less efficient in clearing sight lines since gas can easily cool and replenish the obscuring material (Trebitsch et al. 2019). The decrease in the obscured fraction from high to low redshift may then be explained due to gas consumption by star formation and AGN accretion (Hirschmann et al. 2014). Feedback by AGN-driven winds further reduces the obscuration within and possibly beyond the nuclear region, as spectroscopic observations of high-redshift AGN show that AGN can drive winds with velocities from several hundred up to a thousand kilometers per second, such outflows can reach out to several kiloparsecs beyond the nuclear region (Collet et al. 2016; Nesvadba et al. 2017; Davies et al. 2020).

Another possibility is that the trend in the obscured fraction is driven by the triggering of AGN by major mergers. Galaxy merger simulations suggest that luminous quasars may go through an evolutionary sequence. In this scenario, strong gravitational effects by a merging event funnel large amounts of gas and dust toward the nuclear region, which results in heavily obscured nuclear activity (Hopkins et al. 2006, 2008; Hickox et al. 2009). During the obscured quasar phase, obscuration with a large column density ($\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\gt 22$) as well as near-Eddington limited accretion are induced. This phase is expected to last as long as 10 times the following blowout phase, in which strong winds driven by the AGN blows out the obscuring material leaving an unobscured quasar at the end of the sequence. It is possible that the fraction of AGN triggered by major merger is higher at higher redshifts; observations of merging galaxies in the local universe show that they are heavily obscured compared to those triggered in other processes (Ricci et al. 2017c, 2021).

Lastly, the trend in obscured fraction with redshift may be compatible in the context of radiation pressure on nuclear dust where after exceeding an N_H critical effective Eddington ratio, the dust in the nuclear region is blown out (Fabian et al. 2006, 2008, 2009; Ishibashi & Fabian 2015). In contrast to the merger-driven scenario, mergers are not a prerequisite in order to drive strong outflows. In this model, nuclear obscuration larger than N_H > 5 × 10²¹ cm⁻² occurs from long-lived dust clouds near the AGN while milder obscuration is due to dust lanes outside the SMBH gravitational sphere of influence and independent of the AGN Eddington ratio. At high redshift, accretion activity occurs with higher Eddington ratios than AGN in the local universe (Nobuta et al. 2012; Schulze et al. 2015). The higher average Eddington ratio suggests more AGN are in a blowout phase and less affected by nuclear obscuration. One possibility to explain the increasing trend of the obscured fraction is that the dust abundance in the nuclear region for high-redshift AGN is lower than those in the local universe (Fabian et al. 2009).

CTK-AGN are heavily obscured AGN with $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\gt 24$. Due to the heavy obscuration, the AGN X-ray continuum emission is strongly suppressed; thus, their detection requires the selection of hard X-rays above the rest frame of 10 keV.

Among the high-redshift AGN sample, 53 AGN were detected only in the 2–10 keV band and have lower limits for $\mathrm{log}{N}_{{\rm{H}}}$ larger than $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\geqslant 23$. Considering the band shifting effect, non-detection in the 0.5–2 keV band suggests their AGN X-ray continuum below rest frame 6–8 keV is heavily suppressed. Therefore, the 53 AGN may be considered as CTK-AGN candidates at high redshift.

Assuming the phenomenological AGN model, the survey area function presented in Section 4.5, and using the absorption function in Section 4.4 with f_CTK = 1, the expected number of CTK-AGN detected with $\mathrm{log}{N}_{H}\ ({\mathrm{cm}}^{-2})=\,$24–26 is ${25.61}_{-1.61}^{+0.95}$. This may suggest that roughly half of the 2–10 single band detected AGN may be CTK-AGN candidates. However, the column density determination based on the 0.5–2.0 and 2.0–10.0 keV bands at HR cannot discriminate heavily obscured CTN-AGN and CTK-AGN in the sample as shown in Figure 8. In addition, the adoption of physical torus models, X-ray spectral analysis, or secondary tracers of AGN luminosity are generally required to reliably confirm the CTK nature (Ricci et al. 2017b). A full discrimination of the CTK population is beyond the scope of this paper.

We examine the correspondence between the X-ray obscuration and the SED shapes between UV and IR bands by examining the rest-frame SED of the high-redshift quasars ($\mathrm{log}{L}_{X}\ (\mathrm{erg}\ {{\rm{s}}}^{-1})\gt 44.5$) detected in the 2–10 keV band. For each quasar, the SEDs were constructed by interpolating between the 12-band photometry shifted to rest frame and normalized at 5500 Å. We also include additional IRAC3 (5.8 μm), IRAC4 (8 μm), and MIPS 24 μm bands photometry from the SWIRE data set (Lonsdale et al. 2003). The AGN sample was separated into X-ray unobscured ($\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\lt 22$) and X-ray obscured ($\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\geqslant 22$) AGN based on the column density. The median SED of X-ray unobscured and obscured quasars was constructed by median combining the individual SED of X-ray unobscured and obscured quasars, respectively.

Figure 16 shows the median SED of the high-redshift quasars compared with the type 1 and type 2 QSO SEDs in Polletta et al. (2007) as well as the median SED of high-redshift BL-AGN. The median SED of X-ray unobscured AGN is flat similar to the median SED of the BL-AGN. This is consistent with the expected power-law SED of unobscured AGN. On the other hand, the median SED of X-ray obscured AGN shows a redder UV, optical, and NIR continuum compared to the BL-AGN.

Zoom In Zoom Out Reset image size

Both the X-ray unobscured and obscured AGN show a variety of SED shapes as shown in the 16th and 84th percentile SED distribution. More than 16% of the obscured AGN have a blue UV continuum similar to the median SED of the BL-AGN, while 16% of the unobscured AGN are redder than the median SED of the obscured AGN. This suggests that the correspondence between the UV-optical SED and X-ray obscuration is not strong.

In order to further examine the correspondence between optical properties, UV-optical-NIR color and morphology, and X-ray obscuration, the AGN were separated into four groups based on the rest frame u*– H color and UV morphology. The rest-frame colors were calculated from the SED fitting results of LePhare. The morphology was inferred from the HSC i-band flux ratio between the HSC S20A i-band PSF and cmodel flux. The PSF (cmodel) flux is derived by fitting the PSF (PSF or galaxy) model. If the AGN appears as a point source on the image, then the ratio is expected to approach 1 while extended AGN will have a smaller flux ratio. We consider AGN with a flux ratio larger than 0.95 as point sources.

The distribution of the high-redshift AGN on the rest-frame u*– H color and the flux-ratio plane is shown in Figure 17. Most of the X-ray unobscured AGN have blue rest-frame color and morphology similar to a point source. For the X-ray obscured AGN, approximately 49% show extended morphology while the remaining sources possess morphology consistent with a point source. X-ray obscured AGN also show a broad color distribution where some X-ray obscured AGN have rest-frame optical NIR colors consistent with the X-ray unobscured AGN.

Zoom In Zoom Out Reset image size

The high-redshift AGN sample was divided into four samples based on the color against the morphology plane, and the column density distribution of each sample is shown in Figure 18. Most of the AGN with extended morphology have $\mathrm{log}{N}_{{\rm{H}}}$ consistent with the X-ray obscured AGN with $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\gt 22$. On the other hand, blue point-source AGN (case 1) has a mixture of column density $\mathrm{log}{N}_{{\rm{H}}}$ both consistent with X-ray obscured and unobscured AGN.

Zoom In Zoom Out Reset image size

For each type of morphology, we performed a K-S test between the $\mathrm{log}{N}_{{\rm{H}}}$ distribution for those that are extended and those which are point source. The possibility that the blue and red point source (case 1 and case 2) were drawn from the same parent distribution is rejected at a 0.01 level of significance (${D}_{\max }=0.391$, D_crit = 0.277). However, the number of samples compared between case 1 and case 2 is limited. On the other hand, the possibility that the blue and red extended sources (case 3 and case 4) were drawn from the same parent distribution cannot be rejected at a 0.01 level of significance (${D}_{\max }=0.256$, D_crit = 0.405). This suggests that both red and blue extend sources may contain X-ray obscured AGN.

The largest difference between the cumulative distribution of $\mathrm{log}{N}_{{\rm{H}}}$ of blue extended AGN and red extended AGN comes from the $\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\,=$ 22–23 bin. Obscuration of the largest column densities ($\mathrm{log}{N}_{{\rm{H}}}\ ({\mathrm{cm}}^{-2})\gt 23$) generally occurs on the nuclear scale while moderate obscuration can occur on the kiloparsec scale (Hickox & Alexander 2018). It is possible that red extended X-ray obscured AGN have larger amounts of gas and dust in the host galaxy scale than blue extended X-ray obscured AGN. As a result, the red color can be explained by the larger dust extinction in the ISM. Another possibility is that the blue UV-optical NIR colors found in blue extended X-ray obscured AGN is scattered light from the luminous quasars (Alexandroff et al. 2013; Assef et al. 2016; Alexandroff et al. 2018; Assef et al. 2020) or from ongoing unobscured star formation in the host galaxy and that the AGN is obscured as suggested by the extended morphology of the AGN host galaxy.

In order to examine the relation between the UV spectral properties and X-ray obscuration, the optical spectra of 26 AGN with SDSS spectroscopic data were examined. The morphology of 21 of these AGN is consistent with a point-source object. Three sources have flux ratios >0.90 close to the stellarity threshold. Of the 26 AGN, two AGN have flux ratios consistent with an extended source with one AGN showing a narrow C iv emission line (σ < 1000 km s⁻¹). Among these 26 AGN, two AGN show absorption associated with the C iv emission line. The remaining AGN have spectra consistent with optically unobscured AGN with broad UV emission lines. Broad absorption line quasars (BAL-QSO) are known to be X-ray obscured (Page et al. 2011; Maiolino et al. 2010; Page et al. 2017; Streblyanska et al. 2010); therefore, some of the X-ray obscured AGN with blue UV continuum and point-source morphology may be BAL-QSOs. X-ray obscuration may also be due to warm or ionized absorbers with no dust (Piconcelli et al. 2005; Merloni et al. 2014). This is overall consistent with the concept of the unified model (Antonucci 1993; Urry & Padovani 1995). Detection of broad emission lines and blue continuum suggests that the line of sight toward the broad-line region and the accretion disk is unobscured by dust but may contain ionized or dust-free X-ray absorbers along the line of sight due to the strong UV radiation.

We conclude that the trend in which X-ray unobscured AGN have a flat SED, blue UV-optical color, and point-source morphology while X-ray obscured AGN have a reddened SED and extended morphology is present but the correspondence is not tight. The large variety in the SED shapes may be due to different types of X-ray absorbers, scattered AGN emission, and the variety of host galaxy star formation and dust content.

A large number of high-redshift unobscured quasars are selected based on their optical color and morphology in wide and deep optical imaging surveys. This technique enables us to select faint unobscured AGN but can miss obscured AGN. We examine the relationship between the high-redshift X-ray selected AGN and those with the optical color and morphology criteria as used in Akiyama et al. (2018) and Pouliasis et al. (2022).

Figure 19 shows the color distribution of AGN in the entire X-ray sample compared with those at redshift 3.5–4.1. The color and morphology selection in Akiyama et al. (2018) was constructed to select AGN with point-source morphology at redshift 3.5–4.1. The color selection was designed to minimize contamination by AGN at other redshifts, galaxies at z ∼ 1, and low-mass galactic stars, which are the major contaminants. In this section, point-source morphology is defined based on the adaptive moment measurement (Hirata & Seljak 2003) used in Akiyama et al. (2018) which is available in the HSC-SSP database. We adopt the condition

$$\begin{eqnarray*}&&\displaystyle \frac{{\rm{i}}\_\mathrm{hsmsourcemomentsround}\_\mathrm{shape}11}{{\rm{i}}\_\mathrm{hsmpsfmoments}\_\mathrm{shape}11}\lt 1.1\end{eqnarray*}$$

$$\begin{eqnarray*}&&\displaystyle \frac{{\rm{i}}\_\mathrm{hsmsourcemomentsround}\_\mathrm{shape}22}{{\rm{i}}\_\mathrm{hsmpsfmoments}\_\mathrm{shape}22}\lt 1.1\end{eqnarray*}$$

to define point-source objects. For comparison with the point-source morphology defined based on the flux ratio in Section 6.4, the adaptive moment criteria correspond to objects with flux ratios of ∼0.97. Approximately 16% of AGN at redshift 3.5–4.1 were selected based on color-selection criterion and half of them have point-source morphology. It should be noted that the color criteria are determined for stellar objects brighter than i < 24 mag, and most of the X-ray selected objects are fainter than i > 24 mag.

Zoom In Zoom Out Reset image size

In addition to a modified version of the color selection used in Akiyama et al. (2018), Pouliasis et al. (2022) use the Lyman-break criteria (Ono et al. 2018). The Lyman-break criterion selects 65% of the X-ray selected AGN at redshift 3.5–4.1, including both point-source and extended objects. The remaining objects, which were not selected by the color-selection criteria possess redder colors. We conclude that the AGN selection based on the optical color can miss obscured AGN, which have reddened or host-dominated colors. Moreover, the application of morphological selection excludes obscured AGN, which resides in extended host galaxies.