Piercing through Highly Obscured and Compton-thick AGNs in the Chandra Deep Fields. II. Are Highly Obscured AGNs the Missing Link in the Merger-triggered AGN–Galaxy Coevolution Models?

Junyao Li; Yongquan Xue; Mouyuan Sun; William N. Brandt; Guang Yang; Fabio Vito; Paolo Tozzi; Cristian Vignali; Andrea Comastri; Xinwen Shu; Guanwen Fang; Lulu Fan; Bin Luo; Chien-Ting Chen; Xuechen Zheng

doi:10.3847/1538-4357/abb6e7

1. Introduction

Since the observational establishment that there are tight correlations between the masses of supermassive black holes (SMBHs) and their host-galaxy properties (such as stellar velocity dispersion) in the local universe, how such small-scale SMBHs coevolve with their large-scale host galaxies has become one of the most fundamental problems in understanding the evolution of galaxies (see, e.g., Kormendy & Ho (2013) for a review). Merger-triggered coevolution models (e.g., Sanders et al. 1988; Di Matteo et al. 2005; Hopkins et al. 2006), in which the gas-rich major merger induces both intense star formation and obscured active galactic nucleus (AGN) activity while the subsequent AGN feedback eventually sweeps out the obscuring materials and shuts down the growth of both the SMBH and stellar populations, provide an attractive explanation to how the central AGN communicates with and influences its host galaxy.

Many studies have been devoted to searching for the possible connections between AGN luminosity, obscuration, and host-galaxy properties, such as stellar mass (M_*), star formation rate (SFR), and merger signatures, to test the merger-driven evolutionary models (e.g., Lutz et al. 2010; Mainieri et al. 2011; Schawinski et al. 2012; Chen et al. 2013; Stanley et al. 2015; Donley et al. 2018). However, how AGN activities are triggered and the exact role that mergers/AGNs play in regulating SMBH/galaxy growth are still under debate. The merger fractions are found to be generally low in various AGN populations (typically ≲20%; e.g., Silverman et al. 2011; Kocevski et al. 2012; Schawinski et al. 2012; Lackner et al. 2014; Villforth et al. 2014; Hewlett et al. 2017), even for those obscured quasars (e.g., Zhao et al. 2019) or fast-accreting AGNs (e.g., Villforth et al. 2017; Marian et al. 2019) where we may expect to see a higher incidence of merger signatures (but see Treister et al. 2012). A positive correlation between galaxy-wide star formation and AGN activities has been reported in several works, at least for the luminous populations (e.g., Lutz et al. 2010; Hatziminaoglou et al. 2010; Shao et al. 2010; Rovilos et al. 2012; Rosario et al. 2012; Chen et al. 2013; Dai et al. 2018), but others find a flat relationship (e.g., Stanley et al. 2015; Suh et al. 2017; Schulze et al. 2019) or suggest that SMBH accretion is probably linked to a complex combination of galaxy properties including ${M}_{* }$ , SFR, and morphology (e.g., Rodighiero et al. 2015; Yang et al. 2017; Fornasini et al. 2018; Ni et al. 2019; Yang et al. 2019), especially that the time-averaged black hole accretion rate (BHAR) appears to be only correlated with bulge growth (Yang et al. 2019). The suppression of star formation at high AGN luminosities has been reported only in a few works (e.g., Page et al. 2012; Barger et al. 2015), while Harrison et al. (2012) pointed out that such observed negative AGN feedback may simply be caused by low source number statistics.

Moreover, analyses of the link between AGN obscuration and host-galaxy properties have also presented mixed results. While Lanzuisi et al. (2017) claimed that the hydrogen column density ( ${N}_{{\rm{H}}}$ ) is strongly connected with M_* but not SFR (see also Rodighiero et al. 2015), Lutz et al. (2010) and Chen et al. (2015) suggested a possible correlation between obscuration and SFR indicators. Other studies found no correlation between AGN obscuration and host properties (e.g., Shao et al. 2010; Rosario et al. 2012).

Several factors may be responsible for the contradictory results (see, e.g., Section 3.1 of Xue (2017) and also Section 5 of Brandt & Alexander (2015)), including the limited sample size (e.g., Harrison et al. 2012), the different sample-selection methods (e.g., X-ray versus IR), the adoption of different indicators to trace AGN (e.g., hardness ratio versus ${N}_{{\rm{H}}}$ ) and galaxy properties, how the undetected sources are treated via stacking (e.g., Mullaney et al. 2015), whether the AGN contamination is properly removed through decomposition when calculating the star formation luminosity (especially when performing stacking analyses; see, e.g., Lutz et al. (2010), Rosario et al. (2012), and Ramasawmy et al. (2019)), as well as the influence of AGN variability and the usage of different binning strategies while analyzing the correlation between two parameters that vary on different timescales (e.g., Hickox et al. 2014; Neistein & Netzer 2014; Volonteri et al. 2015; Lanzuisi et al. 2017).

Furthermore, the lack of correlation between AGN and host-galaxy properties may arise because we are looking at the "inappropriate" AGN populations (e.g., Kocevski et al. 2015; Donley et al. 2018). Cosmological simulations suggest that most of the SMBH growth is expected to happen during a phase of heavy obscuration (e.g., Hopkins et al. 2006, 2008), traced by high ${N}_{{\rm{H}}}$ values in the X-ray band. Therefore, highly obscured AGNs (i.e., having ${N}_{{\rm{H}}}\gt {10}^{23}\ {\mathrm{cm}}^{-2}$ ), which are predicted to represent a critical phase in coevolution models where the heavily dust-enshrouded environment, the enhanced star formation activity and active SMBH accretion all happen "together" via mergers (e.g., Springel et al. 2005), may be the "right" AGN population to examine such evolutionary models.

Indeed, some studies have found that the X-ray-selected, most heavily obscured Compton-thick (CT; defined as N_H ≥ 10²⁴ cm⁻²) AGNs exhibit enhanced merger signatures relative to less obscured AGNs (e.g., Kocevski et al. 2015; Koss et al. 2016; Lanzuisi et al. 2018). However, Schawinski et al. (2012) found that 90% of their heavily obscured quasar candidates are hosted in disk galaxies without showing any disturbed signatures, conflicting with other studies.

In addition, the total merger fractions for X-ray-selected highly obscured AGN samples (≈20%–30%; e.g., Kocevski et al. 2015) are found to be significantly lower than that for IR-selected luminous quasars (≈60%–80%; e.g., Fan et al. 2016b; Donley et al. 2018), and their star formation activities (e.g., Georgantopoulos et al. 2013; Lanzuisi et al. 2015) also seem to be more silent than IR-selected dust-obscured AGNs (e.g., Fan et al. 2016a), further raising questions about whether highly obscured AGNs selected from various diagnostics are triggered by different mechanisms or situate in different evolutionary phases.

In this study, we focus on the X-ray-selected highly obscured AGNs, which present the cleanest sample compared to other selection methods (e.g., Brandt & Alexander 2015; Xue 2017) and also ensure the most direct measurements of AGN activity (X-ray luminosity; ${L}_{{\rm{X}}}$ ) and obscuration ( ${N}_{{\rm{H}}}$ ). By systematically analyzing the multiwavelength data for a large dedicated X-ray-selected highly obscured AGN sample (Li et al. 2019b, hereafter Paper I) in the deepest Chandra Deep Fields surveys (CDFs; for a review, see Xue (2017)), we aim at comprehensively exploring (1) the AGN obscuration properties; (2) whether the growth of highly obscured AGNs is isolated in a small nuclear region or somehow linked with host galaxies; (3) the role of merger in igniting highly obscured SMBH accretion; and (4) whether such AGNs are experiencing a blow-out phase that may eventually cause them to evolve to unobscured AGNs, in order to examine whether highly obscured AGNs are the missing link in the merger-triggered SMBH–galaxy coevolution models.

This paper is organized as follows. In Section 2, we describe our X-ray-selected highly obscured AGN sample and the compilation of the multiwavelength data to construct their broadband spectral energy distributions (SEDs). In Section 3, we describe our SED-fitting method to derive AGN and galaxy properties. In Section 4, we present the distributions of our X-ray AGNs on various optical/IR/X-ray color–color diagrams and their implications for AGN obscuration. In Section 5, we discuss the analyses of star formation activity of our AGN hosts, the connections between AGN properties and their host-galaxy growth, the role that mergers play in triggering highly obscured SMBH accretion, and whether highly obscured AGNs are sweeping out the surrounding materials. In Section 6, we summarize the primary conclusions emerging from this work. Throughout this paper, we adopt flat cosmological parameters with H₀ = 70.0 km s⁻¹ Mpc⁻¹, Ω_M = 0.30, and Ω_Λ = 0.70. We define Compton-thin (CN) AGNs as having ${N}_{{\rm{H}}}\lt {10}^{24}\ {\mathrm{cm}}^{-2}$ , and AGNs with ${10}^{23}\ {\mathrm{cm}}^{-2}\lt {N}_{{\rm{H}}}\lt {10}^{24}\ {\mathrm{cm}}^{-2}$ are called highly obscured CN AGNs. The remaining AGNs with ${N}_{{\rm{H}}}\lt {10}^{23}\ {\mathrm{cm}}^{-2}$ are referred to as less obscured AGNs.

2. Multiwavelength Data

One of the primary goals of this work is to characterize the host-galaxy properties of highly obscured AGNs. The most commonly used method to derive galaxy parameters, such as M_* and SFR, is through fitting their SEDs. Among extragalactic surveys, CDFs are among the most extensively investigated fields that enable us to gather a wealth of multiwavelength data from the ultraviolet (UV) to far-infrared (FIR) regimes and compile broadband SEDs for sources of interest (e.g., Gao et al. 2019; Guo et al. 2020). Here, we describe the multiwavelength data sets for our sample.

2.1. X-Ray Data

In Paper I, we systematically analyzed the X-ray spectral and variability properties of a sample of 436 highly obscured AGNs (including 102 CT AGN candidates) selected in the 7 Ms CDF-S (Luo et al. 2017) and 2 Ms CDF-N (Xue et al. 2016) surveys, which are the two deepest Chandra surveys to date. The mean redshift for this sample is 1.88, with 191 sources having spectroscopic redshifts and 245 sources having high-quality photometric redshifts (see Sections 2 and 4.4 of paper I). We performed detailed X-ray spectral modeling and obtained crucial AGN properties such as ${N}_{{\rm{H}}}$ , the observed ( ${L}_{{\rm{X}},\mathrm{obs}}$ ) and intrinsic ( ${L}_{{\rm{X}}}$ ) 2–10 keV luminosities, and fluxes in the rest frame. All the relevant AGN X-ray information is taken from Paper I, and we refer the readers to that work for details of X-ray spectral fitting.

For the purpose of comparison, our analyses also include 492 less obscured AGNs with $\mathrm{log}\,{L}_{{\rm{X}}}\gt 42\ \mathrm{erg}\ {{\rm{s}}}^{-1}$ identified in Paper I. Note that the X-ray spectral fitting model we used in Paper I (i.e., MYTorus; Murphy & Yaqoob 2009) does not allow ${N}_{{\rm{H}}}$ to vary below ${10}^{22}{\mathrm{cm}}^{-2}$ . To derive a column density value for those X-ray unobscured sources, we refit their X-ray spectra by replacing the absorption (MYTZ), reflection (MYTS), and emission line (MYTL) models of MYTorus with the commonly adopted wabs, pexrav, and Gauss models. The new spectral fitting results are consistent with the previous MYTorus-based results on the classification of less obscured and highly obscured AGNs. For sources with ${N}_{{\rm{H}}}\lt {10}^{20}\ {\mathrm{cm}}^{-2}$ , we set their ${N}_{{\rm{H}}}$ values to ${10}^{20}\ {\mathrm{cm}}^{-2}$ .

Note that the depths and sky coverages of multiwavelength surveys significantly drop at the outskirts of the CDF-S and CDF-N. Therefore, we restrict our analyses to the 294 highly obscured AGNs and 250 less obscured AGNs that lie within the central GOODS-S and GOODS-N fields, to ensure reliable SED fitting results (see Figure 1, where the distributions of z, ${L}_{{\rm{X}}}$ and ${N}_{{\rm{H}}}$ of our sample are also shown).

**Figure 1.** Top: highly obscured AGN sample (294 sources) within the GOODS fields (blue region) used in this work (blue points). Sources in Paper I that lie outside the GOODS fields are excluded (orange points). Bottom: distributions of z, ${L}_{{\rm{X}}}$ , and ${N}_{{\rm{H}}}$ for the highly obscured (red) and less obscured (blue) AGN samples used in this work.
Download figure:
Standard image High-resolution image

**Figure 1.** Top: highly obscured AGN sample (294 sources) within the GOODS fields (blue region) used in this work (blue points). Sources in Paper I that lie outside the GOODS fields are excluded (orange points). Bottom: distributions of z, ${L}_{{\rm{X}}}$ , and ${N}_{{\rm{H}}}$ for the highly obscured (red) and less obscured (blue) AGN samples used in this work.
Download figure:
Standard image High-resolution image

2.2. UV and Optical Data

The UV data are taken from the GALEX DR6 catalog.²¹ For the CDF-S, our optical data include U-, B-, V-, R-, I-, and Z-band photometry from the MUSYC survey (Gawiser et al. 2006); the F606W and F814W photometry of the Hubble Space Telescope (HST) from the CANDELS/3D-HST catalog (Skelton et al. 2014); and the F435W, F775W, and F850LP data from the CANDELS multiwavelength catalogs (Guo et al. 2013). We also supplement these data with 18-band Subaru narrowband photometry compiled in Hsu et al. (2014). The optical data in the CDF-N are mainly from Yang et al. (2014), who collected images from Capak et al. (2004) and Ouchi et al. (2009) and presented point-spread function–matched photometry in the U, B, V, R, I, and z' bands in the H-HDF-N. The HST F435W, F775W, and F850LP data are adopted from the GOODS v2.0 catalog (Giavalisco et al. 2004).

2.3. NIR and MIR Data

We combine HST and Spitzer data with ground-based near-infrared (NIR) photometry to construct the NIR to mid-infrared (MIR) SEDs. For the CDF-S, the F098M, F105W, F125W, F140W, and F160W data are collected from Skelton et al. (2014) and Guo et al. (2013). The Spitzer Infrared Array Camera (IRAC) 3.6 μm, 4.5 μm, 5.8 μm, 8.0 μm, MIPS 24 μm, and 70 μm data are adopted from the SIMPLE survey (Damen et al. 2011) and the GOODS-Herschel catalog (Elbaz et al. 2011). We also utilize J₁-, J₂-, J₃-, H_s-, H_l-, K-, and deep K_s-band photometry from the ZFOURGE catalog (Straatman et al. 2016). For the CDF-N, the Spitzer IRAC photometry as well as the J-, H-, K_s-, and H_k-band data are taken from Yang et al. (2014). The detailed description of these data can be found in Table 1 of Yang et al. (2014). The F125W, F140W, F160W data are gathered from Skelton et al. (2014). The Spitzer IRS 16 μm and 24 μm data are taken from Liu et al. (2018a).

2.4. FIR Data

The CDFs had been observed by the PACS and SPIRE instruments aboard the Herschel Space Observatory at FIR wavelengths of 100 μm, 160 μm, 250 μm, 350 μm, and 500 μm. For the CDF-S, we combine the GOODS-Herschel survey (Elbaz et al. 2011) and the HerMES survey (Oliver et al. 2012) to obtain FIR data, which are calculated by adopting the Spitzer MIPS 24 μm positions as priors. For the CDF-N, we use the state-of-the-art "super-deblended" FIR and submillimeter (SCUBA 850 μm data from the James Clerk Maxwell Telescope) photometry presented in Liu et al. (2018a). This advanced super-deblend technique can significantly improve the accuracy of the measured photometry for confused sources.

Following Stanley et al. (2015), for part of the FIR nondetected sources that do not have flux upper limits provided by the catalogs, we use the 100 μm and 160 μm residual maps to derive them.²² For each nondetected source, we randomly extract 1000 aperture photometry measurements in the source-free vicinity (100'') of the source optical position, and calculate the 99.7th percentile of the measured flux distribution as the 3σ upper limit value; see Section 4.3 of Boquien et al. (2019) for details of how CIGALE handles upper limits.

2.5. Construction of Broadband SEDs

For the optical, NIR, and MIR-FIR catalogs, we adopt 0 farcs 5, 1'', and 2'' as the matching radii to cross-match with our X-ray sources using the coordinates of their multiwavelength counterparts (mostly optical ones) provided by the X-ray catalogs, respectively, and combine the matched multiwavelength data to construct the broadband SEDs. We adopt larger matching radii for IR catalogs due to the lower spatial resolution of IR images. When multiple associations are found within the matching radius, we adopt the closest one as the counterpart. The intrinsic (i.e., absorption-corrected) rest-frame 2–10 keV flux is used to represent the X-ray SED. Among the total sample, 164 sources have at least one solid >3σ detection in the aforementioned five Herschel bands. For the remaining sources that lack Herschel detections, their SFRs may not be well-constrained (e.g., Gao et al. 2019). We will discuss the influence of this issue in Section 3.

3. SED Fitting Method and Results

To derive the host-galaxy properties for our sample, we perform multiwavelength SED fitting using X-CIGALE (Yang et al. 2020)—a new release of the SED fitting code CIGALE (Boquien et al. 2019). X-CIGALE has a few important improvements compared with CIGALE. First, it incorporates a new X-ray module that allows us to take advantage of the unique information of AGN intrinsic power provided by X-ray data and to fit SEDs from X-ray to infrared wavelengths. Second, it implements SKIRTOR (Stalevski et al. 2012), a two-phase clumpy torus model where the torus is illuminated by an anisotropic disk to account for AGN emission, which is more realistic than the previous smooth torus model (Fritz et al. 2006) assumed in CIGALE and has been favored by recent simulations and observations (e.g., Stalevski et al. 2012; Ichikawa et al. 2012; Xu et al. 2020).

We adopt the delayed star formation history (SFH), which has shown good performance when recovering the intrinsic galaxy parameters as verified via simulations (e.g., Ciesla et al. 2015). The BC03 stellar population models (Bruzual & Charlot 2003) are adopted to produce galaxy SEDs by assuming the Chabrier initial mass function (IMF; Chabrier 2003), which are then attenuated by the Calzetti et al. (2000) attenuation law and reradiated in IR using the Dale et al. (2014) dust templates. The modified SKIRTOR model (Stalevski et al. 2012; Duras et al. 2017; Yang et al. 2020), which consists of the direct disk radiation in the form of power laws and the reradiation of the clumpy torus surrounding the central source, is adopted to model AGN emission from UV to IR wavelengths. Specifically, the disk SED is modeled as

$\begin{eqnarray}\lambda {L}_{\lambda }=\left\{\begin{array}{ll}{\lambda }^{2} & 8\leqslant \lambda \lt 50\ [\mathrm{nm}]\\ {\lambda }^{0.8} & 50\leqslant \lambda \lt 125\ [\mathrm{nm}]\\ {\lambda }^{-0.5} & 125\leqslant \lambda \lt {10}^{4}\ [\mathrm{nm}]\\ {\lambda }^{-5} & \lambda \geqslant {10}^{4}\ [\mathrm{nm}]\end{array}\right..\end{eqnarray} \tag{ 1 }$

The attenuation of the torus is treated separately from that of the galaxy component. The output torus radiation is calculated based on the 3D radiative transfer code SKIRT (Baes et al. 2011), which is dependent on the assumed geometric structure and density profile of the clumpy materials as well as the inclination angle. The X-ray SED is modeled as a cutoff power law with the photon index being fixed to 1.8 during the fitting, which is the mean value for our sample derived through X-ray spectral fitting in Paper I. The X-ray emission is connected to other wavelengths via the α_ox − L_2500Å relation expressed as α_ox = −0.137 log L_2500Å + 0.2638 (Just et al. 2007). Following Yang et al. (2020), we adopt $| {\rm{\Delta }}{\alpha }_{\mathrm{ox}}|$ , which represents the deviation from the observed α_ox − L_2500Å relation to be 0.2, corresponding to ≈2σ scatter of the α_ox − L_2500Å relation. The summary of the main parameter ranges adopted in the fitting is presented in Table 1.

Table 1. SED Fitting Models and Parameter Spaces Adopted in X-CIGALE

Parameter	Value
Stellar population synthesis model: Bruzual & Charlot (2003)
Initial mass function	Chabrier
star formation history	Delayed τ model
E-folding time of the main stellar population model in Myr	100, 158, 251, 398, 631, 1000, 1584, 2512, 3981, 6309, 10000
Age of the oldest stars in the galaxy in Myr	100, 158, 251, 398, 631, 1000, 1584, 2512, 3981, 6309, 10000
Metallicity	0.02

Galactic dust attenuation: Calzetti et al. (2000)
E(B − V) lines	0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0

Galactic dust emission: Dale et al. (2014)
Power-law slope α	1.5, 2.0, 2.5

Torus model: SKIRTOR (Stalevski et al. 2012; Yang et al. 2020)
Average edge-on optical depth at 9.7 μm	7.0
Angle between equatorial axis and line of sight	30, 70
Half-opening angle of the torus	40
AGN fraction (agn_frac)	0.01, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4
	0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 0.99

X-ray model
Photon index	1.8
Maximum Δα_ox	0.2

Note. See Boquien et al. (2019) and Yang et al. (2020) for model details. We adopt the default values in X-CIGALE for parameters that are not listed in this table. It bears mentioning that, given that our sources are robustly selected as X-ray AGNs with ${L}_{{\rm{X}}}\gt {10}^{42}\ \mathrm{erg}\ {{\rm{s}}}^{-1}$ (see Paper I), we therefore artificially require an AGN component during SED fitting by prohibiting the AGN fraction parameter (defined as the AGN contribution to the total IR luminosity) from being zero, which would allow us to measure an MIR luminosity for each AGN instead of having a zero value. The relevant results in Sections 4.2 and 4.3 are not affected by this choice, as utilizing 271 out of 294 sources that have best-fit agn_frac > 0 while allowing it to take a value of zero would yield the same conclusions. The influence of this forced lower-limit AGN contribution (i.e., 1%) to the total IR luminosity (thus SFR) is also subtle and does not materially affect our SFR-related analyses.

Download table as: ASCII Typeset image

An example of our SED fitting results is displayed in Figure 2. The AGN MIR luminosity is represented by the rest-frame 6 μm luminosity derived from the decomposed AGN component. The galaxy stellar mass (i.e., bayes.stellar_mass) and SFR (i.e., bayes.sfh.sfr) are adopted from CIGALE outputs.

**Figure 2.** Best-fitting SED for CDF-N XID 511 at z = 1.02 (top panel) and fitting residuals (bottom panel), defined as (data-model)/data. Different model components are labeled by different colors. Note that the attenuated stellar emission is largely overlapped with the model spectrum at UV-optical wavelengths and the unattenuated stellar emission at IR wavelengths. Observed fluxes are shown in open squares while the predicted model fluxes at filter wavelengths are shown as filled circles. SFR, stellar mass, E(B-V), and AGN fraction for this source are 7.3 ${M}_{\odot }\,{\mathrm{yr}}^{-1}$ , 10^10.9 ${M}_{\odot }$ , 0.9, and 5%, respectively.
Download figure:
Standard image High-resolution image

Note that the degeneracy between AGN and stellar components encountered during SED decomposition is potentially relevant when the FIR data are absent, and the constraints on SFR become poorer in this situation (e.g., Gao et al. 2019). To validate the usage of sources without solid Herschel detections in our analyses, we carefully test the impact of the lack of FIR data as follows. For Herschel-detected sources, we remove all their Herschel data points and refit the SEDs. We then compare the best-fit SFRs obtained from the Herschel-included fitting to that from the Herschel-excluded fitting, as shown in Figure 3. It can be seen that the SFRs estimated from FIR-data-excluded SEDs are in good agreement with that derived from the whole SEDs. This result shows that, benefiting from utilizing X-ray data in the fitting (which provides a unique insight into the intrinsic AGN power), our good-photometric-coverage optical-to-MIR data are able to provide good constraints on SFR estimates even without resorting to FIR data.

**Figure 3.** Comparison of the SFRs obtained through including (x-axis) and excluding (y-axis) Herschel fluxes in the SED fitting for Herschel-detected sources. Dashed blue lines mark the ±0.5 dex regions that deviate from the one-to-one correlation.
Download figure:
Standard image High-resolution image

4. Revisiting the Optical/IR-selection Methods

We first explore the dependences of SED shapes on AGN physical properties, specifically the X-ray luminosity and obscuring column density. We divide our highly obscured sample into nine ${L}_{{\rm{X}}}$ and ${N}_{{\rm{H}}}$ bins (see Figure 4) and calculate the median composite (AGN + stellar) SED in each bin. The individual SEDs are normalized at rest-frame 1 μm before calculating the median SED, and the results are displayed in Figure 5.

**Figure 4.** Scatter plot of $\mathrm{log}\,{L}_{{\rm{X}}}$ vs. $\mathrm{log}\,{N}_{{\rm{H}}}$ . The ${L}_{{\rm{X}}}$ - ${N}_{{\rm{H}}}$ space is divided into nine bins (low/mid/high- ${L}_{{\rm{X}}}$ vs. low/mid/high- ${N}_{{\rm{H}}}$ ) which are annotated and used to calculate the median SEDs in Figure 5.
Download figure:
Standard image High-resolution image

**Figure 4.** Scatter plot of $\mathrm{log}\,{L}_{{\rm{X}}}$ vs. $\mathrm{log}\,{N}_{{\rm{H}}}$ . The ${L}_{{\rm{X}}}$ - ${N}_{{\rm{H}}}$ space is divided into nine bins (low/mid/high- ${L}_{{\rm{X}}}$ vs. low/mid/high- ${N}_{{\rm{H}}}$ ) which are annotated and used to calculate the median SEDs in Figure 5.
Download figure:
Standard image High-resolution image

**Figure 5.** Composite (stellar + AGN) SEDs for highly obscured AGNs in different ${L}_{{\rm{X}}}$ and ${N}_{{\rm{H}}}$ bins (i.e., S1–S9) defined in Figure 4. Top (Bottom): Comparison of SED shapes at a given ${L}_{{\rm{X}}}$ ( ${N}_{{\rm{H}}}$ ) bin. Each gray curve shows the M82 galaxy template that is attenuated using a specified extinction value (without adding the dust reradiation component in FIR).
Download figure:
Standard image High-resolution image

Comparing the results in different ${N}_{{\rm{H}}}$ and ${L}_{{\rm{X}}}$ bins, we find that the dependence of the composite SED shape on ${N}_{{\rm{H}}}$ is not as sensitive as that on ${L}_{{\rm{X}}}$ at optical-to-MIR wavelengths, which is traced by the large differences between low-luminosity and luminous sources in a given ${N}_{{\rm{H}}}$ bin (e.g., S1 versus S3). We also show the M82 starburst galaxy template being attenuated by three different extinction values in Figure 5 for comparison (the reradiation at FIR is not included). The composite SED shapes for luminous sources are similar to those for typical IR-bright power-law AGNs (e.g., Donley et al. 2012); but for low-luminosity objects, the prominent NIR bump makes their NIR-to-MIR SEDs more similar to those of galaxies whose emission is dominated by dust-obscured star formation (e.g., Riguccini et al. 2015). This overall similarity makes it challenging to identify low-luminosity, highly obscured AGNs using pure SED diagnostics.

Several works have been devoted to using optical and IR colors to select obscured AGN candidates, such as IR-excess methods (e.g., Daddi et al. 2007; Alexander et al. 2008; Luo et al. 2011), WISE-color selection methods (e.g., Tsai et al. 2015; Fan et al. 2016a; Glikman et al. 2018), and IR-to-optical flux-ratio diagnostics (e.g., Fiore et al. 2008); and the studies based on these selection criteria have yielded remarkable insights into our understanding of the obscured AGN population. However, the optical/IR SED-based methods may be biased against low-luminosity AGNs. In the following sections, we will discuss several selection methods in detail. We do not intend to directly quantify the completeness and reliability of each method; instead, we mainly focus on what implications we can deduce by comparing the properties of sources selected using different diagnostics in order to better understand the highly obscured AGN population.

4.1. Can IRAC Colors Effectively Identify Luminous Highly Obscured AGNs?

Among the IR-AGN selection methods, IRAC color is a powerful tool to select large samples of luminous AGN candidates (Stern et al. 2005; Lacy et al. 2007, hereafter L07; Donley et al. 2012, hereafter D12). The most promising aspect of this method compared to X-ray selections (e.g., Xue et al. 2011, 2016; Luo et al. 2017; Xue 2017) is its ability to recover the most heavily obscured sources that are often not detected in X-rays (e.g., 62% of IRAC-selected AGNs do not have X-ray counterparts in D12, which is attributed to heavy X-ray obscuration).

In Figure 6, we plot 247 highly obscured AGNs that are detected/covered in all the four IRAC bands on the IRAC color–color diagram, as well as the color evolutionary tracks at different redshifts calculated from the median composite SEDs in Figure 5. The color evolutionary tracks at z > 1 are located well within the L07 wedge, suggesting that the L07 criterion should be able to efficiently identify highly obscured AGNs. In contrast, almost all the color tracks at z < 3 avoid the refined D12 wedge, suggesting that the X-ray selected highly obscured AGNs will generally be missed by the D12 selection criterion.

When showing (in the right panel of Figure 6) the fractions of our sources being identified as AGNs by the L07 and D12 criteria as a function of ${L}_{{\rm{X}}}$ in two ${N}_{{\rm{H}}}$ bins (corresponding to highly obscured CN and CT sources, respectively), we find that, at $\mathrm{log}\,{L}_{{\rm{X}}}\gt 44.5\ \mathrm{erg}\ {{\rm{s}}}^{-1}$ , the L07 and D12 criteria can recover a substantial fraction of luminous, X-ray-selected highly obscured AGNs with $\overline{\mathrm{log}\,{N}_{{\rm{H}}}}\sim 23.5\ {\mathrm{cm}}^{-2}$ . Such a value is in good agreement with the average column density ( $\mathrm{log}\,{N}_{{\rm{H}}}\sim 23.5\pm 0.4\ {\mathrm{cm}}^{-2}$ ) derived through stacking X-ray-undetected IRAC-selected AGNs in D12 using shallower X-ray data.

However, at $\mathrm{log}\,{L}_{{\rm{X}}}\lt 44.5\ \mathrm{erg}\ {{\rm{s}}}^{-1}$ , the selected source fraction using the D12 criterion dramatically drops even within the range of $44\ \mathrm{erg}\ {{\rm{s}}}^{-1}\lt \mathrm{log}\,{L}_{{\rm{X}}}\lt 44.5\ \mathrm{erg}\ {{\rm{s}}}^{-1}$ , which suggests that it is still incomplete in selecting highly obscured AGNs even for the luminous population (e.g., Kirkpatrick et al. 2017). The $\overline{\mathrm{log}\,{M}_{* }}$ and $\overline{\mathrm{log}\,\mathrm{SFR}}$ values for the D12-missed luminous highly obscured AGNs are 10.9 ${M}_{\odot }$ and 1.3 ${M}_{\odot }$ /yr, which are lower than the D12-selected sources with $\overline{\mathrm{log}\,{M}_{* }}=11.2\ {M}_{\odot }$ and $\overline{\mathrm{log}\,\mathrm{SFR}}=1.7\ {M}_{\odot }$ /yr. Therefore, we suppose that the host-galaxy contamination should not be the main reason responsible for missing a large population of luminous highly obscured AGNs. The average redshift for the missed sources ( $\overline{z}\sim 2.3$ ) being lower than that of the selected sources ( $\overline{z}\sim 3.3$ ) may partly explain the reduced selection efficiency, as can be seen from the color evolutionary tracks. However, we notice that the average redshift for the missed sources is similar to that of IRAC-selected AGNs in D12 (z ∼ 1.8 and z ∼ 2.1 for X-ray-detected and -nondetected sources, respectively). Therefore, the reason that these X-ray-luminous highly obscured AGNs are missed is likely that they are intrinsically fainter in MIR. This may be because the lower dust contents and/or covering factors (CFs) of the tori make their SEDs more similar to that of star-forming galaxies, as verified by their $\overline{\mathrm{log}\,{L}_{6\mu {\rm{m}}}}$ value (∼44.1 $\mathrm{erg}\ {{\rm{s}}}^{-1}$ ) being lower than that of the selected sources (∼45.1 $\mathrm{erg}\ {{\rm{s}}}^{-1}$ ) while the average X-ray luminosities for the two populations are similar ( $\overline{\mathrm{log}\,{L}_{{\rm{X}}}}\sim 44.4\ \mathrm{erg}\ {{\rm{s}}}^{-1}$ ).

For low-luminosity bins, the D12 criterion misses the majority of our sources since the MIR SEDs of low-luminosity AGNs are largely contaminated by the host-galaxy emission (see Figure 5); the more relaxed L07 criterion maintains a relatively high completeness but suffers large contamination from distant star-forming and starburst galaxies (e.g., Donley et al. 2012).

4.2. Is the High Ratio of ${f}_{24\mu {\rm{m}}}/{f}_{R}$ an Efficient Method to Select Highly Obscured AGNs?

Because of large obscuration in highly obscured AGNs, the bulk of UV-optical photons are absorbed and re-emitted in the IR with a peak at MIR. In addition, obscured AGNs tend to have red colors (e.g., Brusa et al. 2005). Consequently, a large ratio of MIR (e.g., Spitzer 24 μm) to optical (e.g., R-band) flux, combined with a red color (e.g., R − K), is expected to be a good tracer of high-level obscuration.

Fiore et al. (2008) (hereafter F08) applied the ${f}_{24\mu {\rm{m}}}/{f}_{R}\gt 1000$ and R − K > 4.5 (in Vega magnitudes, corresponding to 2.86 in AB magnitudes) criteria to the GOODS-MUSIC catalog (Grazian et al. 2006) in order to select the "missing" highly obscured AGN candidates at z ∼ 1.5–2.5 that complement X-ray selections. For the 22 X-ray-detected sources in the 1 Ms CDF-S (Giacconi et al. 2002), the hardness-ratio analysis indicates that they are obscured AGNs with ${N}_{{\rm{H}}}\gt {10}^{22}\ {\mathrm{cm}}^{-2}$ . For the 111 X-ray-undetected sources, the combined stacking analysis and Monte Carlo simulation show that ∼80% of them are possibly highly obscured AGNs. In the era of the 7 Ms CDF-S, with the additional 6 Ms exposure that significantly improves the detectability of heavily obscured sources that are hidden in the previous 1 Ms CDF-S data, we are able to investigate this method in more detail.

In Figure 7, we plot our highly obscured AGNs on the ${f}_{24\mu {\rm{m}}}/{f}_{R}$ versus R − K (in AB magnitudes) diagram using the fluxes predicted at filter wavelengths (i.e., red filled circles in Figure 2), as well as the color evolutionary tracks similar to those in Figure 6. The choice of using model-predicted fluxes instead of observed fluxes here is made in order to enlarge the sample being investigated, i.e., the whole sample can be plotted and we can include each source, even the 101 sources not covered in all bands, in the red, green, or blue populations defined by the shaded regions in Figure 7. We note that using actual observed fluxes to derive colors does not affect our conclusion qualitatively, although the exact values of colors will be slightly different.

The expected correlation between ${f}_{24\mu {\rm{m}}}/{f}_{R}$ and R − K color can be clearly seen (e.g., Fiore et al. 2008), and our sources indeed have a much redder mean color (i.e., ${\rm{\Delta }}\overline{R-K}=1.9$ ) compared to the remaining sources in the SIMPLE survey (Damen et al. 2011).

There are 46 (16%) sources located in the red region defined by F08 (i.e., ${f}_{24\mu {\rm{m}}}/{f}_{R}\gt 1000$ and R − K > 2.86), with 40 of them at z > 2 and the remainder at 1 < z < 2, indicating that this criterion can indeed select heavily obscured AGNs. However, the average redshift for our "red" sources ( $\overline{z}=2.7$ ) is significantly higher than that of the highly obscured AGN candidates selected in F08, which peaks at z = 1.5–2.0. The very low fraction of our sources residing in the red region and the large redshift discrepancy suggest either large incompleteness of this method (Comastri et al. 2011; Brightman & Ueda 2012) or an essential difference between X-ray- and IR-selected populations (e.g., Hickox et al. 2009).

Note that the $\overline{\mathrm{log}\,\mathrm{SFR}}$ of red sources is slightly lower (i.e., ${\rm{\Delta }}\overline{\mathrm{log}\,\mathrm{SFR}}\sim -0.2$ dex) than that of blue ones (i.e., having ${f}_{24\mu {\rm{m}}}/{f}_{R}\lt 1000$ and R − K < 2.86) with matched redshifts, hence the increased ${f}_{24\mu {\rm{m}}}$ of red sources is not primarily caused by the enhanced star formation, but should instead be related to the central AGN. Even if we only consider the most luminous sources (i.e., with $\mathrm{log}\,{L}_{{\rm{X}}}\gt 44\ \mathrm{erg}\ {{\rm{s}}}^{-1}$ ) in order to avoid host contamination of the observed colors, most of them (53/76) still avoid the red region.

To understand the differences between red, blue, and green (i.e., having ${f}_{24\mu {\rm{m}}}/{f}_{R}\lt 1000$ and R − K > 2.86) sources, we annotate their average source properties in Figure 7. It can be seen that, at similar redshifts, the average ${N}_{{\rm{H}}}$ and ${L}_{{\rm{X}}}$ for the three source populations are roughly the same, but red sources have $\overline{{L}_{6\mu {\rm{m}}}}$ significantly higher than those of blue and green sources. Aside from the diverse galaxy contributions to the observed colors, another explanation for the widely distributed colors of our sources in a given redshift bin could be that the dust contents and CFs of the tori for blue and green sources are smaller than those of red ones, resulting in weaker reprocessed MIR emission and smaller ${f}_{24\mu {\rm{m}}}/{f}_{R}$ . Alternatively, if a significant portion of the heavy X-ray obscuration is contributed by dust-free materials, such as broad-line region (BLR) gas and/or disk wind (e.g., Burtscher et al. 2016; Liu et al. 2018b; Ichikawa et al. 2019), the UV-optical continuum will not be significantly attenuated, leading to smaller values of ${f}_{24\mu {\rm{m}}}/{f}_{R}$ and R − K. It is also possible that the interstellar medium (ISM) may contribute significantly to X-ray absorption even up to ${N}_{{\rm{H}}}\gt {10}^{23}\ {\mathrm{cm}}^{-2}$ for high-redshift gas-rich galaxies (e.g., Gilli et al. 2014; Shu et al. 2018; Circosta et al. 2019; D'Amato et al. 2020); if so, since the dust temperature in the ISM is much lower than that in the torus, the reprocessed emission will peak at longer wavelengths (e.g., FIR-to-submm), thus the ${f}_{24\mu {\rm{m}}}/{f}_{R}$ value may not be that large. Indeed, the z = 4.75 CT AGN (Gilli et al. 2011) reported in the 4 Ms CDF-S (XID 403; Xue et al. 2011), whose ISM in the central starburst region (r_half ∼ 0.9 ± 0.3kpc) is able to produce ${N}_{{\rm{H}}}\sim (0.3\mbox{--}1.1)\times {10}^{24}\ {\mathrm{cm}}^{-2}$ as revealed by ALMA observations (Gilli et al. 2014), does have a very low value of ${f}_{24\mu {\rm{m}}}/{f}_{R}$ = 168, especially considering that its very high redshift is supposed to make it easier to fulfill the F08 criteria (see the color evolutionary tracks in the right panel of Figure 7).

In conclusion, we find that the heaviest X-ray obscuration is not equivalent to having extremely large ${f}_{24\mu {\rm{m}}}/{f}_{R}$ and the reddest color, possibly owing to the diverse properties of obscuring materials (e.g., different CFs, gas/dust contents, and ${N}_{{\rm{H}}}$ distributions), complex origins of the X-ray obscuration along our sightline (e.g., X-rays absorbed by dust-free BLR gas, disk wind, dusty torus, and/or ISM), and galaxy contamination to the observed colors.

4.3. Can the Value of ${L}_{6\mu {\rm{m}}}/{L}_{{\rm{X}},\mathrm{obs}}$ Be Used As a Reliable Indicator of ${N}_{{\rm{H}}}$ ?

Since the AGN MIR emission produced by the absorption and reradiation of UV-to-optical photons is largely unaffected by dust attenuation, whereas X-ray photons will be significantly absorbed when ${N}_{{\rm{H}}}$ reaches the highly obscured regime, a large ratio of the MIR luminosity to observed X-ray luminosity ( ${L}_{6\mu {\rm{m}}}/{L}_{{\rm{X}},\mathrm{obs}}$ ) has been widely adopted as an indicator of heavy obscuration (e.g., Alexander et al. 2008; Del Moro et al. 2013, 2016; Rovilos et al. 2014; Corral et al. 2016).

In Figure 8, we show the dependence of ${L}_{6\mu {\rm{m}}}/{L}_{{\rm{X}},\mathrm{obs}}$ on ${N}_{{\rm{H}}}$ for our sample. We confirm that there is a positive correlation between the two parameters (with Spearman's ρ = 0.40 and p ≪ 0.001), albeit with large dispersion. Considering the theoretical argument proposed by Yaqoob & Murphy (2011) that ${L}_{6\mu {\rm{m}}}/{L}_{{\rm{X}},\mathrm{obs}}$ is more sensitive to the torus CF and the incident X-ray continuum shape, rather than ${N}_{{\rm{H}}}$ , it is possible that highly obscured AGNs with small CFs and hard spectral shapes may have lower ${L}_{6\mu {\rm{m}}}/{L}_{{\rm{X}},\mathrm{obs}}$ values (see the sources in regions B and D in Figure 8 that lie below the best-fit line) and high-CF, less obscured AGNs with soft X-ray spectra may have ${L}_{6\mu {\rm{m}}}/{L}_{{\rm{X}},\mathrm{obs}}$ values as large as CT AGNs (see the sources in region A versus those in region D).

These statements are supported by the result that, when we plot (in Figure 8) the $\overline{\mathrm{log}\,{L}_{6\mu {\rm{m}}}/{L}_{{\rm{X}},\mathrm{obs}}}$ and $\overline{\mathrm{log}\,{N}_{{\rm{H}}}}$ values for the z > 1 red, green, and blue populations defined in Figure 7 and Section 4.2, it can be clearly seen that red sources show the highest $\overline{{L}_{6\mu {\rm{m}}}/{L}_{{\rm{X}},\mathrm{obs}}}$ at a given ${N}_{{\rm{H}}}$ , consistent with the scenario that they are deeply buried by plentiful dusty materials. In contrast, for blue and green sources, the dust contents and CFs might be lower, resulting in smaller values of ${L}_{6\mu {\rm{m}}}/{L}_{{\rm{X}},\mathrm{obs}}$ even though the levels of their line-of-sight (LOS) X-ray obscuration are indistinguishable from red sources. Therefore, we conclude that a simple ${L}_{6\mu {\rm{m}}}/{L}_{{\rm{X}},\mathrm{obs}}$ value alone is not sufficient to identify CT AGNs (also see Georgantopoulos et al. 2011).

Even so, the ${L}_{6\mu {\rm{m}}}/{L}_{{\rm{X}},\mathrm{obs}}$ value may still be useful to distinguish highly obscured and less obscured AGNs. In the bottom panel of Figure 8, we show the ${L}_{{\rm{X}},\mathrm{obs}}$ versus ${L}_{6\mu {\rm{m}}}$ relation at three ${N}_{{\rm{H}}}$ ranges by including less obscured AGNs in the plot. The majority (∼60%) of highly obscured AGNs are separated from less obscured ones, while sources with $\mathrm{log}\,{N}_{{\rm{H}}}\lt 22\ {\mathrm{cm}}^{-2}$ are mixed with those with $22\ {\mathrm{cm}}^{-2}\lt \mathrm{log}\,{N}_{{\rm{H}}}\lt 23\ {\mathrm{cm}}^{-2}$ due to the fact that the rest-frame 2–10 keV flux is not significantly absorbed in the Compton-thin regime. We use the linear support vector machine algorithm built into the Python package scikit-learn to derive the boundary line between highly obscured and less obscured AGNs. The optimal boundary is shown as a black line, parameterized as

$\begin{eqnarray}&&{L}_{{\rm{X}},\mathrm{obs}}=0.93\times {L}_{6\mu {\rm{m}}}+2.2.\end{eqnarray} \tag{ 2 }$

This boundary line serve as a complementary method to the hardness ratio criterion we presented in Section 4.3 of Paper I to select highly obscured AGNs.

5. Are Highly Obscured AGNs the Missing Link in the Merger Models?

5.1. Do Highly Obscured AGNs Show Enhanced Star Formation Activities?

To evaluate the star-forming activity of our highly obscured AGN sample in the context of the general galaxy population, we construct a control nonactive galaxy sample from Santini et al. (2015) (hereafter S15) in which the SED-fitting results for 34,929 galaxies from ten independent teams adopting different model configurations are available. The X-ray AGNs identified in the 7 Ms CDF-S catalog are excluded. Following Yang et al. (2017), we adopt the median values of ${M}_{* }$ and SFR reported from teams 2a_τ, 6a_τ, 11a_τ, 13a_τ, and 14a in the following analyses, all of which assumed the same BC03 stellar templates and Chabrier IMF as in this paper.

Note that, although S15 did not consider the AGN component in the SED fitting, their results may still provide reliable mass estimates for highly obscured AGNs with moderate luminosities, as their rest-frame optical-to-NIR SEDs (which are the most important to constrain ${M}_{* }$ ) are largely dominated by the galaxy component (e.g., Luo et al. 2010; Xue et al. 2010). Therefore, we compare our ${M}_{* }$ estimates with S15 for the 82 common sources in the two works with redshift difference Δz < 0.05. The derived $\overline{{\rm{\Delta }}\mathrm{log}\,{M}_{* }}$ between the two works is 0.06 ± 0.02 dex, suggesting that there is no significant systematic bias induced by the different SED-fitting approaches.

In Figure 9(a), we plot our highly obscured AGNs in the SFR versus M_* plane. Also shown are less obscured AGNs with a roughly matched redshift distribution (see Figure 1) and normal galaxies from S15 at similar redshifts. The distributions of SFR of the two AGN populations suggest that they are mainly hosted by star-forming galaxies, and there is no noticeable enhancement of star-forming activity in highly obscured AGNs versus less obscured AGNs (e.g., Suh et al. 2019; Zou et al. 2019).

To further control the redshift and M_* dependence of SFR, we divide the M_* − z space (Figure 9(b)) into a series of subgrids with ΔM_* = 0.2 dex and Δz = 0.2. In each subgrid, the number of highly obscured AGNs is denoted as N_i and we randomly select N_i highly obscured AGNs and N_i normal galaxies while allowing duplication. By repeating this procedure in each subgrid, new highly obscured AGN and normal galaxy samples with matched z and M_* distributions can be constructed. Note that the fraction of galaxies hosting an AGN dramatically increases with ${M}_{* }$ (e.g., Xue et al. 2010; Yang et al. 2017, 2018); thus, at the highest-mass end, we may not be able to find a sufficient number of control galaxies that do not contain an active SMBH, as is the case for our $\mathrm{log}\,{M}_{* }\gt 11.2$ ${M}_{\odot }$ sources (see Figure 9(b)). Therefore, we restrict our sampling procedure to $\mathrm{log}\,{M}_{* }\lt 11.2\ {M}_{\odot }$ , which accounts for 85% of our sample. We perform the above procedures 1000 times and show one example of the distributions of M_* and z for a randomly selected AGN sample and a control galaxy sample in Figure 9(c) and (d). As can be seen, the M_* and z distributions have been controlled well. The random samples vary every time we repeat the sampling procedure. For each sampling, we calculate the 20th, 50th, and 80th percentiles from the SFR distributions of the randomly selected AGNs and control galaxies. The average SFR ( $\overline{\mathrm{log}\,\mathrm{SFR}}$ ) at each percentile is calculated by averaging the values of the 1000 random samples, and the results are summarized in Table 2. The respective 1σ uncertainty is calculated from the 84th percentile and the 16th percentile of the corresponding resampled parameter distribution.

Table 2. Comparison of Star Formation Properties Between Highly Obscured AGNs and Their M_*- and z-controlled Normal Galaxies

Sample	$\overline{\mathrm{log}\,{\mathrm{SFR}}_{\mathrm{agn},50\mathrm{th}}}$	$\overline{\mathrm{log}\,{\mathrm{SFR}}_{\mathrm{gal},50\mathrm{th}}}$	$\overline{\mathrm{log}\,{\mathrm{SFR}}_{\mathrm{agn},20\mathrm{th}}}$	$\overline{\mathrm{log}\,{\mathrm{SFR}}_{\mathrm{gal},20\mathrm{th}}}$	$\overline{\mathrm{log}\,{\mathrm{SFR}}_{\mathrm{agn},80\mathrm{th}}}$	$\overline{\mathrm{log}\,{\mathrm{SFR}}_{\mathrm{gal},80\mathrm{th}}}$
Total	${0.87}_{-0.02}^{+0.02}$	${0.90}_{-0.13}^{+0.11}$	${0.20}_{-0.15}^{+0.12}$	$-{0.63}_{-0.14}^{+0.17}$	${1.40}_{-0.07}^{+0.04}$	${1.70}_{-0.05}^{+0.06}$
high ${L}_{{\rm{X}}}$	${0.91}_{-0.03}^{+0.02}$	${0.92}_{-0.18}^{+0.17}$	${0.49}_{-0.11}^{+0.06}$	$-{0.73}_{-0.16}^{+0.20}$	${1.31}_{-0.07}^{+0.04}$	${1.75}_{-0.11}^{+0.06}$
low ${L}_{{\rm{X}}}$	${0.75}_{-0.08}^{+0.09}$	${0.87}_{-0.17}^{+0.17}$	${0.00}_{-0.14}^{+0.06}$	$-{0.46}_{-0.30}^{+0.17}$	${1.40}_{-0.07}^{+0.04}$	${1.65}_{-0.06}^{+0.07}$

Download table as: ASCII Typeset image

The $\overline{\mathrm{log}\,{\mathrm{SFR}}_{50\mathrm{th}}}$ for highly obscured AGN hosts is found to be consistent with that of control galaxies ( $\bar{{\rm{\Delta }}{\rm{log}}\,{{\rm{SFR}}}_{50{\rm{th}}}}\,=-{0.03}_{-0.12}^{+0.12}$ ). Such a result also holds for both low- ( $\bar{{\rm{\Delta }}{\rm{log}}\,{{\rm{SFR}}}_{50{\rm{th}}}}=-{0.11}_{-0.19}^{+0.15}$ ) and high-luminosity ( $\bar{{\rm{\Delta }}{\rm{log}}\,{{\rm{SFR}}}_{50{\rm{th}}}}\,=-{0.02}_{-0.16}^{+0.20}$ ) AGNs if we split the sample into two subsamples based on the median ${L}_{{\rm{X}}}$ at each redshift grid. The $\overline{\mathrm{log}\,{\mathrm{SFR}}_{80\mathrm{th}}}$ for highly obscured AGN hosts appears to be lower than that of control galaxies with $\bar{{\rm{\Delta }}{\rm{log}}\,{{\rm{SFR}}}_{80{\rm{th}}}}=-{0.30}_{-0.07}^{+0.08}$ . These results suggest that the star-forming activity of highly obscured AGNs is not enhanced with respect to normal star-forming galaxies.

However, a deficiency of quiescent hosts among highly obscured AGNs can be clearly seen with $\bar{{\rm{\Delta }}{\rm{log}}\,{{\rm{SFR}}}_{20{\rm{th}}}}\,\approx {0.83}_{-0.20}^{+0.17}$ . As a result, the SFR distribution of highly obscured AGNs appears to be less diverse and more main-sequence, like than that of normal galaxies (e.g., Bernhard et al. 2019), suggesting that the sufficient cold gas supply that is responsible for sustaining star formation may also be an important factor in triggering highly obscured AGNs. On the other hand, the indistinguishable $\overline{\mathrm{log}\,{\mathrm{SFR}}_{50\mathrm{th}}}$ for highly obscured AGN hosts to that of control galaxies suggests that, unlike the significant populations of IR-selected ultraluminous IR galaxies and hot dust-obscured galaxies, which are generally believed to hold both highly obscured AGN (e.g., Vito et al. 2018) and enhanced starburst activity as a consequence of mergers (e.g., Farrah et al. 2003; Fan et al. 2016a, 2016b), there is no evidence supporting the notion that the presence of X-ray-selected highly obscured AGNs is more frequently connected to violent, possibly merger-driven starburst activities (e.g., Georgantopoulos et al. 2013; Lanzuisi et al. 2015; Suh et al. 2017).

5.2. Are AGN Activity and Obscuration Linked with Host-galaxy Properties in Highly Obscured AGNs?

Many studies have explored the correlations between AGN and host-galaxy properties (e.g., ${L}_{{\rm{X}}}$ , ${N}_{{\rm{H}}}$ versus M_*, SFR) in a variety of redshift ranges, but there has been considerable debate about whether AGN activity and obscuration are linked with galaxy-wide star formation (e.g., Lutz et al. 2010; Shao et al. 2010; Harrison et al. 2012; Page et al. 2012; Stanley et al. 2015; Lanzuisi et al. 2017; Dai et al. 2018; Schulze et al. 2019), as well as whether M_* or host-galaxy compactness is a more fundamental factor that governs the average black hole accretion rate (BHAR) (e.g., Yang et al. 2017, 2018; Fornasini et al. 2018; Ni et al. 2019). In particular, Yang et al. (2019) presented an attractive scenario in which the SMBH only coevolves with the galaxy bulge as traced by the significant correlation between $\overline{\mathrm{BHAR}}$ and $\overline{\mathrm{SFR}}$ in bulge-dominated galaxies, while for non-bulge-dominated galaxies, the $\overline{\mathrm{BHAR}}$ is not linked with $\overline{\mathrm{SFR}}$ but instead is predominantly determined by ${M}_{* }$ .

However, we note that, in Yang et al. (2017) (Y17) and Yang et al. (2019) (Y19), the $\overline{\mathrm{BHAR}}$ (calculated from $\overline{{L}_{{\rm{X}}}}$ by averaging ${L}_{{\rm{X}}}$ for both X-ray-detected and X-ray-undetected galaxies) is derived by assuming a wabs × zwabs × powerlaw model, which, as shown in Section 4.1 of Paper I, is inappropriate for highly obscured AGNs because it will significantly underestimate ${L}_{{\rm{X}}}$ , owing to the negligence of the Compton-scattering process. Although their main results will not be influenced by this issue, since highly obscured AGNs do not appear to be the dominant population in Y19, and the use of the X-ray band in their analysis also minimizes this effect (see Section 3.5.1 of Y17), it is currently unclear whether highly obscured sources follow the same trend as the general AGN population in Y19. Therefore, it is crucial to extend the aforementioned works to the highly obscured regime.

Given the fact that M_* and SFR are positively correlated through the galaxy main-sequence relation and both of them increase with increasing redshifts owning to observational bias or/and actual galaxy evolution, a simple bivariate analysis (e.g., ${L}_{{\rm{X}}}$ versus SFR or ${M}_{* }$ ) is not able to reveal the leading factor that may predominantly govern the fueling and obscuration environments of SMBH growth. To overcome this issue, we simultaneously perform multivariate linear regression and a Spearman partial correlation test using the R packages lm.fit and ppcor of AGN parameters on all three variables, ${M}_{* }$ , SFR, and z, that describe how AGN activity and obscuration depend on ${M}_{* }$ (SFR) at given SFR ( ${M}_{* }$ ) and z, thereby enabling us to break the degeneracies.

The linear regression result for our highly obscured AGN sample using ${L}_{{\rm{X}}}$ as a direct tracer of AGN activity is

$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}\,{L}_{{\rm{X}}} & = & (0.21\pm 0.06)\times \mathrm{log}\,{M}_{* }\\ & & +\ (0.11\pm 0.05)\times \mathrm{log}\,\mathrm{SFR}\\ & & +\ (0.28\pm 0.04)\times z+(40.72\pm 0.60).\end{array}\end{eqnarray} \tag{ 3 }$

There is a statistically significant positive correlation between ${L}_{{\rm{X}}}$ and ${M}_{* }$ (3.7σ),²³ and a positive but less significant correlation between ${L}_{{\rm{X}}}$ and SFR (2.4σ). The regression result is further confirmed by the Spearman partial correlation test that ${L}_{{\rm{X}}}$ has a stronger correlation with ${M}_{* }$ (ρ = 0.29 at the 5.1σ confidence level) than with SFR (ρ = 0.13 at the 2.4σ confidence level), which appears to be in support of the scenario proposed by Y19 that $\overline{\mathrm{BHAR}}$ (equivalent to ${L}_{{\rm{X}}}$ in our analysis) is mainly linked with ${M}_{* }$ instead of SFR for non-bulge-dominated galaxies (see Section 5.3 for the result that 71% of the sources of our morphology sample have a significant disk or irregular morphology). Such enhanced AGN activity in massive galaxies may possibly be related to the greater gravitational potential, which makes it easier to fuel the central SMBH with gas in the vicinity of the nuclear region. Furthermore, we examine whether ${L}_{{\rm{X}}}$ for our bulge-dominated galaxies (i.e., the 51 SPHs identified in Section 5.3) traces SFR as proposed in Y19. This time, we do not find any statistically robust correlation between ${L}_{{\rm{X}}}$ and ${M}_{* }$ (1.0σ) or SFR (1.9σ) using Spearman partial correlation tests, perhaps because the small sample size does not properly average over all galaxies to assess duty cycle effects and thus precludes us from finding any significant trend.

Note that, even if ${M}_{* }$ is controlled, there still remains a somewhat weaker positive trend between ${L}_{{\rm{X}}}$ and SFR, in agreement with previous studies that reported a positive correlation between ${L}_{{\rm{X}}}$ and SFR (e.g., Dai et al. 2018; Lanzuisi et al. 2018) and higher average X-ray luminosities in starburst galaxies (e.g., Rodighiero et al. 2015; Grimmett et al. 2019), which could be explained by a common cold gas supply for both SF and SMBH accretion.

We also look for trends indicating whether AGN obscuration is correlated with host-galaxy properties. The linear regression result for ${N}_{{\rm{H}}}$ is

$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}\,{N}_{{\rm{H}}} & = & (-0.04\pm 0.05)\times \mathrm{log}\,{M}_{* }\\ & & +\ (0.05\pm 0.04)\times \mathrm{log}\,\mathrm{SFR}\\ & & -\ (0.02\pm 0.03)\times z+(24.2\pm 0.5),\end{array}\end{eqnarray} \tag{ 4 }$

and the Spearman correlation test yields correlation coefficients consistent with zero for both ${M}_{* }$ (ρ = −0.04) and SFR (ρ = 0.06). The lack of correlation is also confirmed by calculating the average ${N}_{{\rm{H}}}$ in different M_* or SFR bins, because in either case, we find a flat trend within 1σ uncertainties (e.g., Lutz et al. 2010; Shao et al. 2010; Rosario et al. 2012). The independence of LOS obscuration with regard to host-galaxy properties is naturally expected in the AGN unification model (Antonucci 1993), suggesting that the high ${N}_{{\rm{H}}}$ values for a significant fraction of our sources likely result from high inclination angles—instead of being caused by an intensively dusty environment as a consequence of violent mergers, where an enhancement in SFR is expected when the absorbing column density reaches the highly obscured regime (e.g., Hopkins et al. 2008). The lack of a significant correlation with total stellar mass suggests that, for the bulk of X-ray selected highly obscured AGNs, the absorbing materials responsible for the CT-level obscuration are probably confined in the nuclear region (e.g., Buchner & Bauer 2017).

This finding appears to be inconsistent with some studies that reported a positive correlation between obscuration and M_*, at least at low column densities (Buchner & Bauer 2017; Lanzuisi et al. 2018; Zou et al. 2019). To alleviate the limitation of the narrow ${N}_{{\rm{H}}}$ range being explored, we also include less obscured AGNs while performing partial correlation tests. However, no correlation is detected for either M_* or SFR when the full ${N}_{{\rm{H}}}$ range is considered. This discrepancy is likely due to the different natures of X-ray and optical obscuration (e.g., Shimizu et al. 2018; Xu et al. 2020) as well as sample differences. It is also possible that the narrow M_* range (∼10¹⁰–10¹¹ ${M}_{\odot }$ ) of our sample prevents us from finding any significant trend, thus wider surveys with similar depths are required to probe highly obscured AGNs in more massive galaxies and extend our current analyses.

5.3. Are Highly Obscured AGNs Mainly Triggered by Mergers?

In order to investigate the relevance of mergers in triggering highly obscured AGNs, we cross-match our sample with the Huertas-Company et al. (2015) catalog, which provides morphology classifications for galaxies with H-band magnitude <24.5 in the five CANDELS fields based on high-resolution HST images and deep-learning techniques. The classification algorithm is trained on the GOODS-S visual classification results (Kartaltepe et al. 2015) and has a very high accuracy. Note that all galaxies in Huertas-Company et al. (2015) are classified based on H-band images, thus we are investigating the rest-frame NIR images for low-redshift sources and rest-frame optical images for high-redshift sources. However, since Kartaltepe et al. (2015) showed that only a small fraction of their sources (i.e., 84 out of 7634 galaxies) have very different classifications between V-band and H-band images, we argue that this morphological k-correction will not significantly influence our results.

Since most of our z > 3 sources do not have measured morphology information in Huertas-Company et al. (2015), we exclude them from the morphology analysis. We use the CANDELS counterpart coordinates for our highly obscured AGNs given by the X-ray source catalogs to perform cross-matching. A total of 226 sources are matched using a 0 farcs 5 matching radius (hereafter the morphology sample). For each galaxy, five parameters are assigned to describe their morphology: f_sph, f_disk, f_irr, f_ps, and f_unc, which represent the respective possibilities that a galaxy is spheroidal, disky, irregular, point-like, and unclassifiable. Among the morphology sample of 226 sources, 191 have a set of the above morphology parameters being derived (hereafter the f-measured sources) and we divide them into four groups (Huertas-Company et al. 2015):

1.
pure bulges (SPH): f_sph > 2/3, f_disk < 2/3 and f_irr < 1/10;
2.
disks (DISK): f_disk > 2/3 and f_irr < 1/10;
3.
irregular disks (DISKIRR): f_disk > 2/3, f_sph < 2/3 and f_irr > 1/10; and
4.
irregulars/mergers (IRR): f_disk < 2/3, f_sph < 2/3 and f_irr > 1/10.

Motivated by Kocevski et al. (2015) (see their Section 3), we do not distinguish late-type and early-type disks (i.e., we merge the DISK and DISKSPH groups in Huertas-Company et al. (2015) into DISK) to reduce the possible contamination from the AGN to the bulge component. This consideration is also motivated by the fact that the disk components are easily destroyed in a major merger event (Hopkins et al. 2009); therefore, as long as a significant undisturbed disk is observable, it is less possible that the galaxies have experienced violent mergers.

Figure 10 presents the distributions of morphology type for our morphology sample in four redshift bins. Among the 191 f-measured galaxies, 173 (91%) of them have been classified as one of the four types, including 51 SPHs, 76 DISKs, 30 DISKIRRs, and 16 IRRs. The $\overline{z}$ and $\overline{\mathrm{log}\,{L}_{{\rm{X}}}}$ of the classified sources are 1.56 and 43.5 $\mathrm{erg}\ {{\rm{s}}}^{-1}$ , respectively. Most of these 173 sources (61%) have a significant disk component (f_disk > 2/3, i.e., DISK + DISKIRR; see also Schawinski et al. (2012)), while only 27% of them exhibit irregular signatures (i.e., DISKIRR + IRR).

For the 18 unclassified sources, 61% of them (11/18) have f_irr > 0.1, with $\overline{{f}_{\mathrm{irr}}}=0.31$ , $\overline{{f}_{\mathrm{disk}}}=0.48$ , and $\overline{{f}_{\mathrm{sph}}}=0.76$ (the remaining f_irr ≤ 0.1 sources have $\overline{{f}_{\mathrm{irr}}}=0.05$ , $\overline{{f}_{\mathrm{disk}}}=0.47$ , and $\overline{{f}_{\mathrm{sph}}}=0.40$ ). Visual inspection of their images confirms that some of them do exhibit irregular morphologies, and thus it is possible that these galaxies are experiencing mergers and are transforming their morphology from being disk-dominated to bulge-dominated. If we simply treat all the 11 unclassified f_irr > 0.1 sources as IRRs, then 57 out of the 191 (30%) sources show some level of irregularities (i.e., DISKIRR+IRR). This optimistic fraction is in general agreement with Kocevski et al. (2015), who proposed that ∼22% of X-ray-selected highly obscured AGNs exhibit merger or interaction signatures.

Such a small irregular fraction suggests that major mergers cannot be the leading mechanism that triggers highly obscured SMBH accretion, especially considering the fact that an irregular disk morphology does not necessarily imply galaxy interactions.²⁴ The majority (61%) of the classified sources having a significant disk component (i.e., 76 DISKs + 30 DISKIRRs) can be considered as a further argument against the major-merger scenario, which disfavors the probability that the likely time lag between the merger and the later onset of nuclear activity (e.g., Emonts et al. 2006; Hopkins 2012) prevents us from finding merger signatures that have already faded, as it is not very likely for the disk structure to survive after experiencing the violent merger process for ordinary galaxies (e.g., Hopkins et al. 2009).

However, it has been argued that major mergers may only be important in triggering the most luminous AGNs and/or the most obscured CT AGNs (e.g., Treister et al. 2012; Kocevski et al. 2015; Chang et al. 2017a). To examine such arguments, we perform a similar morphological analysis on less obscured AGNs and show the DISKIRR and IRR fractions in different ${N}_{{\rm{H}}}$ and ${L}_{{\rm{X}}}$ bins in Table 3. The 1σ errors are calculated using the method of Cameron (2011). We confirm the trend presented in Kocevski et al. (2015) that the irregular fraction increases with ${N}_{{\rm{H}}}$ . However, we note that the average ${L}_{{\rm{X}}}$ is 0.8 dex higher for the CT population, and an elevated irregular fraction in the high- ${L}_{{\rm{X}}}$ bin is also detected. To make a fair comparison, we construct ${L}_{{\rm{X}}}$ - and z-controlled CT and CN AGN subsamples using the same method described in Section 5.1, and the result is plotted in the left panel of Figure 11. It can be seen that, although the sample is very limited (only 29 CT AGNs have ${L}_{{\rm{X}}}$ - and z-matched CN counterparts with $\overline{\mathrm{log}\,{L}_{{\rm{X}}}}=43.7$ , thus we do not further divide the CN subsample into different ${N}_{{\rm{H}}}$ bins), both DISKIRR and IRR fractions are higher in CT AGNs than their low- ${N}_{{\rm{H}}}$ counterparts, suggesting that the observed difference in ${N}_{{\rm{H}}}$ should not be caused only by inclination effects. The elevated irregular fraction in the most heavily obscured CT population implies that mergers indeed play a part in boosting the nuclear obscuration (e.g., Ricci et al. 2017a; Koss et al. 2018).

**Figure 11.** Comparison of morphology classification distributions between ${L}_{{\rm{X}}}$ - and z-matched CT and CN AGNs (left), between ${M}_{* }$ - and z-matched low- ${L}_{{\rm{X}}}$ (i.e., ${L}_{{\rm{X}}}$ smaller than the median X-ray luminosity value at each redshift grid) highly obscured AGNs and nonactive galaxies (middle), and between ${M}_{* }$ - and z-matched high- ${L}_{{\rm{X}}}$ (i.e., ${L}_{{\rm{X}}}$ larger than the median X-ray luminosity value at each redshift grid) highly obscured AGNs and nonactive galaxies (right), respectively. The 1σ error bars are derived from the 16th and 84th percentiles of the 1000 resampled distributions.
Download figure:
Standard image High-resolution image

Table 3. DISKIRR and IRR Fractions for the Combined Highly Obscured and Less Obscured AGN Samples in Different ${N}_{{\rm{H}}}$ and ${L}_{{\rm{X}}}$ Bins

Fraction	Total	${N}_{{\rm{H}}}\lt {10}^{23}\ {\mathrm{cm}}^{-2}$	${10}^{23}\ {\mathrm{cm}}^{-2}\lt {N}_{{\rm{H}}}\lt {10}^{24}\ {\mathrm{cm}}^{-2}$	${N}_{{\rm{H}}}\gt {10}^{24}\ {\mathrm{cm}}^{-2}$	${L}_{{\rm{X}}}\lt {10}^{44}\ \mathrm{erg}\ {{\rm{s}}}^{-1}$	${L}_{{\rm{X}}}\gt {10}^{44}\ \mathrm{erg}\ {{\rm{s}}}^{-1}$
f_DISKIRR	${15}_{-2}^{+2} \%$	${12}_{-2}^{+3} \%$	${16}_{-3}^{+4} \%$	${22}_{-5}^{+7} \%$	${13}_{-2}^{+2} \%$	${26}_{-5}^{+7}$ %
f_IRR	${11}_{-1}^{+2} \%$	${8}_{-2}^{+3} \%$	${8}_{-2}^{+3} \%$	${12}_{-3}^{+6} \%$	${10}_{-2}^{+2} \%$	${11}_{-3}^{+6} \%$

Download table as: ASCII Typeset image

However, we notice that, in either population (high ${N}_{{\rm{H}}}$ or high ${L}_{{\rm{X}}}$ ), the DISKIRR and IRR classes still only occupy a small fraction of the total sample. To better understand the importance of mergers in triggering highly obscured AGNs in the context of galaxy evolution, we split our morphology sample into two subsamples based on the median ${L}_{{\rm{X}}}$ in each redshift grid and construct the ${M}_{* }$ - and z-controlled normal galaxy sample for each population from the S15 GOODS-S data set. Their morphology classification results are then compared with our highly obscured AGNs using the same classification criteria. The middle and right panels of Figure 11 show the average morphology classification distributions for the 1000 randomly selected control-galaxy samples and our highly obscured AGNs in the low- ${L}_{{\rm{X}}}$ and high- ${L}_{{\rm{X}}}$ regimes, respectively. For the low-luminosity highly obscured population, the irregular fraction is indistinguishable from that of control galaxies in terms of both DISKIRR fraction ( ${20}_{-3}^{+3} \%$ versus ${16}_{-3}^{+3} \%$ ) and IRR fraction ( ${12}_{-3}^{+3} \%$ versus ${13}_{-3}^{+3} \%$ ), while for the high-luminosity highly obscured population, the IRR fraction ( ${8}_{-2}^{+2} \%$ ) is similar to that of control galaxies ( ${9}_{-3}^{+3} \%$ ) but the DISKIRR fraction increases from ${12}_{-3}^{+3} \%$ for control galaxies to ${19}_{-3}^{+3} \%$ for luminous highly obscured AGNs. This, together with the result that the irregular fraction increases with ${N}_{{\rm{H}}}$ , suggests that galaxy interactions are more relevant in triggering the most luminous and the most heavily obscured (X-ray-selected) AGNs; however, it may still play a limited role, as even for such extreme populations, the disturbed hosts are still in the minority.

Therefore, although mergers do have the ability to trigger highly obscured AGNs (e.g., Ricci et al. 2017a; Goulding et al. 2018; De Rosa et al. 2018; Pfeifle et al. 2019), our result here argues that the statement that the majority of X-ray-selected highly obscured AGNs are triggered by mergers is not true. The large fractions of undisturbed disk and spheroid hosts for our highly obscured sample shown in Figure 11 implies that secular processes (e.g., galactic disk instabilities) should be the predominant triggering mechanism of the current highly obscured SMBH accretion by fueling cold gas streams to the central AGNs stochastically (e.g., Kocevski et al. 2012, 2015; Schawinski et al. 2012; Chang et al. 2017b).

5.4. Are Highly Obscured AGNs Experiencing a Blow-out Phase?

Blueshifted absorption lines detected in AGNs have been widely interpreted as an indicator of outflowing materials (e.g., Gibson et al. 2009; Filiz et al. 2013; Li et al. 2019a), which have long served as crucial ingredients in the evolutionary models that are responsible for transforming AGN types (i.e., from type 2 to type 1 or X-ray obscured to unobscured) and quenching star formation; see, e.g., King & Pounds (2015) for a review. A possible way to study whether our highly obscured AGNs are experiencing a blow-out phase and may eventually become unobscured AGNs is to investigate whether they are in the "forbidden" outflow region in the ${N}_{{\rm{H}}}$ versus ${\lambda }_{\mathrm{Edd}}$ diagram as defined in Fabian et al. (2008). To obtain ${\lambda }_{\mathrm{Edd}}$ , X-ray luminosity is converted into bolometric luminosity ( ${L}_{\mathrm{bol}}$ ) using the Lusso et al. (2012) luminosity-dependent conversion factor in the form of

$\begin{eqnarray}&&\mathrm{log}\,{L}_{\mathrm{bol}}\,/\,{L}_{{\rm{X}}}=0.230\,x+0.050\,{x}^{2}+0.001\,{x}^{3}+1.256,\end{eqnarray} \tag{ 5 }$

where x = log L_bol − 12. The black hole mass is obtained by scaling the total stellar mass using the relation proposed by Sun et al. (2015) that is parameterized as log M_BH/ ${M}_{* }$ = −2.85. The ${M}_{\mathrm{BH}}$ calculated from this method will inevitably suffer large uncertainties, but we note that it is still possible to obtain the average ${M}_{\mathrm{BH}}$ information for our sample, as various studies have observed a positive correlation between ${M}_{\mathrm{BH}}$ and ${M}_{* }$ , albeit with large dispersions (e.g., Merloni et al. 2010; Sun et al. 2015; Reines & Volonteri 2015). The ${\lambda }_{\mathrm{Edd}}$ is then calculated as ${L}_{\mathrm{bol}}/{L}_{\mathrm{Edd}}$ , where ${L}_{\mathrm{Edd}}=1.26\,\times \,{10}^{38}\ ({M}_{\mathrm{BH}}\,/\,{M}_{\odot })$ erg s⁻¹. The typical uncertainty in ${\lambda }_{\mathrm{Edd}}$ is ∼0.8 dex, which is estimated by combining the typical uncertainties (i.e., the median value of the whole sample) from the ${L}_{{\rm{X}}}$ (∼0.25 dex) and the M_* (∼0.05 dex) measurements, as well as the scatters of the ${L}_{{\rm{X}}}$ − ${L}_{\mathrm{bol}}$ (∼0.2 dex; Lusso et al. 2012) and ${M}_{\mathrm{BH}}$ − M_* (∼0.3 dex; Sun et al. 2015) relations.

The ${N}_{{\rm{H}}}$ versus ${\lambda }_{\mathrm{Edd}}$ relation is plotted in Figure 12, where the distribution of ${\lambda }_{\mathrm{Edd}}$ is also shown. It can be seen that the bulk of our sources are now in the long-lived region (e.g., Raimundo et al. 2010; Ricci et al. 2017b). As pointed out by Liu & Zhang (2011), for an effective Eddington ratio $\lambda =A{\lambda }_{\mathrm{Edd}}$ below a critical value of 7/18, where A is the boost factor of radiation pressure owing to the presence of dust, the gravitational force will always defeat radiation pressure in all directions, hence the dusty materials will not be blown out. Adopting the boost factor A calculated by Fabian et al. (2008), we derive the critical ${\lambda }_{\mathrm{Edd}}$ at which a radiatively driven outflow can be launched as a function of ${N}_{{\rm{H}}}$ . We display the results in Figure 12. The majority of our sources are located in the left side of the critical curves, suggesting that their current obscuring materials are long-lived in all directions, given their instantaneous ${\lambda }_{\mathrm{Edd}}$ . Given that these sources have smaller $\overline{\mathrm{log}\,{\lambda }_{\mathrm{Edd}}}$ (−1.5), $\overline{\mathrm{log}\,\mathrm{sSFR}}$ (−9.6 M_⊙/yr) and larger $\overline{\mathrm{log}\,{M}_{* }}$ (10.7 ${M}_{\odot }$ ) than those sources on the right side of the critical curves with the corresponding values of −0.2, −8.6/yr and 10.2 ${M}_{\odot }$ , respectively, they may have already evolved to a less active phase in terms of both star formation and BHAR, thus their obscuring materials are likely to survive from the feedback events and they will remain X-ray-obscured, instead of transforming to X-ray-unobscured AGNs. However, simulations have shown that ${\lambda }_{\mathrm{Edd}}$ can vary by orders of magnitude on timescales of ∼10⁵–10⁶ yrs (e.g., Novak et al. 2011; Yuan et al. 2018). Therefore, we cannot rule out the possibility that, once a new cycle of significant SMBH accretion is triggered (e.g., through mergers), these obscuring materials may be cleared out and an unobscured AGN may be revealed along our LOS.

6. Conclusions and Discussion

In this work, by using a large sample of X-ray-selected ${N}_{{\rm{H}}}\gt {10}^{23}\ {\mathrm{cm}}^{-2}$ AGNs (294 sources at z = 0–5) and the wealth of spectroscopic and photometric data available in CDFs, we have carried out a systematic multiwavelength study of highly obscured AGNs that is supplementary to our previous X-ray spectral and long-term variability analyses (see Paper I), aiming at examining whether highly obscured AGNs are the missing link in the merger-triggered SMBH–galaxy coevolutionary models.

Specifically, we investigate the distributions of our X-ray sources on various optical/IR/X-ray color–color diagrams in order to explore the AGN obscuration properties. We also perform detailed multicomponent SED decomposition using the X-CIGALE code (Yang et al. 2020) in order to obtain crucial host-galaxy parameters (e.g., ${M}_{* }$ and SFR). The inclusion of X-ray data and the use of the clumpy torus model in SED fitting allow us to better constrain the AGN power and thereby obtain more reliable host-galaxy properties. We explore potential correlations among M_* and SFR with the direct tracers of AGN radiation power and obscuration (i.e., the X-ray luminosity and column density derived from X-ray spectral fitting) in order to search for possible connections between the growths of SMBHs and their host galaxies. The morphology classification is then performed for the purpose of evaluating the importance of mergers in triggering highly obscured SMBH accretion. Finally, we present our analysis of whether highly obscured AGNs are sweeping out the surrounding obscuring materials, which may ultimately make them unobscured AGNs. The primary conclusions emerging from this work are summarized as follows.

1.
The Donley et al. (2012) IRAC color–color diagram can successfully identify a substantial fraction of highly obscured AGNs at $\mathrm{log}\,{L}_{{\rm{X}}}(\mathrm{erg}\ {{\rm{s}}}^{-1})\gt 44.5$ . However, the identification fraction dramatically drops to ≲20% even for sources with $44.0\lt \mathrm{log}\,{L}_{{\rm{X}}}(\mathrm{erg}\ {{\rm{s}}}^{-1})\lt 44.5$ . The low identification efficiency for X-ray luminous ( $\mathrm{log}\,{L}_{{\rm{X}}}\gt 44.0\ \mathrm{erg}\ {{\rm{s}}}^{-1}$ ) highly obscured AGNs is likely due to fact that the IRAC-missed sources have lower dust contents and/or torus CFs, and are therefore intrinsically fainter in MIR ( $\overline{\mathrm{log}\,{L}_{6\mu {\rm{m}}}}=44.1\ \mathrm{erg}\ {{\rm{s}}}^{-1}$ ) than the IRAC-selected ones ( $\overline{\mathrm{log}\,{L}_{6\mu {\rm{m}}}}=45.1\ \mathrm{erg}\ {{\rm{s}}}^{-1}$ ), given that the average X-ray luminosities for the two populations are the same ( $\overline{\mathrm{log}\,{L}_{{\rm{X}}}}=44.4\ \mathrm{erg}\ {{\rm{s}}}^{-1};$ see Section 4.1).
2.
A large fraction (84%) of our X-ray highly obscured AGNs will not be selected as heavily obscured candidates using the flux ratio of ${f}_{24\mu {\rm{m}}}/{f}_{R}\gt 1000$ and R − K > 4.5 criteria proposed by Fiore et al. (2008), even for the most luminous (i.e., ${L}_{{\rm{X}}}\gt {10}^{44}\,\mathrm{erg}\ {{\rm{s}}}^{-1}$ ) population (i.e., 70% being missed). This result suggests that the heaviest X-ray obscuration is not equivalent to the extremely large MIR-to-optical flux ratios and the reddest colors, possibly owing to the diverse structures of obscuring materials (e.g., different CFs, gas/dust contents, and ${N}_{{\rm{H}}}$ distribution), complex origins of the LOS X-ray obscuration (e.g., dust-free BLR gas, dusty torus, disk wind, and/or ISM) as well as galaxy contamination of the observed colors (see Section 4.2).
3.
The 50th percentile of the SFR distribution for highly obscured AGN hosts is similar to that of ${M}_{* }$ - and z-controlled normal galaxies ( $\overline{{\rm{\Delta }}\mathrm{log}\,{\mathrm{SFR}}_{50\mathrm{th},\mathrm{agn}-\mathrm{galaxy}}}\,=-{0.03}_{-0.12}^{+0.12}$ ), but the 20th percentile for the former is characteristically larger than that of the latter ( $\overline{{\rm{\Delta }}\mathrm{log}\,{\mathrm{SFR}}_{20\mathrm{th},\mathrm{agn}-\mathrm{galaxy}}}={0.83}_{-0.20}^{+0.17}$ ). As a result, the SFR distribution of highly obscured AGNs is narrower than normal galaxies (i.e., there is a lack of quiescent hosts), suggesting that the sufficient cold gas supply that is responsible for maintaining star formation may also be an important source in fueling highly obscured SMBH accretion. Furthermore, the 80th percentile of the SFR distribution of highly obscured AGNs is not enhanced relative to control galaxies ( $\overline{{\rm{\Delta }}\mathrm{log}\,{\mathrm{SFR}}_{80\mathrm{th},\mathrm{agn}-\mathrm{galaxy}}}\,=-{0.30}_{-0.07}^{+0.08}$ ), suggesting that the presence of X-ray-selected highly obscured AGNs is not more frequently connected to starburst activities (see Section 5.1).
4.
The multivariate linear regression and Spearman partial correlation analyses among ${L}_{{\rm{X}}}$ , ${M}_{* }$ , SFR, and z suggest that highly obscured SMBH accretion (traced by ${L}_{{\rm{X}}}$ ) is more fundamentally related to ${M}_{* }$ , which might be a result of the greater gravitational potentials of massive hosts, being consistent with previous studies on non-bulge-dominated galaxies (see conclusion 6). However, highly obscured SMBH accretion still remains a positive trend with SFR after M_* is controlled, albeit at a weaker significance level, suggesting that, at a given ${M}_{* }$ , galaxies with sufficient gas contents that are able to fuel higher SFRs are also more likely to trigger highly obscured AGNs (see Section 5.2).
5.
We find no correlation between ${N}_{{\rm{H}}}$ and either ${M}_{* }$ or SFR, consistent with the prediction from the AGN unification model (Antonucci 1993). Such a result suggests that, for a significant fraction of our sources, their high ${N}_{{\rm{H}}}$ values likely result from high inclination angles, instead of being caused by an intensively dusty environment as a consequence of violent mergers where an enhancement in SFR is expected when the absorbing column densities reach the highly obscured regime (e.g., Hopkins et al. 2008). The lack of a significant correlation with total stellar mass suggests that for the bulk of X-ray selected highly obscured AGNs, the absorbing materials responsible for the CT-level obscuration are probably confined in the nuclear region.
6.
To examine whether highly obscured SMBH accretion is mainly triggered by mergers, we cross-match our sample with the Huertas-Company et al. (2015) galaxy morphology catalog. We find that 61% of them have a significant disk component, while only 27% of them exhibit irregular signatures (9% IRR + 18% DISKIRR). The incidence of disturbed morphologies increasing with both ${L}_{{\rm{X}}}$ and ${N}_{{\rm{H}}}$ , which supports the scenario that mergers are more relevant in triggering the most luminous and the most heavily obscured (X-ray-selected) AGNs. However, mergers may still only play a limited role, as even for such extreme populations, the disturbed hosts are still in the minority (see Section 5.3).
7.
The majority of our sources are located in the stable long-lived region in the ${N}_{{\rm{H}}}$ versus ${\lambda }_{\mathrm{Edd}}$ plane defined in Fabian et al. (2008). The fact that long-lived sources have, on average, smaller $\overline{{\lambda }_{\mathrm{Edd}}}$ and $\overline{\mathrm{sSFR}}$ and higher $\overline{{M}_{* }}$ with respect to sources that have matched $\overline{{N}_{{\rm{H}},\mathrm{LOS}}}$ but lie in or close to the outflow region suggests that they may have already evolved to a less active phase in terms of both star formation and BHAR, thus their obscuring materials are likely to survive from the feedback events and they will remain X-ray-obscured. However, we cannot rule out the possibility that, once a new cycle of significant SMBH accretion activity is being triggered, these obscuring materials may be cleared out and an unobscured AGN will be revealed along our LOS (see Section 5.4).

Overall, our findings suggest that the majority of the X-ray-selected highly obscured AGNs within the luminosity and redshift ranges examined here are not the "missing link" (i.e., dust-enshrouded AGNs likely having enhanced star-forming activities, disturbed host-galaxy morphologies, and strong outflows that may eventually make themselves unobscured) in the merger-triggered SMBH–galaxy coevolution and AGN type-transition models, given the combined evidence of complex origins of their heavy LOS X-ray obscuration, the star formation activity being similar to that of normal star-forming galaxies, the lack of correlation between absorbing column densities and SFR, the large fraction of sources having undisturbed hosts, and the fact that the majority of them are far away from the outflow region defined in the ${N}_{{\rm{H}}}$ versus ${\lambda }_{\mathrm{Edd}}$ plane.

However, although our work has provided deep insights into the elusive highly obscured AGN population by taking advantage of the deepest Chandra surveys, the current X-ray data are still insufficient. The small sky coverages of CDFs severely limit the source statistics toward the highest-luminosity end, and the current detection limit for Chandra is insufficient to probe CT AGNs that have ${N}_{{\rm{H}}}$ as high as 10²⁵ ${\mathrm{cm}}^{-2}$ . Therefore, we are still missing the most luminous and/or the most heavily obscured AGNs, which may be most important to understand the merger models (see Section 5.3). To improve the current situation, wider X-ray surveys with similar or deeper depth are necessary. This requirement is expected to be fulfilled by future Athena (Nandra et al. 2013) and Lynx (Gaskin et al. 2019) surveys, which can reach detection limits close to (i.e., Athena) or even deeper (i.e., Lynx) than the current Chandra deep surveys while also achieving much larger survey areas. In concert with the improved constraints of dust properties and star formation and AGN activities at IR wavelengths by future MIR spectroscopy from, e.g., the James Webb Space Telescope (e.g., Kirkpatrick et al. 2017) and the SPICA mission (e.g., Roelfsema et al. 2018), we will be able to better address the issues related to current studies and have a more profound understanding of the highly obscured AGN population as well as their role in galaxy evolution.

We thank the referee for helpful suggestions. J. Y. L. and Y. Q. X. acknowledge support from the National Natural Science Foundation of China (NSFC-11890693, 11421303), the CAS Frontier Science Key Research Program (QYZDJ-SSW-SLH006), and the K. C. Wong Education Foundation. M. Y. S acknowledges support from the National Natural Science Foundation of China (NSFC-11973002). W. N. B. acknowledges support from Chandra X-ray Center grant GO8-19076X, NASA grant 80NSSC19K0961, and Penn State ACIS Instrument Team Contract SV4-74018 (issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060). The Chandra ACIS team Guaranteed Time Observations (GTO) utilized were selected by the ACIS Instrument Principal Investigator, Gordon P. Garmire, currently of the Huntingdon Institute for X-ray Astronomy, LLC, which is under contract to the Smithsonian Astrophysical Observatory via Contract SV2-82024. F. V. acknowledges financial support from CONICYT and CASSACA through the Fourth Call for Tenders of the CAS-CONICYT Fund, and CONICYT grant Basal-CATA AFB-170002. P. T. acknowledges financial contribution from the agreement ASI-INAF n.2017-14-H.0. L. L. F. acknowledges support from the National Natural Science Foundation of China (NSFC-11822303 and 11773020) and Shandong Provincial Natural Science Foundation, China (ZR2017QA001, JQ201801).

Piercing through Highly Obscured and Compton-thick AGNs in the Chandra Deep Fields. II. Are Highly Obscured AGNs the Missing Link in the Merger-triggered AGN–Galaxy Coevolution Models?

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Multiwavelength Data

2.1. X-Ray Data

2.2. UV and Optical Data

2.3. NIR and MIR Data

2.4. FIR Data

2.5. Construction of Broadband SEDs

3. SED Fitting Method and Results

4. Revisiting the Optical/IR-selection Methods

4.1. Can IRAC Colors Effectively Identify Luminous Highly Obscured AGNs?

4.2. Is the High Ratio of ${f}_{24\mu {\rm{m}}}/{f}_{R}$ an Efficient Method to Select Highly Obscured AGNs?

4.3. Can the Value of ${L}_{6\mu {\rm{m}}}/{L}_{{\rm{X}},\mathrm{obs}}$ Be Used As a Reliable Indicator of ${N}_{{\rm{H}}}$ ?

5. Are Highly Obscured AGNs the Missing Link in the Merger Models?

5.1. Do Highly Obscured AGNs Show Enhanced Star Formation Activities?

5.2. Are AGN Activity and Obscuration Linked with Host-galaxy Properties in Highly Obscured AGNs?

5.3. Are Highly Obscured AGNs Mainly Triggered by Mergers?

5.4. Are Highly Obscured AGNs Experiencing a Blow-out Phase?

6. Conclusions and Discussion

Footnotes

Piercing through Highly Obscured and Compton-thick AGNs in the Chandra Deep Fields. II. Are Highly Obscured AGNs the Missing Link in the Merger-triggered AGN–Galaxy Coevolution Models?

Article metrics

Permissions

Share this article

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Multiwavelength Data

2.1. X-Ray Data

2.2. UV and Optical Data

2.3. NIR and MIR Data

2.4. FIR Data

2.5. Construction of Broadband SEDs

3. SED Fitting Method and Results

4. Revisiting the Optical/IR-selection Methods

4.1. Can IRAC Colors Effectively Identify Luminous Highly Obscured AGNs?

4.2. Is the High Ratio of {f}_{24\mu {\rm{m}}}/{f}_{R} an Efficient Method to Select Highly Obscured AGNs?

4.3. Can the Value of {L}_{6\mu {\rm{m}}}/{L}_{{\rm{X}},\mathrm{obs}} Be Used As a Reliable Indicator of {N}_{{\rm{H}}}?

5. Are Highly Obscured AGNs the Missing Link in the Merger Models?

5.1. Do Highly Obscured AGNs Show Enhanced Star Formation Activities?

5.2. Are AGN Activity and Obscuration Linked with Host-galaxy Properties in Highly Obscured AGNs?

5.3. Are Highly Obscured AGNs Mainly Triggered by Mergers?

5.4. Are Highly Obscured AGNs Experiencing a Blow-out Phase?

6. Conclusions and Discussion

Footnotes

4.2. Is the High Ratio of ${f}_{24\mu {\rm{m}}}/{f}_{R}$ an Efficient Method to Select Highly Obscured AGNs?

4.3. Can the Value of ${L}_{6\mu {\rm{m}}}/{L}_{{\rm{X}},\mathrm{obs}}$ Be Used As a Reliable Indicator of ${N}_{{\rm{H}}}$ ?