Piercing Through Highly Obscured and Compton-thick AGNs in the Chandra Deep Fields: I. X-ray Spectral and Long-term Variability Analyses
Junyao Li, Yongquan Xue, Mouyuan Sun, Teng Liu, Fabio Vito, William N. Brandt, Thomas M. Hughes, Guang Yang, Paolo Tozzi, Shifu Zhu, Xuechen Zheng, Bin Luo, Chien-Ting Chen, Cristian Vignali, Roberto Gilli, Xinwen Shu
DDraft version April 9, 2019
Typeset using L A TEX twocolumn style in AASTeX62
Piercing Through Highly Obscured and Compton-thick AGNs in the
Chandra
Deep Fields: I. X-ray Spectral andLong-term Variability Analyses
Junyao Li,
1, 2
Yongquan Xue,
1, 2
Mouyuan Sun,
1, 2
Teng Liu, Fabio Vito,
4, 5
William N. Brandt,
6, 7, 8
Thomas M. Hughes,
1, 2, 9, 10
Guang Yang,
6, 7
Paolo Tozzi, Shifu Zhu,
6, 7
Xuechen Zheng,
1, 2, 12
Bin Luo,
13, 14, 15
Chien-Ting Chen, Cristian Vignali,
17, 18
Roberto Gilli, and Xinwen Shu CAS Key Laboratory for Research in Galaxies and Cosmology, Department of Astronomy, University of Science and Technology ofChina, Hefei 230026, China; [email protected], [email protected] School of Astronomy and Space Science, University of Science and Technology of China, Hefei 230026, China Max-Planck-Institut f¨ur extraterrestrische Physik, Giessenbachstrasse 1, D-85748 Garching bei M¨u nchen, Germany Instituto de Astrofisica and Centro de Astroingenieria, Facultad de Fisica, Pontificia Universidad Catolica de Chile, Casilla 306,Santiago 22, Chile Chinese Academy of Sciences South America Center for Astronomy, National Astronomical Observatories, CAS, Beijing 100012, China Department of Astronomy & Astrophysics, 525 Davey Lab, The Pennsylvania State University, University Park, PA 16802, USA Institute for Gravitation and the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA Department of Physics, The Pennsylvania State University, University Park, PA 16802, USA CAS South America Center for Astronomy, China-Chile Joint Center for Astronomy, Camino El Observatorio Instituto de F´ısica y Astronom´ıa, Universidad de Valpara´ıso, Avda. Gran Breta˜na 1111, Valpara´ıso, Chile Istituto Nazionale di Astrofisica (INAF) – Osservatorio Astrofisico di Firenze, Largo Enrico Fermi 5, I-50125 Firenze Italy Leiden Observatory, Leiden University, PQ Box 9513, NL-2300 RA Leiden, the Netherlands School of Astronomy and Space Science, Nanjing University, Nanjing 210093, China Key Laboratory of Modern Astronomy and Astrophysics (Nanjing University), Ministry of Education, Nanjing, Jiangsu 210093, China Collaborative Innovation Center of Modern Astronomy and Space Exploration, Nanjing 210093, China Astrophysics Office, NASA Marshall Space Flight Center, ZP12, Huntsville, AL 35812 Dipartimento di Fisica e Astronomia, Alma Mater Studiorum, Universit`a degli Studi di Bologna, Via Gobetti 93/2, I-40129 Bologna,Italy INAF – Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Gobetti 93/3, I-40129 Bologna, Italy Department of Physics, Anhui Normal University, Wuhu, Anhui, 241000, China
ABSTRACTWe present a detailed X-ray spectral analysis of 1152 AGNs selected in the
Chandra
Deep Fields(CDFs), in order to identify highly obscured AGNs ( N H > cm − ). By fitting spectra withphysical models, 436 (38%) sources with L X > erg s − are confirmed to be highly obscured,including 102 Compton-thick (CT) candidates. We propose a new hardness-ratio measure of theobscuration level which can be used to select highly obscured AGN candidates. The completenessand accuracy of applying this method to our AGNs are 88% and 80%, respectively. The observedlog N − log S relation favors cosmic X-ray background models that predict moderate (i.e., betweenoptimistic and pessimistic) CT number counts. 19% (6/31) of our highly obscured AGNs that haveoptical classifications are labeled as broad-line AGNs, suggesting that, at least for part of the AGNpopulation, the heavy X-ray obscuration is largely a line-of-sight effect, i.e., some high-column-densityclouds on various scales (but not necessarily a dust-enshrouded torus) along our sightline may obscurethe compact X-ray emitter. After correcting for several observational biases, we obtain the intrinsic N H distribution and its evolution. The CT-to-highly-obscured fraction is roughly 52% and is consistentwith no evident redshift evolution. We also perform long-term ( ≈
17 years in the observed frame)variability analyses for 31 sources with the largest number of counts available. Among them, 17sources show flux variabilities: 31% (5/17) are caused by the change of N H , 53% (9/17) are caused bythe intrinsic luminosity variability, 6% (1/17) are driven by both effects, and 2 are not classified dueto large spectral fitting errors. a r X i v : . [ a s t r o - ph . GA ] A p r Li, Xue, Sun et al.
Keywords: galaxies: active — galaxies: nuclei — quasars: general — X-rays: galaxies INTRODUCTIONHighly obscured active galactic nuclei (AGNs), whichare defined as AGNs with hydrogen column density ( N H )larger than 10 cm − , are believed to represent a crucialphase of active galaxies. According to our knowledge ofco-evolution of supermassive black holes (SMBHs) andtheir host galaxies (for reviews, see, e.g., Alexander &Hickox 2012; Kormendy & Ho 2013) and the hierar-chical galaxy formation model, AGN activity may betriggered in a dust enshrouded environment, into whichgas inflows due to either internal (e.g., disk instabili-ties; e.g., Hopkins & Hernquist 2006) or external (e.g.,major mergers; e.g., Di Matteo et al. 2005) processesboth fuel and obscure the SMBH accretion, resultingin short-lived heavily obscured AGNs (e.g., Fiore et al.2012; Morganti 2017). The subsequent AGN feedbackprocess may blow out the obscuring material and leaveout an unobscured optically bright quasar. Comparedwith unobscured AGNs, AGNs in the highly obscuredphase tend to have smaller BH masses, higher Edding-ton ratios ( λ Edd ; e.g., Lanzuisi et al. 2015) and largermerger fractions (e.g., Kocevski et al. 2015; Ricci et al.2017a), which may indicate a fast growth state of cen-tral SMBHs (e.g., Goulding et al. 2011). Moreover, thecosmic X-ray background (CXB) synthesis models alsorequire a sizable population of highly obscured AGNs,or even Compton-thick (CT) AGNs ( N H (cid:38) cm − ;see, e.g., Comastri 2004; Xue 2017; Hickox & Alexan-der 2018 for reviews), to reproduce the peak of CXB at20 −
30 keV (e.g., Gilli et al. 2007, but see Treister et al.2009). Therefore, the study of highly obscured AGNsacross cosmic epochs is vital for our understanding ofthe AGN triggering mechanism, SMBH growth, AGNenvironment and the origin of CXB.Thanks to the powerful penetrability of high energyX-ray photons, X-ray observations provide a great win-dow to uncover the mysterious veil of these heavily ob-scured sources that are likely missed in optical surveys.In the past twenty years or so, the deep X-ray sur-veys conducted by
Chandra (e.g., Alexander et al. 2003;Xue et al. 2011, 2016; Luo et al. 2017),
XMM-Newton (e.g., Ranalli et al. 2013),
Swift /BAT (e.g., Baumgart-ner et al. 2013) and
NuSTAR (e.g., Lansbury et al. 2017)have provided relatively unbiased AGN samples thanksto their unprecedented depths and sensitivities, whichallow us to identify a significant population of heav-ily obscured AGNs (e.g., Risaliti et al. 1999; Bright-man & Ueda 2012; Ricci et al. 2015) using either X-raycolor, spectral analysis and/or stacking technique (e.g., Alexander et al. 2011; Iwasawa et al. 2012; Georgan-topoulos et al. 2013; Brightman et al. 2014; Corral et al.2014; Del Moro et al. 2016; Koss et al. 2016). Moreover,the combination of mid-infrared (MIR), optical and X-ray data provides additional methods to select heavilyobscured systems, such as MIR excess (e.g., Daddi 2007;Alexander et al. 2008; Luo et al. 2011), high 24 µm tooptical flux ratio (e.g., Fiore et al. 2009) and high X-rayto optical flux ratio (e.g., Fiore et al. 2003).Among various methods, X-ray spectroscopy pro-vides the most direct and unambiguous way to mea-sure the column density of the obscuring materials.Several previous studies have focused on deriving theintrinsic N H distribution corrected for the survey bi-ases. Tozzi et al. (2006) presented a N H distributionthat has an approximately log-normal shape peaking at ∼ cm − and with an excess at ∼ cm − . Liuet al. (2017) (hereafter L17) reported a similar N H dis-tribution that peaks in a higher N H range, due to theinclusion of more sources with low X-ray luminosities( L X ) and high redshifts than that of Tozzi et al. (2006),which are expected to have relatively high N H values.However, both works only focused on bright AGNs, andneither was dedicated to or optimized for investigatinghighly obscured sources. In particular, L17 excluded CTAGNs in their work and only focused on the Compton-thin population. Hence the absorption distribution andevolution of the most deeply buried AGNs are still un-clear, especially at high redshifts (Vito et al. 2018).Therefore, unveiling the apparently faint, CT regimeusing the deepest X-ray survey data is indispensable forfully understanding the entire AGN population.There have also been several attempts to constrain theobscured AGN fraction and CT fraction on the basis ofmodeling the CXB (e.g., Gilli et al. 2007; Akylas et al.2012), the X-ray luminosity function (e.g., Aird et al.2015; Buchner et al. 2015) or X-ray spectral analysis(e.g., Burlon et al. 2011). Although most of the studiessupport the picture of evolved absorption with redshiftand luminosity (e.g., Hasinger 2008; Liu et al. 2017),the value of the CT fraction varies from study to study,ranging from a few percent to ∼ N H cir-cumstances. Therefore, deeper X-ray observations aswell as the physically appropriate spectral models for ighly Obscured and Compton-thick AGNs in CDFs Chandra
Deep Fields(CDFs) surveys (see, e.g., Brandt & Alexander 2015 andXue 2017 for reviews), which consist of the 7 Ms
Chan-dra
Deep Field-South survey (CDF-S; Luo et al. 2017)and the 2 Ms
Chandra
Deep Field-North survey (CDF-N; Xue et al. 2016) along with the 250 ks Extended
Chandra
Deep Field-South survey (E-CDF-S; Xue et al.2016), provide us the most promising data to studyhighly obscured AGNs. In particular, the 7 Ms CDF-S,which is the deepest X-ray survey to date, significantlyimproves the count statistics that allow us to extracthigh-quality X-ray spectra, detect more faint, highly ob-scured sources and perform more robust spectral analy-ses compared with previous 4 Ms analyses (e.g., Bright-man et al. 2014). Furthermore, recent works suggestthat the power spectral density (PSD) break frequencyof AGN light curves might be related to N H variability(Zhang et al. 2017; Gonz´alez-Mart´ın 2018), which makesobscuration a very important factor in investigating thedriving mechanism of AGN variability. Benefiting fromthe very long timespan of the 7 Ms CDF-S data (16.4years in the observed frame), we are able to, for thefirst time, quantify the detailed variability behavior fora large, dedicated sample of highly obscured AGNs, inorder to better understand the location of the obscur-ing materials and their contribution to AGN variability(also see Yang et al. 2016, Y16 hereafter).In this study, we construct the largest dedicated highlyobscured AGN sample in the deepest Chandra surveyswhich enables us to extend the studies of deeply buriedsources to lower luminosities and higher redshifts ingreat details, and present systematic X-ray spectral andlong-term variability analyses to study their evolutionand physical properties. This paper is organized as fol-lows. We describe our data reduction procedure andsample selection in Section 2. In Sections 3 and 4, wepresent detailed X-ray spectral analyses of our sample,focusing on the column density and luminosity distribu-tions of highly obscured AGNs as well as their relations;the number counts of CT AGNs and the constraint onCXB models. The reprocessed components, the cover-ing factor of the obscuring materials and their indica-tions to AGN structures are discussed in Section 5. InSection 6, by correcting for several observational biases,we constrain the intrinsic N H distribution representa-tive for the highly obscured AGN population and studyits evolution across cosmic time. In Section 7, we se-lect a subsample of highly obscured AGNs that havelargest counts available and perform detailed long-termX-ray variability analyses, in order to find out the vari- log (0.5-7 keV Net Counts) N full samplehighly obscuredvariability sample0 1 2 3 4 5 z N full samplespec- z Figure 1.
Top:
Chandra net counts distributions of thefull sample of 1152 AGNs (blue solid histogram), the highlyobscured sources confirmed with X-ray spectral fitting inSection 4 (orange dashed histogram), and the subsample se-lected to perform variability analyses in Section 7 (solid greenhistogram), respectively. Bottom: Redshift distributions forthe full sample (blue histogram) and the 603 sources withspectroscopic redshifts (purple histogram), respectively. able fraction as well as the main driven mechanism ofvariability.Throughout this paper, we adopt a Galactic columndensity of N H = 8 . × cm − for the CDF-S, and N H = 1 . × cm − for the CDF-N, respectively(Stark et al. 1992). We adopt cosmological parametersof H = 70 . − Mpc − , Ω M = 0 .
30, and Ω Λ = 0 . σ confidence level unless other-wise stated. DATA REDUCTION AND SAMPLE SELECTIONThis work is based on the
Chandra data in the 7 MsCDF-S main-source catalog (Luo et al. 2017) and the2 Ms CDF-N main-source catalog (Xue et al. 2016).Since the exposure in the E-CDF-S is ∼ Li, Xue, Sun et al. tails). The redshift (spectroscopic redshift) complete-nesses in the CDF-S and the CDF-N main source cata-logs are 99.4% (64.8%) and 95.2% (52.4%), respectively;and the corresponding mean 1 σ photometric redshift er-rors are about 0.21 and 0.17, respectively.The source spectra from individual observations wereextracted using the ACIS Extract (AE) software pack-age (Broos et al. 2010). AE generates the point-spreadfunction (PSF) model based on the MARX ray-tracingsimulator and constructs a polygonal extraction re-gion that corresponds to an encircled-energy fractionof ∼ BETTER BACKGROUNDS algorithm. The most significantaspect of the above photometry and spectral extrac-tion procedure, compared to the widely used circular-aperture extraction, is that it can obtain photometryand spectra as accurate as possible and remove the con-tamination from neighboring sources to faint sources tothe greatest extent (see Section 3.2 of Xue et al. 2011 fordetails). This is extremely important for our work sincewe are dealing with the highly obscured AGNs that aregenerally fainter and expected to have limited countsdue to significant obscuration.The spectra eventually used in this work are themerged spectra for which all the individual observationswere matched to a common K s -band astrometric frame(see Section 2.2.1 in Xue et al. 2016 ) and stacked usingthe MERGE OBSERVATIONS algorithm in AE. The corre-sponding response matrix files (rmf) and ancillary re-sponse files (arf) were also generated and combined dur-ing this stage.We construct our sample by selecting sources which(1) are classified as AGN (TYPE=AGN; see Section4.5 in Luo et al. 2017 for details); (2) have 0.5–7 keVnet counts >
20 (FB COUNTS >
20) to allow basic X-ray spectral fitting; and (3) have redshift measurements(ZFINAL ! = −
1) in the main catalogs. The result-ing full sample consists of 1152 sources, with 660 fromthe CDF-S and 492 from the CDF-N, respectively. Thecounts and redshift distributions for the full sample areshown in Figure 1. The median counts are 112, and 53%(33%) of the sources have counts larger than 100 (200).The median redshift is 1.45, and 603 (52%) sources havespectroscopic redshifts. X-RAY SPECTRAL ANALYSIS3.1.
Spectral Fitting Models
We use MYTorus-based models (Murphy & Yaqoob2009) to fit the observed-frame 0.5–7 keV spectra of thefull sample in order to identify heavily obscured sources. Due to limited counts, we do not bin the spectra becauseit may lose some key information of the sources. We usethe Cash statistic (Cstat in XSPEC; Cash 1979) as ourspectral fitting statistic. Cstat has a similar probabilitydistribution to χ statistics and has been proved to bemore appropriate in the low-counts regime. Since Cstatis not appropriate for the background-subtracted spec-tra, we simultaneously fit the source and backgroundspectra, with the latter (with 1642 median counts) be-ing fit with the cplinear model (Broos et al. 2010) thatproperly describes the observed Chandra background bya continuous piecewise-linear function in ten energy seg-ments (an example is shown in Figure 3). In this way, weare able to maximize the usage of information relevantto the sources.We adopt two models to fit the source spectra withdifferent components and degrees of freedom (d.o.f) inorder to find a statistically robust best-fit one:1. The MYTorus baseline model: phabs × (MYTZ × zpow + f ref × MYTS + f ref × gsmooth (MYTL)).2. The soft-excess model: phabs × (MYTZ × zpow + f ref × MYTS+ f ref × gsmooth (MYTL)+ f exs × zpow ).These models include all the typical spectral featuresfound in highly obscured AGNs. The phabs compo-nent models the Galactic photoelectric absorption. TheMYTZ × zpow term represents the zeroth-order trans-mitted power-law continuum across the torus that takesinto account both photoelectric absorption and Comp-ton scattering processes. MYTS and gsmooth (MYTL)stand for the reflection component and the broadenedFe K α , K β fluorescent emission lines, respectively. Asecond power law is added to represent the soft excesscomponent often found in AGN spectra, possibly origi-nating from the zeroth-order continuum being scatteredby the extended Compton-thin (CN) materials in ob-scured AGNs (Guainazzi et al. 2005; Bianchi et al. 2006;Corral et al. 2011). During spectral fitting, the inclina-tion angle θ is fixed at 75 ◦ , which is the average valuewithin a range of 60 ◦ to 90 ◦ where the torus interceptsour line of sight (l.o.s). The nominal normalizations ofMYTS, MYTL and the second power law are set to bethe same as that of the first power law, and we use theconstants f ref and f exs to represent the real normaliza-tions. f exs is allowed to vary between 0 − .
1, and weassume a 200 keV high-energy cutoff throughout. Onething to mention is that the best-fit N H values in MY-Torus represent the equatorial values, and we always usethe l.o.s values N H , l . o . s = (1 − θ ) N H in this work.3.2. Model Selection ighly Obscured and Compton-thick AGNs in CDFs Energy (keV) p h o t o n s c m s k e V (a) pwl= 2.0, log N H = 10 , f exs = 5%transimittedreflectionFe K + Ksoft excess Energy (keV) (b) pwl= 2.0, log N H = 10 , f exs = 1%transimittedreflectionFe K + Ksoft excess Figure 2.
Examples of MYTorus-based models with different parameters used in spectral fitting. The unattenuated continuum(pwl), absorbed continuum (blue curve), transmitted zeroth-order continuum, reflection component, iron emission lines and thesoft excess component are shown in different colors. (a) Model used to generate curve A in Figure 7 with N H , Γ and f exs fixedat 10 cm − , 2.0 and 5%, respectively. (b) Model used to generate curve B in Figure 7 with N H , Γ and f exs fixed at 10 cm − ,2.0 and 1%, respectively. Energy (keV) C o un t s s k e V Figure 3.
An example of fitting the
Chandra backgroundspectrum that has 332 counts using the cplinear model(binned only for illustration purpose).
Several previous works suggested that the inclusion ofthe soft excess and reflection components in the spec-tral models is crucial for us to correctly estimate N H ofhighly obscured sources (e.g., Brightman et al. 2014;Lanzuisi et al. 2015). However, the low counts of manysources do not allow us to apply complex models withfree parameters since it could lead to large degenera-cies. Therefore, we choose the model components andthe parameter spaces according to the following criteria:1. We fix Γ at 1.8 (e.g., Tozzi et al. 2006; Marchesi et al.2016) for the 851 sources with counts less than 300.For the 301 sources with counts larger than 300, weset Γ free to find a best-fit value.2. For all sources, we first fit the spectra with a free f ref .If f ref is less than 10 − , we then fix it to 10 − that isan arbitrary value set as the lower limit for f ref (i.e., indicating a negligible reflection component); other-wise, we fix f ref at 1 that is the default value adoptedin MYTorus.3. To determine whether we should add a soft excesscomponent, we compare the Cstat between modelswith and without considering soft excess emission.The best-fit model is chosen to be the one that has astatistically robust low Cstat value. More specifically,adding a new soft excess component should improvethe Cstat at least for ∆ C = 3 .
84 with 1 more d.o.f.This criterion is based on the fact that ∆ C approx-imately follows the χ distribution, thus ∆ C = 3 . >
95% confidence level(Tozzi et al. 2006; Liu et al. 2017).4. In order to avoid extremely untypical Γ caused by thedegeneracy between Γ and N H , we re-fit the spectraof those sources with Γ pegged at 1.4 (i.e., the lowestvalue permitted by MYTorus) or Γ > λ Edd
AGNs;e.g., Wang et al. 2004; Fanali et al. 2013) by fixingtheir Γ at 1.8. If ∆ C = C fix − C free > f ref as model A (B) and model 1 (2) with f ref fixed at 1.0 as model C (D), as summarized in Table1. In order to ensure the consistency of the models usedfor subsequent bias corrections (see Section 6) as muchas possible, we have to assume fixed parameters in caseof low counts. We justify our models in Appendix A by Li, Xue, Sun et al.
Table 1.
Model definitions used in this workmodel f ref f exs A 10 −
0B 10 − < f exs < .
1C 1 . . < f exs < . validating that the usage of fixed parameters does notsignificantly affect our results. We also compare our fit-ting results with several previous works (see Section 4.1)and those obtained from the Borus model (Balokovi´cet al. 2018; see Section 5.3) and find good consistency.We note that the simple absorbed power-law model(e.g., phabs × zwabs × zpow ; hereafter model z ) has beenwidely used in the literature to obtain rough estimatesof AGN parameters. However, such a model is not ap-propriate for highly obscured AGNs since zwabs onlymodels photoelectric absorption and does not take intoaccount the Compton scattering process, which is par-ticularly important in the high N H regime. Therefore,we also fit the spectra using model z , aiming at directlytesting how accurately such a simple model reproducesthe main spectral parameters. SPECTRAL FITTING RESULTSOn the basis of the spectral fitting results, we findthat 39% (458/1152) of sources in the full sample areidentified to be highly obscured AGNs (hereafter thehighly obscured sample; HOS). The main fitting resultsfor HOS using MYTorus models A–D are presented inTable 2 and the full version is available on line. Thebest-fit models for sources in different count bins and N H ranges are summarized in Table 3. In the follow-ing analyses, we only consider the 436 sources that havethe absorption-corrected, rest-frame 2–10 keV luminos-ity (calculated from the lumin command after deletingall the additional components and only keeping the zpow component) L X > erg s − to avoid possible contam-ination from star-forming galaxies.4.1. Comparisons With Previous Works and the Modelof phabs × zwabs × zpow Thanks to the increased exposure time, the improveddata reduction procedure (see Table 1 in Xue et al. 2016for a summary), the updated redshift measurements andthe usage of physically appropriate spectral fitting mod-els for highly obscured AGNs, we are able to extracthigher-quality spectra and obtain more robust param-eter constraints than previous works in the same field(e.g., Tozzi et al. 2006), especially for faint sources. Be-fore performing further analyses, we first make direct comparisons with several works which presented spec-tral fitting parameters in the CDF-S to understand howdifferent methods influence the spectral fitting results.The works used for comparisons are:1. Liu et al. (2017) (L17): L17 performed an exten-sive X-ray spectral analysis for the bright sources(hard-band counts >
80) in the 7 Ms CDF-S us-ing the wabs × ( plcabs × power − law + zgauss + power − law + plcabs × pexrav × constant ) model.The background spectrum was also fitted with the cplinear model. The redshift information used inthe two works is the same.2. Brightman et al. (2014) (B14): B14 presented X-ray spectral fitting results in the 4 Ms CDF-S usingthe BNtorus model (Brightman & Nandra 2011).112 sources are found common between the twoworks; 24 of them have redshifts different by morethan 0.2 (∆ z > .
2) compared with the valuesadopted in the updated 7 Ms catalog, and are ne-glected in the following comparison.3. Buchner et al. (2015) (B15): The 4 Ms CDF-Sspectra were analyzed using a physically motivatedtorus model and a Bayesian methodology to esti-mate spectral parameters. 114 sources are foundcommon between the two works and 30 of themwith ∆ z > . N H and L X valuesare in general agreement with previous works, despitethat for individual sources, the usage of different modelconfigurations and data may result in large discrepan-cies. The largest distinction happens between B14 andour work that ten highly obscured AGNs in our sampleare reported to be unobscured in B14. To understandthe discrepancy, we re-fit their spectra using the samemodel and method as described in B14 and the resultsagain confirm their heavily obscured nature. Therefore,the large discrepancy could be due to the adopted dataof different depths as well as the different data reductionand spectral extraction methods used in the two works.The comparisons between the results obtained throughMYTorus and model z are shown in Figure 5. Notethat since MYTorus does not allow N H to vary below10 cm − , a number of unobscured sources are thushave best-fit N H = 10 cm − . Therefore, we only con-sider sources with MYTorus log N H > . − in thecomparison.The N H and the observed 0.5–7 keV flux derived bythe two models are generally consistent. But at a given N H , the absorption-corrected intrinsic 0 . − L X ratio between the two models ighly Obscured and Compton-thick AGNs in CDFs Table 2.
X-ray spectral fitting results for highly obscured AGNsXID field RA DEC z ztype HR Γ N H L X counts model(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12)8 CDF-N 188.841072 62.250256 2.794 zphot 0.17 1.80f 24 .
42 45.30 39 A10 CDF-N 188.846392 62.292835 1.173 zphot 0.07 1.80f 23 .
88 43.15 41 C11 CDF-N 188.847062 62.217590 1.498 zphot -0.23 1.80f 24 .
63 44.90 51 D14 CDF-N 188.853852 62.256847 1.652 zphot 0.59 1.80f 23 .
81 44.22 105 A16 CDF-N 188.869802 62.240976 1.732 zphot 0.25 1.80f 23 .
37 43.97 197 A
Notes.
Column 1: source ID in the Luo et al. (2017) and Xue et al. (2016) catalogs. Column 2: field. Columns 3 and 4: rightascension (RA) and declination (DEC) for the X-ray source position. Column 5: ZFINAL redshift. Column 6: redshift type(“zspec”: spectroscopic; “zphot”: photometric). Column 7: hardness ratio. Column 8: photon index (“f”: fixed value).Column 9: line-of-sight column density in units of 10 cm − . Column 10: logarithm of absorption-corrected rest-frame2–10 keV luminosity in units of erg s − . Column 11: 0.5–7 keV background-subtracted net counts. Column 12: best-fit model(see Section 3.1). The full version of this table is available online. Table 3.
The best-fit models for highly obscured sources with 0.5–7 keV net counts < ≤ counts < ≥ ≥ cm − < N H < cm − , and N H ≥ cm − , respectivelycases N A B C D
C+D N N s A s B s C s D s C s +D s N s (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) <
100 counts 210 80 21 89 20 51.9% 79 14 11 38 16 68.4%100 −
300 counts 146 51 7 69 19 60.3% 67 17 5 29 16 67.2% ≥
300 counts 80 14 5 52 9 76.2% 45 4 2 31 8 86.7% ≥ cm − < N H < cm −
328 120 17 166 24 57.9% 136 27 10 80 19 72.8% N H ≥ cm −
108 25 16 44 23 62.0% 55 8 8 18 21 70.9% N H ≥ cm − and count ≥
150 16 3 1 6 6 75.0% 9 1 1 1 6 77.8%total 436 145 33 210 48 58.9% 191 35 18 98 40 72.3%
Notes.
Column 1: cases of four count bins and two N H ranges corresponding to CN and CT AGNs. Column 2: total sourcenumber in each case. Columns 3 −
6: numbers of sources that are best-fitted by models A − D, respectively. Note that models Cand D have f ref fixed at 1.0 while models A and B have negligible f ref thus negligible reflection component and Fe K lines, andwe show the fraction of the sources that have reflection and iron emissions lines in Column 7. Columns 8–13: same as Columns2–7, but for the spectroscopic-redshift subsample. dramatically increase with N H . It is obvious that model z underestimates the intrinsic luminosity due to neglectof the Compton scattering process, and the discrepancyis already evident at N H ∼ cm − (see also Burlonet al. 2011). Therefore such models with only photo-electric absorption taken into consideration (e.g., zwabs , zphabs , ztbabs ) must be used with caution in the highlyobscured regime.4.2. Photon Index Distribution
The photon index distribution for the 62 highly ob-scured sources with free-Γ during the fitting process isshown in Figure 6. The mean value of the distributionis 1 . ± .
04, in agreement with previous X-ray spectralanalyses in the CDF-S (e.g., Tozzi et al. 2006; Liu et al.2017) and COSMOS (e.g., Lanzuisi et al. 2013). Notethat our Γ distribution peaks at Γ = 1.4 instead of themean value 1.8. This is because the photon index is re- stricted within 1 . − . N H and L X using the Spearman rank correlation test.No correlation between Γ and redshift (Spearman’s ρ =0.00, p -value = 0.99) is detected, indicating that theinner disk and corona structures have little evolutionacross cosmic time. The Spearman tests also suggest nosignificant correlation between Γ and N H ( L X ).It is noteworthy that there are 8 sources with Γ = 1.4or Γ > C is 9.0 if we set Γ free). Sources with very flat photon Li, Xue, Sun et al.
22 23 24 25 26 log N H (This work) l o g N H ( L )
41 42 43 44 45 46 log L X (This work) l o g L X ( L )
22 23 24 25 26 log N H (This work) l o g N H ( B )
22 23 24 25 26 log N H (This work) l o g N H ( B ) Figure 4.
Comparing the spectral fitting results with Liu et al. (2017), Brightman et al. (2014) and Buchner et al. (2015) forcommon sources with redshift difference ∆ z < .
22 23 24 25
MYTorus log N H (cm ) z w l o g N H ( c m ) (a) NH_zw
17 16 15 14 13 obs. log (0.5-7 keV flux) o b s . l o g ( . - k e V f l u x ) (b)
17 16 15 14 13 int. log (0.5-7 keV flux) i n t . l o g ( . - k e V f l u x ) (c)
42 44 46 log L X (erg s ) l o g L X ( e r g s ) (d)
22 23 24 25 log N H (cm ) l o g ( L X , M YT o r u s / L X , z w ) (e) Figure 5.
Comparing the spectral fitting results between using the MYTorus ( x -axes) and phabs × zwabs × zpow ( y -axes inall panels but e ) models. a - d : Comparison results of log N H , observed 0.5–7 keV flux, intrinsic 0.5–7 keV flux and intrinsic2–10 keV luminosity, respectively. e . The log-ratio of the intrinsic 2–10 keV luminosity between MYTorus and model z as afunction of the MYTorus N H . ighly Obscured and Compton-thick AGNs in CDFs photon index N Figure 6.
The distribution of the best-fit photon index forthe 62 sources with free Γ. The green and red vertical dashedlines show the median (Γ = 1.77) and mean (Γ = 1.82) values,respectively.
Table 4.
Information of sources with extreme photon in-dexes field ID Γ N H counts ztype zqualCDFS 98 1.40 0.24 +0 . − . +0 . − .
617 zspec insecure172 1.40 0.33 +0 . − .
534 zphot ...243 1.40 0.12 +0 . − .
372 zspec secure249 1.40 0.18 +0 . − .
830 zphot ...597 1.40 0.24 +0 . − .
351 zspec secure760 1.40 0.72 +0 . − .
351 zspec insecureCDF-N 546 1.40 0.11 +0 . − .
456 zspec secure
Notes.
Columns are the same as in Table 2. The zqual columnrepresents the quality of the spectroscopic redshift (Secure or Inse-cure). For sources with insecure spectroscopic redshifts, the X-raysource catalogs may choose to adopt the photometric redshifts asZFINAL (see Section 4.4 in Xue et al. 2011 for details). indexes are often considered to be reflection-dominatedCT candidates and their extremely obscured nature maynot be revealed by our best-fit N H (Georgantopouloset al. 2009, 2011). Moreover, since we are dealing withthe stacked spectra here, sources which possessed largespectral variability may also exhibit untypical spectralshape, as might be the case for XID 249 (see Section7.4). Additionally, we cannot rule out the possibilitythat the extreme photon indexes of some sources arewrongly measured due to insecure photometric redshifts.4.3. Hardness Ratio versus Redshift
For highly obscured AGNs, the significant absorptionand scattering of soft X-ray photons may lead to largehardness ratio (HR), thus the HR may be used as an in-dicator of obscuration level (e.g., Wang et al. 2004). InFigure 7 we show the observed HR (here we adopt thedefinition of HR as (H − S)/(H+S), where H stands forthe observed-frame 2–7 keV count rate and S stands forthe 0.5–2 keV count rate) as a function of redshift for z H a r d n e ss curve A : = 2.0, log N H = 23, f scat = 5%curve B : = 2.0, log N H = 24, f scat = 1%curve C : = 2.0, log N H = 24, f scat = 5% Figure 7.
Hardness ratio of the highly obscured sample asa function of source redshift. CT candidates selected in Sec-tion 4.5 with best-fit N H > cm − and the 1 σ lower limitof N H > × cm − are shown in red while other highlyobscured AGNs are shown in blue. The size of the symbolindicates their best-fit N H value. Less obscured sources withbest-fit N H < cm − are shown in gray. Curves in dif-ferent colors represent different selection curves presented inSection 4.3. highly obscured sources and CT candidates confirmedby the spectral fitting as well as less obscured AGNs(best-fit N H < cm − ) that contaminate our sam-ple. Heavily obscured sources have significantly largerHR than less obscured sources as expected. We alsocalculate the effective photon index Γ eff (obtained byfitting the spectra using XSPEC model phabs × zpow with N H fixed at the Galactic value) for the HOS. 90%of them have Γ eff < eff is only − .
38 for CT AGNs, which again verify their heavilyobscured nature.The anti-correlation between HR and z is due to thatfor high redshift sources, the observed soft band actuallycorresponds to a much harder band in the rest-frameand the photons may not be absorbed significantly. Notethat the HRs of the most heavily obscured CT AGNs arenot always the largest at a given redshift. This can beascribed to the additional soft excess component softensthe CT AGN spectra, and Compton-thin (CN) AGNswith intrinsically flat photon indexes are also plausibleto have larger HRs than CT AGNs.0 Li, Xue, Sun et al.
22 23 24 25020406080 N Candidate (1%)Confirmed (1%)Candidate (5%)Confirmed (5%) log N H (cm ) f r a c t i o n z log L X erg s log f exs Completeness (1%)Completeness (5%)
Figure 8.
Top: Distributions of source properties of the candidate samples selected from different HR-curves (Candidate) andthose confirmed to be heavily obscured through spectral fitting (Confirmed). Bottom: Completeness fractions of the candidatesamples selected by different HR-curves. The 1% and 5% in the legend represent different soft excess fractions assumed whilederiving the HR-curves in Figure 7. Sources without a detected soft-excess component are shown as log f exs = − To test whether a simple HR value can be used toselect heavily obscured AGNs, we calculate the evolutionof HR as a function of redshift for typical X-ray spectralparameters. Since our purpose is to find the critical HRfor highly obscured AGNs as a function of redshift, wesimply consider the threshold condition: a source with N H = 10 cm − . Additionally, we set f ref at 1.0, thehigh-energy cutoff at 200 keV and fix the photon indexΓ of the two power-laws at 2.0. The constant f exs is setto 5% that is typical for AGNs with N H ∼ cm − (e.g., Brightman et al. 2014). The reason why we adoptΓ = 2 . ∼ (cid:46) .
0. Therefore, with this value we expectthat our derived HR curve will be able to successfullyselect out most of highly obscured AGNs. The modelis shown in Figure 2a and the derived model HRs atdifferent reshifts are shown in green in Figure 7 (curveA).It is worth noting that not all CT AGNs have largerHRs than CN AGNs. We also show the HR curve de-rived for a CT AGN with N H = 10 cm − . This time f exs is chosen to be 1% due to the fact that f exs de-creases with increasing N H (see Section 5.2). The modelis shown in Figure 2b and the simulated HR curve is alsoshown in Figure 7 (curve B). It is clear that at low red-shifts, even a very small soft excess fraction can easilydominate the soft X-ray spectra that could make CTAGNs even softer than CN AGNs. Therefore, to avoid missing a large number of high- N H sources at low red-shifts, we use curve B at z ≤ . z > . N H onlyslightly smaller than 10 cm − and lying close to theboundary curve, as expected. The redshift, luminos-ity and soft excess fraction distributions for the selectedcandidates are very similar to sources which are con-firmed to have N H > cm − .In addition, we define the completeness fraction asthe ratio of the amount of highly obscured sources (con-firmed) lie above the selection curve to the total highlyobscured source amount (i.e., all red and blue pointsin Figure 7). The completeness fractions as a functionof several parameters are shown in the second row inFigure 8. The average completeness fraction for the se-lection curve is 88% (382/436) and remains >
90% for ighly Obscured and Compton-thick AGNs in CDFs f exs , the completeness frac-tion significantly drops to < − . z + 0 . z + 0 . z − . , ( z ≤ . − . z + 0 . z − . z + 0 . , ( z > . . (1)We note that the conclusions in the following sectionsremain largely unchanged if we use the HR-selected sam-ple instead, suggesting that a simple HR value can beused to identify highly obscured AGNs and select a rep-resentative highly obscured sample as long as the red-shift and soft excess effects are properly taken into ac-count while deriving the boundary curve.4.4. Redshift and Luminosity Distributions
In Figure 9 we show the redshift distribution of theHOS in the right-bottom panel. 191 sources have spec-troscopic redshifts and 245 sources have photometricredshifts. The peak appears at z ∼ − .
5, which isknown as the peak epoch of both star-formation andblack hole accretion activities (e.g., Aird et al. 2015).The source amount slightly decreases down to z ∼ . z < .
5. Since
Chandra onlypointed at small patches of sky in these two surveys, thesmall number of sources detected at z < . N H ∼ cm − ) sources shifts out of the Chandra soft bandpass at z ∼ .
5, making it more diffi-cult to constrain N H (Vito et al. 2018).The L X distribution for the HOS is displayed in thetop-left panel of Figure 9. The log L X ( L X in unitsof erg s − ) peaks in the range of 43.5–44.0 with themean value of 43 . ± .
03. The luminous sources with log L X > . ± .
06, because ofthe minimum detectable L X for a flux-limited survey sig-nificantly increases with increasing N H (see Figure 16).The correlations among L X , N H and redshift are alsoplotted in Figure 9. The well-known Malmquist biascan be clearly seen in the bottom-left L X - z corner andwe will discuss the influence of this bias to our resultsin Section 6.4.4.5. Observed N H Distribution
The observed distributions of the best-fit N H in sixluminosity bins of the HOS are shown in Figure 10.The small amount of high N H sources observed in thesmaller L X bins will be clearly illustrated in Section6.2 (see Figure 16) that the minimum detectable lu-minosity significantly increases with N H . We show thenormalized cumulative N H distributions in Figure 11a.Sources with best-fit N H > . × cm − accountfor ∼
25% (108/436) of the HOS. The distribution ob-tained in Brightman et al. (2014), which utilized the4 Ms CDF-S data, is shown in the same plot for com-parison. The K-S tests imply that the N H distributionsfrom the CDF-S and the CDF-N are consistent with p -value = 0.96, but are different from that in Brightmanet al. (2014) with p -value (cid:28) . N H = 23.5 – 24.0 cm − while the distribution in Brightman et al. (2014) peaksat log N H = 23.0 – 23.5 cm − .Although we have significantly improved the photonstatistics for each source using the deepest data, weshould caution that most of CT AGNs discovered stillhave counts <
200 and their N H might be poorly con-strained. Therefore we adopt a more conservative cri-terion that a source will be regarded as CT candidateonly if it has the best-fit N H > cm − and the1 σ lower limit on N H is constrained to be greater than5 × cm − . This criterion selects 102 sources, with66 from CDF-S and 36 from CDF-N, accounting for ∼
23% of the HOS, to be CT candidates. Note thatthere are some well-studied CT candidates in the liter-ature that are not classified as CT AGNs in our sampledue to their best-fit N H < cm − , such as XID 328,XID 551 (XIDs 153, 202 in Comastri et al. 2011, respec-tively), XID 539 (XID 403 in Gilli et al. 2011) and XID375 ( BzK N H > × cm − , which shows good consistency with previousstudies, thus we confirm their heavily obscured natureusing deeper Chandra data.2
Li, Xue, Sun et al. N log L X = 43.65 +0.700.75 CT .
22 3 .
62 4 .
02 4 .
42 4 . l o g N H N log N H = 23.60 +0.620.39 log L X z . . . . . log N H z N z = 1.88 +1.010.92 CTvar
Figure 9.
The triangle plot of the correlations among L X , N H and redshift. The outer parts of the triangle show the distributionsof L X , N H and redshift from the top-left to the bottom-right panels, respectively. The vertical dashed lines show the medianvalue for each distribution. In the redshift and L X distribution panels, the solid histograms represent the total sample and thedotted ones show the distributions of CT candidates selected in 4.5. The shaded histogram shows the redshift distribution of ourvariability sample in Section 7. The inner parts of the triangle show the scatter plots and corresponding density maps overlaidwith contours among L X , N H and redshift, respectively. In the L X versus z plane, the green rectangle shows the subsample weselect in Section 6.5 to investigate the evolution of the intrinsic fraction of CT AGNs among highly obscured ones. We note that at f −
10 keV > − erg cm − s − , thenumber of our CT AGNs in the CDF-N (22 sources) ishigher than Georgantopoulos et al. (2009) which pre-sented 10 CT candidates in the same field. This maybecause of Georgantopoulos et al. (2009) mainly focusedon searching reflection-dominated CT AGNs and possi-bly missed some transmission-dominated sources. Wealso fit the 22 CT candidates using the BNtorus model(Brightman & Nandra 2011) and find that 20 sources have best-fit N H > cm − and f −
10 keV > − erg cm − s − , which confirms the MYTorus result.4.6. Number Counts for Compton-thick AGNs
We calculate the observed number counts (log N − log S ) for CT candidates selected in Section 4.5. Theobserved log N − log S is defined as the observed sourceamount divided by survey area. However, as we willdiscuss in Section 6.2, the sky coverage is not a con-stant value across the whole flux and N H ranges, hence ighly Obscured and Compton-thick AGNs in CDFs
23 24 25 log N H (cm ) N log L X < 43
23 24 25 log N H (cm ) N
43 < log L X < 43.543.5 < log L X < 44
23 24 25 log N H (cm ) N
44 < log L X < 44.544.5 < log L X < 45log L X > 45 Figure 10.
The observed N H distributions in six L X bins. The small number of high N H sources in low L X bins is due to thatsuch sources are difficult to detect. This bias will be corrected in Section 6.4. (a) log N H (cm ) N CDF-SCDF-NBrightman 4Ms CDF-S 10 (b) obs. 2-10 keV flux (erg cm s ) N ( > S , d e g ) Gilli+07Akylas+12Ueda+14Balleytane+11CDF-S correctedCDF-N correctedBrightman+12 (4Ms CDFS)Lanzuisi+18 (COSMOS)Kocevski+18 (X-UDS) (c) obs. 2-10 keV flux (erg cm s ) l o g N H ( c m ) CDF-SCDF-N
Figure 11.
Left: Normalized cumulative N H distribution for our highly obscured sample. The 4Ms CDF-S result adopted fromBrightman et al. (2014) is also plotted for comparison. Middle: The observed log N − log S for CT candidates compared withprevious number counts in the 4 Ms CDF-S (Brightman & Ueda 2012), X-UDS (Kocevski et al. 2018) and COSMOS (Lanzuisiet al. 2018) as well as several CXB model predictions (Gilli et al. 2007; Ballantyne et al. 2011; Akylas et al. 2012; Ueda et al.2014). Right: Column density of CT AGNs as a function of observed 2–10 keV flux in the CDF-S (blue) and CDF-N (red),respectively. we assign a weighting factor 1 /ω area (see Section 6.2 fordetails) for each source while calculating their cumula-tive number counts. The corrected results are shown inFigure 11b, and several number counts measurementspresented in previous works in the 4 Ms CDF-S (Bright-man & Ueda 2012), X-UDS (Kocevski et al. 2018) andCOSMOS (Lanzuisi et al. 2018) fields are also displayedfor comparison. Note that Brightman & Ueda (2012)reported the result in the 0.5–8 keV band and we con-vert the 0.5–8 keV flux into 2–10 keV by assuming aΓ = − . f −
10 keV > × − erg cm − s − , possibly caused by the lownumber statistics and the cosmic variance due to thesmall sky coverage of CDFs (e.g., Harrison et al. 2012).While at the faint end, both CDF-S and CDF-N numbercounts flatten due to the limited detection ability, andthe CDF-N number counts are significantly smaller thanthe CDF-S, because the three times shallower exposurelimits its detectability of faint CT sources (see Figure11c). Our results are also in agreement with the 4 MsCDF-S result at the faint end, the COSMOS result atthe bright end, and the X-UDS result at moderate fluxes.4 Li, Xue, Sun et al.
We also compare the observed log N − log S with sev-eral CXB model predictions in z = 0 − , Akylas et al.(2012) and Ueda et al. (2014) CXB models, the lumi-nosity range used to derive the predicted number countsis log L X = 42 − . − , similar to our sample;while for the Ballantyne et al. (2011) model, the pre-dicted result is presented in log L X = 41 . −
48 erg s − .As can be seen from Figure 11b, our observed log N − log S prefers the moderate CT number counts as pre-dicted by Akylas et al. (2012) and Ueda et al. (2014),while other models more or less overestimate or under-estimate the number counts.In summary, the observed parameter distributions andrelationships presented in this section provide a basicdescription of the highly obscured AGN population andare crucial for distinguishing various CXB models. How-ever, these observed distributions, in particular, the ob-served N H distribution, are influenced by several biases.We will discuss more details about these biases and re-construct the intrinsic N H distribution representative forthe highly obscured AGN population in Section 6. THE ORIGIN OF HEAVY X-RAYOBSCURATION5.1.
X-ray Highly-Obscured Broad-Line AGNs
We collect optical classification results for our CDF-S sample by cross-matching their optical/NIR/IR/radiocounterpart positions presented in the X-ray source cat-alogs (CP RA and CP DEC) with the Silverman et al.(2010) E-CDF-S optical spectroscopic catalog using a0.5 (cid:48)(cid:48) matching radius. Among the 55 matched sources,31 sources have been classified. To our surprise, 19%(6/31) of them are labeled as broad-line AGN (BLAGN).The detailed information of these sources are summa-rized in Table 5.Most of the six sources have sufficient counts and re-liable redshift measurements to constrain N H . Threesources have insecure spectroscopic redshifts althoughthey are labeled as BLAGN, and the 7 Ms main cata-log chooses to adopt the photometric redshift as ZFI-NAL for two of them. This could be due to that theyhave low S/N optical spectra or only one emission lineis available for sources within particular redshift ranges,which makes it difficult to conclusively determine the ∼ gilli/count.html/. We assume ahigh- z declined LF (Vito et al. 2018). http://indra.astro.noa.gr/xrb.html. We assume a 40% CTfraction and the default 4.5% reflection fraction. ∼ yueda/xrb2014.html line nature, thus giving an insecure redshift. For threesources with insecure redshifts, we set redshift as a freeparameter in the spectral fitting and obtain the best-fitX-ray redshift and corresponding N H . All of them arestill best-fitted by N H > cm − . The redshift rangein which the source will remain X-ray highly obscuredis also listed. In general, the N H estimates for the sixsources are robust. This result suggests that the heavyX-ray absorption in a fraction of sources is largely a l.o.seffect caused by some compact clumpy clouds obscuringthe central X-ray emitting region, but the global cov-ering factor (CF) for the high- N H materials is limitedthus the BLR is not blocked; or maybe the heavy X-rayobscuration is produced by the BLR itself.5.2. Soft Excess Fraction Dependences
There are 80 sources in our sample that require a sec-ond power law to fit the soft excess component. Theorigin of this component is still a puzzle. Many differentmodels, e.g., warm Comptonization or blurred ionizedreflection from the disk, partially ionized absorption in awind, or power-law continuum in other directions scat-tered into the line of sight from large-scale Compton-thin matters, are proposed to explain the diverse situa-tions found in different sources (e.g., Boissay et al. 2016).However, for highly obscured AGNs, the scattered mech-anism is preferred since either blurred reflection or warmComptonization from the accretion disk will be signifi-cantly attenuated by the obscuring materials.We show the f exs , which represents the relative nor-malization of the soft component with respect to theintrinsic power law, as a function of N H in Figure 12.We find a significant anti-correlation between the twoparameters: Spearman’s ρ = − . p -value = (cid:28) . f exs is 3.7% for the total sample. The scat-tered fraction decreases from 6.1% to 4.3% and 1.3% for N H < × cm − , 5 × cm − < N H < × cm − and N H > × cm − , respectively.To verify that this anti-correlation is not caused by theparameter degeneracy that sources with a large soft ex-cess fraction and a high N H might be misclassified as low N H , low f exs and low intrinsic luminosity sources owingto the low S/N, we assume that CT AGNs have exactlythe same f exs as CN AGNs (here we simply choose f exs =5%). Then we generate fake spectra using model D fora sample of simulated sources, which have redshifts and N H uniformly distributed between z = 0 − z = 0 . − < log N H <
25 cm − with an interval of ∆log N H = 0 .
25. For each source wesimulate 100 spectra with 0.5–7 keV counts ∼
50, 100and 150, respectively. We then fit the fake spectra andperform Spearman rank correlation test for the output ighly Obscured and Compton-thick AGNs in CDFs Table 5.
Information of X-ray highly obscured BLAGNs in the CDF-S
XID model f ref Oclass z ztype zqual counts log N H (cm − ) z -range z x log N H , x (cm − )(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12)882 A (cid:28) (cid:28) z > . z > . Notes.
Column 4: optical classification result from Silverman et al. (2010). 31 X-ray highly obscured AGNs in our CDF-S sample haveoptical classification results and 6 of them are identified as BLAGNs. Column 10: the redshift range for sources with insecure redshifts tobe determined as being X-ray highly obscured. Column 11: the best-fit X-ray redshift by setting redshift as a free parameter duringspectral fitting. Column 12: best-fit log N H when adopting the X-ray best-fit redshift. Other columns have the same meaning as Table 2. N H and f exs . The Spearman’s ρ ≈ . f exs . Moreover, Ueda et al. (2015) and Kawa-muro et al. (2016) found that the [O III] and [O IV] tohard X-ray luminosity ratios are lower in low scatter-ing fraction sources, inferring that the low f exs AGNsare buried in small opening angle torus. Therefore,the anti-correlation between f exs and N H indicates thathigh N H sources might preferentially reside in high CFtoruses. This is also consistent with the results fromBrightman & Ueda (2012) and Lanzuisi et al. (2015),who found a similar anti-correlation between f exs and N H and explained it as highly obscured AGNs being alsoheavily geometrically obscured (but see the discussion inthe next section).5.3. Reprocessed Components and Covering Factor
Unlike local CT AGNs in which the prominent neutralFe lines and reflection hump (hereafter the reprocessedcomponents) are prevalently detected in the high-qualityhard X-ray spectra and can be served as unambiguoussignatures of being Compton-thick (e.g., Tanimoto et al.2018; Zhao et al. 2019; Marchesi et al. 2019), there are41% of our highly obscured sources and 38% of the CTcandidates having negligible f ref (see Table 3), respec-tively. It should be noted that the low detection rateof reprocessed components is not in conflict with theirhighly obscured nature. The poor S/N for low-countsources and the narrow Chandra spectral coverage makeit challenging to detect narrow iron lines and to distin-guish the transmission and reflection components, whichmay cause a significant overestimation of the fractions. N H (10 cm ) f e x s Figure 12.
The soft excess fraction f exs as a function of N H . The blue points show the individual sources and thered points show the binned results. The negative correla-tion between f exs and N H indicates that a portion of high N H sources might have higher “torus” covering factors thatmakes the soft excess photons hard to escape, under the as-sumption that the excess originates from the scattered-backcontinuum in highly obscured AGNs. But we note that it is possible that highly obscuredAGNs may have weak reprocessed components if theglobal CFs of high- N H materials are low or the cen-tral AGN is geometrically fully buried by CT materialssuch that even the reflected photons cannot escape (e.g.,Brightman et al. 2014). To distinguish different scenar-ios, we use the Borus model (Balokovi´c et al. 2018) ,which allows CF to vary freely, to fit the spectra forall sources with counts > tor , which is defined as the cosine valueof the opening angle θ tor measured from the symmetry In the Borus model, CF tor = 1 . tor = 0 . ∼ mislavb/download/). Li, Xue, Sun et al. CF tor N (a) f ref f ref = 1 photon index (MYTorus) p h o t o n i n d e x ( B o r u s ) log N H (MYTorus) l o g N H ( B o r u s )
42 43 44 45 log L X (MYTorus) l o g L X ( B o r u s ) k e v ( p h o t o n s c m s k e V ) ID : 746 (b) total, CF=0.1transmissionreflectiontotal, CF=1.0transmissionreflection12
CF=0.1
Energy (keV) CF=1.0 r e s i d u a l Figure 13. (a) The covering factor derived through the Borus model (Balokovi´c et al. 2018) and the comparisons betweenphoton index, N H , and L X obtained from the Borus and MYtorus models. The photon indexes in this plot are the originalresults without fixing Γ at 1.8 for some extreme cases. (b) The spectral fitting result for CDF-S XID 746 using the Borus model.Both the fully buried model (CF tor = 1 .
0; Cstat = 942.3) and the weak torus model (CF tor = 0 .
1; Cstat = 941.8) can wellreproduce the data. The contribution from the reflection component to the total spectrum in both models are largely negligible. axis towards the equatorial plane, is shown in Figure13a. We also show the comparisons of the main spectralfitting parameters with those obtained from the MY-Torus model and find good consistency.For sources with negligible reprocessed components inthe MYTorus model ( f ref (cid:28) tor ≈ . f ref = 1.0 are mostly best fitted bya highly-covered model (0 . < CF tor < .
0) that havetorus covering angles (90 ◦ − θ tor ) between 65 ◦ − ◦ .We emphasize that though the l.o.s and the global torus N H are linked in the spectral fitting, those fully-buriedsources do not need to be covered by very high N H ma-terials in all directions, since the best-fit N H will morelikely converge to the values determined by the l.o.s com-ponent because the photoelectric absorption is a muchstronger spectral feature. Moreover, we find that theaverage CF tor for f exs < .
05 sources is 0.60; while for f exs > .
05 sources, the CF tor is 0.32. This might pro-vide evidence for that the lower soft-excess fraction insome sources is caused by a geometrically more buriedstructure (see Section 5.2).However, we caution here that the current S/N andspectral coverage do not allow us to constrain CF fromX-ray spectral analysis. If we fix CF tor at 0.1 for sourcesthat are best fitted by a fully-buried model, we may stillobtain acceptable fits with Cstat values only slightly in-creasing (an example is shown in Figure 13b). Conclu-sively disentangling several scenarios is beyond the scopeof this work, and multi-wavelength data are needed tofurther shed light on the nature of the obscuring mate-rials (Li et al., in preparation). INTRINSIC N H DISTRIBUTION FOR HIGHLYOBSCURED AGNSIn Section 4.5, we present the observed N H distri-bution which is affected by several sample incomplete-ness. In order to derive the intrinsic N H distribution, weshould take into account the errors on best-fit N H andcorrect for several survey biases (Liu et al. 2017). Toperform such corrections, we restrict our analysis to the ighly Obscured and Compton-thick AGNs in CDFs < (cid:48) , for the sakeof avoiding large background contamination, extremelysmall sky coverage (see Section 6.2) and limited de-tectable fraction (see Section 6.4).6.1. Errors on Best-fit N H To consider the uncertainty of the spectral fitting,we perform a resampling procedure to the best-fit N H .Given the asymmetric errors, we assume that the errorson N H obey the “half-gaussian” distribution. We firstgenerate 1000 N H for each source with the σ of the half-gaussian distribution equals to the lower 1 σ error, andthen generate another 1000 N H with σ equals to the up-per 1 σ error. The mean value µ is set to the best-fit N H .The resampled N H distribution (calculated by averag-ing the 2000 resampled distributions and also correctedfor the sky coverage effect; see Section 6.2) is shown inthe shaded region of Figure 18. It has an extended taildown to N H < cm − which contains only 4.0% ofthe resampled data, suggesting that even if consideringthe spectral fitting errors, most of our sources are stillconsistent with being highly obscured.6.2. Sky Coverage Effect
Due to
Chandra ’s instrument features, the point-source PSF size increases and the effective exposuretime dramatically decreases toward large off-axis angles.The sensitivities of detecting faint sources reduce promi-nently at the outskirt of the field, leading to small skycoverage while the observed flux is low. Therefore, acorrection must be made.We calculate the energy flux to count rate conversionfactor (ECF) by assuming the soft-excess model (modelD) with Γ, f ref and f exs fixed at 1.8, 1.0 and 1%, respec-tively. Combining ECF and the exposure map, we builda sensitivity map that represents the flux limits corre-sponding to the 20-count cut. Using this map, we calcu-late the sky coverage as a function of observed 0.5–7 keVflux, N H and redshift as shown in Figure 14. Then wemeasure the sky coverage for each source based on theirobserved flux, N H and redshift obtained from spectralfitting. We define ω area as the ratio between the sourcesky coverage and the maximum sky coverage of the two Chandra surveys (484.2 arcmin and 447.5 arcmin inthe CDF-S and the CDF-N, respectively). To correctfor the sky coverage effect, we simply weigh each sourceby 1 /ω area while resampling the observed N H distribu-tion. To avoid extremely large weights, we apply an-other cut for the 22 sources with sky coverage less than100 arcmin , setting their weighting factors to the me-dian factor 1.3, since these sources lie close to the countcut and we cannot rule out the possibility that their extremely small sky coverage may be a result from in-appropriate spectral modeling.6.3. Eddington Bias
Considering the measurement error of net counts,sources with intrinsically low counts may exceed ourcount cut (20 photons), and sources with high countsmay be missed in our sample. Since the number of faintand bright sources are not equal, this error leads tothe Eddington bias. To correct for this bias, we con-volve the cumulative count distribution for AGNs inthe two
Chandra survey catalogs with the count errorsusing Poisson sampling in order to obtain the intrinsiccount distribution. We follow the ‘pseudo-deconvolved’method proposed in L17 (see Section 5.1.3 of L17 for de-tails), by shifting the observed curve leftward with thevalue equal to the displacement between the observedand convolved curves. We then obtain the deconvolvedcurve that represents the intrinsic cumulative count dis-tribution.The difference between the deconvolved and the ob-served curves can be treated as the number of sourcesmissed or mis-included (depending on the adopted countcut) in our sample. As shown in Figure 15, at a 20-photon count cut, we miss 11.6 and 6.5 sources in theCDF-S and CDF-N, respectively. This is due to the factthat though sources with counts <
20 may outnumberthe bright sources, such faint sources are not detected inthe source catalogs. By controlling the sample size (i.e.,the cumulative source number in the y -axes) to be thesame, we obtain the effective net counts of 22.4 and 21.1for the CDF-S and CDF-N, respectively. Thus the 20-photon count cut will be replaced by the new effectivecount cut while performing other corrections.6.4. Malmquist Bias
So far we have only obtained the corrected observed N H distribution for our sample. In order to derivethe intrinsic N H distribution for the highly obscuredAGN population, we should also consider the unde-tectable part of the Chandra survey, which is know asthe Malmquist bias that sources with lower luminositieswill be missed by a flux-limited survey. More specifi-cally, given the count rate limit of the survey (the ef-fective count cut in our situation), whether a source isdetectable or not depends on its redshift, spectral shape,column density, intrinsic luminosity as well as its posi-tion in the field, since large off-axis angles lead to sig-nificant lower sensitivities.To quantify this bias, we generate fake X-ray spec-tra to find out the minimum L X ( L cut , corresponding toour effective count cut) in different redshift and N H bins,8 Li, Xue, Sun et al.
Observed 0.5-7 keV flux (10 erg cm s ) s k y c o v e r a g e ( a c r m i n ) CDF-S z = 1.5, log N H = 23.0 z = 1.5, log N H = 23.6 z = 1.5, log N H = 24.2 z = 2.5, log N H = 23.0 z = 2.5, log N H = 23.6 z = 2.5, log N H = 24.2 Observed 0.5-7 keV flux (10 erg cm s ) s k y c o v e r a g e ( a c r m i n ) CDF-N z = 1.5, log N H = 23.0 z = 1.5, log N H = 23.6 z = 1.5, log N H = 24.2 z = 2.5, log N H = 23.0 z = 2.5, log N H = 23.6 z = 2.5, log N H = 24.2 Figure 14.
Sky coverage as a function of the observed 0.5–7 keV flux corresponding to the count cut 20 in the CDF-S andCDF-N, shown for three different N H values at redshifts, respectively. s o u r c e nu m b e r CDF-S observedconvolveddeconvolved
20 30 40 50 60 70 80 90 100
Full band net counts r e s i d u a l
50 150 250 counts e rr o r f r a c t i o n s o u r c e nu m b e r CDF-N observedconvolveddeconvolved
20 30 40 50 60 70 80 90 100
Full band net counts r e s i d u a l
50 150 250 counts e rr o r f r a c t i o n Figure 15.
Illustration of the Eddington bias in the CDF-S and CDF-N, respectively. In each panel/field, the red curve showsthe smoothed cumulative count distribution for all AGNs with counts >
20 and off-axis angle < (cid:48) . The blue curve is obtainedby Poisson-sampling the observed count distribution by considering the measurement errors. For sources without count errormeasurements in the catalogs, we fit the mean error fraction f err ( f err = count error / count) in a given count range (shown asorange points) as a function of full-band net counts using sources with errors; the fitted curve is shown in the inset. Thenwe simply assign an error to sources without errors, which is given by their total net counts times the corresponding f err . Weapply a pseudo-deconvolved method to obtain the deconvolved intrinsic count distribution by shifting the blue curve leftward orrightward (based on the location on the right or left side of the node of the two curves) with the value equal to the displacementbetween the blue and red curves, as shown as the green curve. From the residual between the observed and deconvolved curves,we can see that we missed 11.6 and 6.5 sources (shown as red crosses in the bottom panels), which correspond to the effectivecount cuts of 22.4 and 22.1 (shown as red crosses in the top panels) in the CDF-S and CDF-N, respectively. ighly Obscured and Compton-thick AGNs in CDFs z l o g L X ( e r g s ) CDF-S log N H = 23.0, @3log N H = 24.0, @3log N H = 24.5, @3log N H = 23.0, @6log N H = 24.0, @6log N H = 24.5, @6log N H = 23.0, @8log N H = 24.0, @8log N H = 24.5, @8 z l o g L X ( e r g s ) CDF-N log N H = 23.0, @3log N H = 24.0, @3log N H = 24.5, @3log N H = 23.0, @6log N H = 24.0, @6log N H = 24.5, @6log N H = 23.0, @8log N H = 24.0, @8log N H = 24.5, @8 Figure 16.
Detection boundary curves as a function of redshift for different N H and off-axis angles, corresponding to theeffective count cut 22.4 in the CDF-S and 21.1 in the CDF-N, respectively. l o g L X ( e r g s )
41 42 43 44 45 46 z d / d l o g L X M p c d e x MiyajiGeorgakakis l o g L X ( e r g s )
41 42 43 44 45 46 z d / d l o g L X M p c d e x Aird l o g L X ( e r g s )
41 42 43 44 45 46 z d / d l o g L X M p c d e x Ueda
Figure 17.
Left: The combined X-ray LFs taken from Miyaji et al. (2015) at z < z > − < log N H <
24 cm − ) AGN LFs from Aird et al. (2015). Right: The Compton-thin AGN(log N H <
24 cm − ) LFs from Ueda et al. (2014). Li, Xue, Sun et al. which can be detected at a specific position in the field.We divide the two CDFs into 10 annuli, with each cor-responding to an interval of 1 (cid:48) . While running XSPECsimulations, we use the averaged exposure time, rmf fileand arf file calculated from all sources in each annulus.Model D is used in the simulation with the parametersset as the same as in Section 6.2.We show some of the detectable boundary curves inFigure 16. We define the detectable region as the areabetween the boundary curve and the maximum log L X =46 erg s − . As we can clearly see, sources with high N H and large off-axis angles have significantly small de-tectable regions. We follow Equation 3 in T06, F ( N H ) dN H = f ( N H ) dN H (cid:90) z max dVdz dz × (cid:90) L max L cut( N H ,z ) N ( L X , z ) dL X , (2)where f ( N H ) represents the intrinsic N H distribution, F ( N H ) represents the observed N H distribution, and N ( L X , z ) is the X-ray LF. The detectable source frac-tion for each N H bin relative to the log N H = 23 cm − bin in a given redshift interval can be calculated as: f = (cid:82) L max L cut( N H ,z ) N ( L X , z ) dL X (cid:82) L max L cut(23 , z) N ( L X , z ) dL X . (3)Thus for the bin we choose, f = (cid:40) , if 23 . − ≤ log N H < .
25 cm − ,< , if log N H ≥ .
25 cm − . (4)Note that T06 derived the above Equation 2 by assum-ing that f ( N H , L X , z ), which represents the possibilityof detecting a source with column density in the range of N H to N H + d N H for a given L X and z , varying slowlyas L X and z change. There is evidence that the obscuredfraction of AGNs evolves with both luminosity and red-shift. However, we do not consider such a complicatedevolution and simply extract f ( N H , L X , z ) from the in-tegral. Therefore, by multiplying 1 /f to the observed N H distribution F ( N H ) dN H in each bin, sources withdifferent N H are corrected to cover the same observ-able space with respect to the log N H = 23 − .
25 cm − bin and the Malmquist bias is corrected. To avoid thecase of the corrected distribution being dominated bysources with extremely small detectable fractions, weapply a cut that while f < /f will be simply set to 1/0.25. Adopting a lower (0.2)or higher (0.3) detectable fraction cut will lead to ∼ z < z> − < N H <
24 cm − ). We do not use theirLFs for CT AGNs, since they obtained a systematicallylower CT fraction than other works and our result inSection 4.6 prefers higher CT number counts (but seethe discussion in Section 6.4 of Aird et al. (2015) thatthis low CT fraction should not be ruled out) and thespectral model they used may not be appropriate forCT AGNs (see Section 4.1). The third LF we used isfrom Ueda et al. (2014). The three LF models are shownin Figure 17. We note that the adopted obscured AGNLFs may not be suitable for CT AGNs, but an extensiveexploration of this issue is beyond the scope of this work.6.5. Intrinsic N H Distribution
Figure 18 shows the intrinsic N H distributions in fiveredshift bins after correcting for all the aforementionedbiases. We define the intrinsic CT-to-highly-obscuredsource fraction f CH as f CH = N (log N H ≥ N (log N H ≥ , (5)and list its evolution with redshift in Table 6.As we can see, different LFs give consistent results.It is obvious that the corrected N H distributions arequite different from the observed ones, especially in theCT regime. This highlights the importance of carefullycorrecting for the observational biases in order to un-cover the underlying distributions. Though the intrin-sic N H distributions are different between high-redshiftand low-redshift bins, the CT-to-highly-obscured frac-tion is in accordance with no evident redshift evolutiongiven the uncertainties (i.e., f CH ≈ .
52 for all redshiftbins; see Table 6). Note that Buchner et al. (2015)also claimed that the CT fraction is constant acrosscosmic time, but the CT fraction in their work is de-fined as the number of CT AGNs to the total AGNpopulation. To avoid possible biases induced by thefact that different redshift bins are sensitive to differentluminosity ranges, we also select two subsamples with43 . − < log L X < . − at z = 1 − z = 2 −
3, respectively, which have similar L X distribu-tions (see Figure 9), and perform the same corrections.The results are listed in Table 6 and are still consistentwith no redshift evolution.Vito et al. (2018) presented the N H distribution forthe high-redshift ( z >
3) AGNs also in the CDFs. Their ighly Obscured and Compton-thick AGNs in CDFs log N H (cm ) N z =0.0-0.8 observed allobserved sub ( <10 )corrected: Miyaji & Georgakakiscorrected: Airdcorrected: Uedaresampled 22.5 23.0 23.5 24.0 24.5 25.0 25.5 log N H (cm ) N z =0.8-1.6 log N H (cm ) N z =1.6-2.4 log N H (cm ) N z =2.4-3.0 log N H (cm ) N z =3.0-5.0
23 24 25 log N H cm VitoLi log N H (cm ) N z =0.0-5.0 Figure 18. N H distributions for our HOS in five redshift bins and the entire redshift range ( z = 0 − N H distribution;the black histogram shows the observed N H distribution for sources with off-axis angle < (cid:48) ; the blue shaded area shows theresampled N H distribution that takes into account the spectral fitting errors and sky coverage effect; and the red and blue solidhistograms and purple dashed histogram show the intrinsic N H distributions after correcting for the undetectable parts of thetwo surveys using the Ueda et al. (2014), Aird et al. (2015), as well as Miyaji et al. (2015) & Georgakakis et al. (2015) LFs,respectively. The comparison of the normalized N H distribution of highly obscured AGNs at z > z = 3 ∼ z = 0 − N H distribution shown in blue is resampled into the shaded histogram by considering N H errors. analyses were based on modeling the X-ray spectra usingthe wabs × zwabs × zpow model for the z > z = 3 . − . N H bin for which we have to make assumptions regard-ing the unknown AGN population. The series of cor-rections can be simply understood as multiplying anevolved weighting factor to the observed distribution, which have several limitations. First of all, as mentionedbefore, different redshift bins sample different luminos-ity ranges. Several authors have found that the obscuredfraction changes with luminosity (Buchner et al. 2015;Lusso et al. 2013), thus the N H distribution might beluminosity-dependent. However, in Equation 2, we ig-nore this effect to simplify the calculation. Second, forthe extremely buried sources ( N H > cm − ), or in-trinsically very faint sources, the corrections cannot bemade since there are hardly any such sources observeddue to the limit of the survey and the restrictions im-posed by our sample-selection count cut. Therefore, thepart of the N H distribution of sources that are under thedetection limit is actually still missed in our final result.Furthermore, our identification method of CT AGNs isbased on the absorption turnover imprinted in the spec-2 Li, Xue, Sun et al.
Table 6.
The evolution of the observed ( f obs ) and the cor-rected intrinsic ( f int ) CT-to-highly-obscured fraction withredshift. The values listed below are calculated by averagingthe results from the three LF models. z f obs , CH f int , CH log L X (erg s − )(1) (2) (3) (4)Total sample0.0 − . ± .
07 42.90.8 − . ± .
05 43.51.6 − . ± .
05 43.72.4 − . ± .
06 43.93.0 − . ± .
09 44.10.0 − . ± .
06 43.6Subsample1.0 − . ± .
05 43.72.0 − . ± .
05 43.8
Notes.
Top table: the results for the total sample. Bottomtable: the results for two subsamples with43 . − < log L X < . − at z = 1 . − . z = 2 . − .
0, which have similar L X distributions (see Figure 9). tra. For a source with N H > cm − , the turnoveroccurs at rest-frame ∼
20 keV. This implies that for low-redshift sources, the characteristic feature of CT AGNsmay not be observable in the
Chandra spectra of limitedenergy range, which is coupled with the uncertaintiesinduced by the photometric redshifts, leading to greatchallenges in correctly measuring the obscuration level.It is plausible that N H > cm − sources may con-tribute significantly to heavily obscured AGN popula-tion (Risaliti et al. 1999), but such sources are missedin our sample and the corrections at this part is beyondour attainment. Therefore, the N H distribution at thehighest N H bin displayed in Figure 18 is highly uncer-tain and incomplete so that the CT-to-highly-obscuredfractions in Table 6 should be better treated as lowerlimits. SPECTRAL VARIABILITY ANALYSESX-ray variability is a ubiquitous feature among AGNsthat can provide useful information about AGN prop-erties and relevant underlying physical processes (e.g.,Vaughan et al. 2003; McHardy et al. 2006; Gonz´alez-Mart´ın & Vaughan 2012; Soldi et al. 2014; Zheng et al.2017). In this section, we aim to study the main drivingmechanism that causes the variability of highly obscuredAGNs by investigating the variability of some main spec-tral parameters, such as the luminosity, column densityand reflection strength.To perform detailed long-term spectral variabilityanalyses of highly obscured AGNs and expand the sam-ple size, we select our sources from the HOS based on only one additional criterion, i.e., the net 0.5–7 keVcounts in the total 7 Ms exposure in the CDF-S or inthe total 2 Ms exposure in the CDF-N should be largerthan 700 or 900, respectively. We bin data from neigh-boring observations as one epoch to enhance the S/N.For sources in the CDF-S that have counts > < counts < ∼ combine spectra tool in CIAOto generate the source spectra, background spectra, rmfand arf files in each epoch.7.1. Method
We use the same spectral fitting models as describedin Section 3.1 to perform spectral variability analyses,except for removing the constants f ref and f exs in themodel and untying the normalizations of the intrinsicpower law and other spectral components. We simulta-neously fit the spectra in each epoch using the best-fitmodel for each source obtained in Section 4, but thistime we allow norm ref to vary freely to search for possi-ble variation in the Compton-scattered continuum. Con-sidering the relatively low counts in each epoch, the fit-ting strategies are adopted as follows:1. The uncertainties on Γ can cause large degeneracies ifwe set it free in all epochs. Given that some sourcesmay have large variability in Γ, we first let Γ in eachepoch vary freely and obtain C free . Then for the firsttwo adjacent epochs (i.e., epoch1 and epoch2), we linkΓ and measure the ∆ C with respect to C free . If ∆ C> . o . f = 1), Γ is considered to be differentbetween the two epochs at >
95% confidence level,and we will set Γ free; otherwise, we will link Γ. Thenfor the second epoch pair (i.e., epoch2 and epoch3),we link their Γ and compare the Cstat value with thelast step. After traversing all the epoch pairs, if novariability is detected, Γ is linked together.2. The reflection and soft excess components are oftenconsidered to be produced in large-scale clouds, e.g.,torus, for highly obscured AGNs. Since the times-pan in the rest-frame is (1+ z ) times less than in the ighly Obscured and Compton-thick AGNs in CDFs Table 7.
Observational epochs, variability sample and bin-ning information
CDF-S Start time End time CDF-N Start time End timeI 1999 Oct 2000 Dec I 1999 Nov 2000 NovII 2007 Sep 2007 Nov II 2001 Feb 2001 NovIII 2010 Mar 2010 Jul III 2002 Feb 2002 FebIV 2014 Jun 2015 JanXID z N M counts N obs I II III IV(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)CDF-S 4 epochs49 2.394 4 A 2345 98 0.93 0.96 1.96 2.7081 3.309 4 A 3579 102 0.93 0.96 1.96 2.8898 1.412 4 B 3689 102 0.93 0.96 1.96 2.88186 2.810 4 A 2532 102 0.93 0.96 1.96 2.88214 3.740 4 A 2863 102 0.93 0.96 1.96 2.88236 2.562 4 A 1578 102 0.93 0.96 1.96 2.88308 1.097 4 B 2484 102 0.93 0.96 1.96 2.88458 2.291 4 C 2093 102 0.93 0.96 1.96 2.88746 3.064 4 A 2144 102 0.93 0.96 1.96 2.88760 3.350 4 B 1974 102 0.93 0.96 1.96 2.88876 3.470 4 A 4191 102 0.93 0.96 1.96 2.88933 1.654 4 A 1792 102 0.93 0.96 1.96 2.88CDF-S 3 epochs63 0.737 3 B 848 29 0.45 0.56 1.0491 2.256 3 A 756 29 0.45 0.56 1.04752 0.733 3 A 866 55 0.80 1.27 1.69785 1.600 3 C 1345 22 0.31 0.36 0.6073 2.509 3 A 879 102 1.89 1.96 2.88240 1.185 3 C 785 102 1.89 1.96 2.88249 0.735 3 A 830 102 1.89 1.96 2.88328 1.536 3 C 1432 102 1.89 1.96 2.88367 0.604 3 C 948 102 1.89 1.96 2.88399 1.730 3 C 1077 102 1.89 1.96 2.88551 3.700 3 C 756 102 1.89 1.96 2.88621 1.213 3 A 784 102 1.89 1.96 2.88658 1.845 3 A 839 102 1.89 1.96 2.88733 2.404 3 D 973 102 1.89 1.96 2.88818 2.593 3 C 711 102 1.89 1.96 2.88840 1.220 3 C 1321 102 1.89 1.96 2.88846 2.483 3 A 974 102 1.89 1.96 2.88CDF-N 3 epochs66 0.959 3 A 934 20 0.49 0.87 0.58143 1.727 3 A 1516 20 0.49 0.87 0.58
Notes.
Top: Definition of the observational epochs. Bottom: Vari-ability sample and the binning information. Column (3): number ofbinning epochs. Column (4): best-fit model. Column (5): total 0.5–7 keV net counts. Column (6): number of individual observations ofthe source during the total CDF-S 7 Ms or CDF-N 2 Ms exposure.Columns (7 – 10): exposure time in units of Ms in each epoch. observed frame and the typical redshift of our sam-ple is relatively high, it is reasonable to assume thatin such a short timescale, the large-scale componentsmay not vary dramatically. To better constrain thenormalization of the intrinsic power law, we tie thenormalizations of the reflection component and thesoft excess component in each epoch, respectively, un-less setting them free leading to ∆
C > χ .4. The N H and normalization of the intrinsic power lawalways vary freely.The simultaneous fitting yields the best-fit model pa-rameters Γ, norm ref and norm scat . Then we fit the spec-tra in each epoch with Γ i , norm ref , i and norm scat , i fixedat the best-fit value obtained in the simultaneous fit-ting to obtain the N H , observed flux, intrinsic flux andthe corresponding errors in each epoch. The spectralfitting parameters of each epoch are listed in Table 8.No significant photon index variability is detected forall sources.7.2. L X , N H and Observed Flux Variability We identify L X , N H and observed flux variable sourcesbased on the classical χ test. For illustration, we ex-plain how we determine L X -variable sources. First, weuse the cf lux model in XSPEC to calculate the ob-served and absorption-corrected 0 . − . − L . − =4 πd L f . − , int (1+ z ) Γ int − , where the “int” subscriptrepresents the intrinsic value. Since the photon indexand redshift are all fixed during spectral fitting of eachsingle epoch, the error on L . − is totally attributedto the error of the intrinsic flux f . − , int . Thus the χ of the rest-frame L . − actually equals the χ ofthe observed-frame f . − , int . If a source does notvary, i.e., the intrinsic luminosity (or N H , flux) is con-stant, assuming that the errors obey the Gaussian dis-tribution, the χ value calculated from χ = N (cid:88) i =1 ( x i − x ) σ x,i (6)should obey the χ distribution with ( N −
1) d.o.f, where x represents the tested parameter ( L X , N H or flux), σ x,i means the 1 σ error in the i th epoch, and N is the num-ber of epochs. The χ test results are listed in Table 9.4 Li, Xue, Sun et al.
Table 8.
Multi-epoch spectral fitting results for the variability sample
XID Γ N a H1 N H2 N H3 L bX1 L X2 L X3 flux c flux flux CDF-S63 1.54 0.06 +0 . − . +0 . − . +0 . − . +1 . − . +1 . − . +0 . − .
73 1.81 0.64 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
91 2.08 0.16 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
240 1.40 1.80 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
249 1.51 0.50 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
328 2.12 1.91 − . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
367 2.07 0.42 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
399 1.40 0.24 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
551 1.84 1.56 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
621 2.11 0.18 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
658 1.67 0.15 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
733 1.84 0.21 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
752 2.24 0.15 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
785 1.97 0.54 +0 . − . +0 . − . +0 . − . +2 . − . +2 . − . +1 . − .
818 1.97 0.14 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
840 1.60 0.58 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
846 1.56 0.14 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . CDF-N66 1.92 0.34 +0 . − . +0 . − . +0 . − . +1 . − . +0 . − . +0 . − .
143 1.46 0.15 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . XID Γ N H1 N H2 N H3 N H4 L X1 L X2 L X3 L X4 flux flux flux flux CDF-S49 1.66 0.19 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
81 1.83 0.17 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
98 0.98 0.18 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +1 . − . +0 . − . +0 . − .
186 1.79 0.22 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
214 1.87 0.26 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
236 1.62 0.15 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
308 1.69 0.21 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
458 1.83 0.19 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
746 1.68 0.41 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
876 1.87 0.15 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
933 1.71 0.31 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − .
760 1.49 0.81 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . Notes. a: in units of 10 cm − . b: 2 −
10 keV intrinsic luminosity, in units of 10 erg s − . c: observed 0.5–7 keV flux, in units of10 − erg cm − s − . When the source with 3 (4) binning epochs satisfies χ > >
98% con-fidence level. We also check the L X - and N H -variableidentification results using the Akaike information cri-terion (AIC) method as described in Y16 that does notneed to assume the Gaussian errors (see Section 3.2.3of Y16 for details). The ∆AIC values for each source The choice of this confidence level is because Y16 showedthat it roughly corresponds to ∆AIC = 4, so that we can directlycompare the χ results with those obtained using the AIC method. are also listed in Table 9. Sources with ∆AIC > . L X and N H variability while using the χ test andthe AIC method, respectively. We decide to use the χ test results in the following analysissince for the two discrepant sources: CDFS XID 933 and CDFNXID 66, the AIC results show neither L X nor N H variability, whichis inconsistent with their flux-variable nature. ighly Obscured and Compton-thick AGNs in CDFs log L (erg s ) -0.2-0.10.00.10.20.30.4 n xv intrinsic n xv observed Figure 19.
Normalized excess variance ( σ ) of the ob-served (top) and intrinsic (bottom) 0.5–7 keV flux versusthe intrinsic 0.5–7 keV luminosity for the variability sam-ple. The original values without fixing negative σ at 0 aredisplayed. The error bars are calculated using Equation 8. Based on the χ test results, sources that show ob-served flux variability account for 55 ±
13% (17/31) ofthe entire sample. The resulting L X and N H variablesource fractions are 29 ±
10% (9/31) and 19 ±
8% (6/31),respectively. We do not find any correlation between L X and N H variability patterns, suggesting that the mainreason that causes the variation of N H is likely not thechange of ionization parameter induced by the variable L X , but is likely the occultation of the clumpy cloudsmoving in/out of our l.o.s.7.3. Variability Amplitude Estimation
We use the normalized excess variance σ (nxv) toestimate the variability amplitude. σ is defined as σ = 1 N (cid:104) x (cid:105) N (cid:88) i =1 [( x i − (cid:104) x (cid:105) ) − ( δx i ) ] , (7)where N is the number of binning epochs (3 or 4), x i and δx i are the best-fit parameters (observed flux, L X , N H )and their 1 σ errors, (cid:104) x (cid:105) is the unweighted mean of thecalculated parameter. The error on σ is estimated as s D / ( (cid:104) x (cid:105) √ N ), where s = 1 N − N (cid:88) i =1 [( x i − (cid:104) x (cid:105) ) − ( δx i ) − σ (cid:104) x (cid:105) ] . (8) Table 9.
The Chi-squared value χ , the ∆AIC value ∆, andthe normalized excess variance σ of the observed 0.5–7 keVflux, intrinsic rest-frame L . − and N H . XID χ flux χ L X ∆ L X χ N H ∆ N H σ flux σ L X σ N H CDF-S 4 Epochs49 53.3 21.9 11.9 17.3 11.3 0.074 0.056 0.07381 69.1 23.7 13.9 3.1 − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − Notes.
Negative values of σ are set to 0. Sources with χ > . χ > . > . >
98% confidence level (Y16).
We calculate σ and corresponding errors on the in-trinsic rest-frame 0.5–7 keV luminosity ( σ ,L ), ob-served 0.5–7 keV flux ( σ ,F ) and N H ( σ ,N ) for eachsource. The results are listed in Table 9. To checkwhether our calculation of σ is affected by the lim-ited counts (e.g., Zheng et al. 2017; hereafter Z17), weperform Spearman rank correlation tests and find nosignificant correlation between σ ,F , σ ,N , σ ,L and counts with Spearman’s ρ = 0 . p -value = 0.70; ρ = − . p -value = 0.56 and ρ = 0 . p -value = 0.08,respectively.We also calculate the fractional root-mean-square(frms) variability amplitude, which is defined as (cid:112) σ for sources having σ > Li, Xue, Sun et al. that our sample lacks sources that display large vari-ability amplitudes. The mean frms values of the totalsample for the observed flux, intrinsic luminosity and N H are 12%, 10% and 6%, respectively. While forsources having positive nxv, the mean (maximum) frmsvalues are 17% (29%), 19% (34%), and 16% (39%),respectively.In addition, we do not find any correlation between σ ,F , σ ,L and L X (Figure 19) whereas other stud-ies suggested a negative trend, i.e., luminous AGNsshow relatively small variations (Y16; Z17; Paolillo et al.2017). We perform a Spearman correlation test forthe 25 common sources using the σ ,F derived inZ17 which is obtained from the light curve analysis,the result is still consistent with no correlation between L X and σ ,F (Spearman’s ρ = − . p -value = 0.13).As pointed out by Allevato et al. (2013), the calcula-tion of σ is biased at low counts and uneven cadenceand should be restricted to large samples. Therefore,the limited source number and the broad redshift spanmay prevent us from detecting such a correlation. Butbeyond that, since Y16 and Z17 mainly focused on theAGNs with largest counts available in the CDF-S ratherthan highly obscured AGNs, the variability behavior andits underlying physical drivers of our sources and thosein Y16 and Z17 may be intrinsically different.Recently, Gonz´alez-Mart´ın (2018) found a new X-rayvariability plane for AGNs that links the characteris-tic break frequency of the PSD ( f break ) with bolomet-ric luminosity and N H . This makes N H play a non-negligible role in understanding AGN variability. As wewill show in the next section, although Z17 argued thatthe N H variation does not contribute significantly to thetotal variability in their sample, the N H variability doesplay a substantial role for highly obscured AGNs. Nev-ertheless, the small σ values indicate that the state ofhighly obscured AGNs does not vary dramatically whileconcerning the average source properties on a 17-yeartimescale in the observed-frame.7.4. Detailed Variability Analysis
There are some sources in our sample that show signif-icant variability of the observed flux, intrinsic luminosityor N H . In this section, we perform detailed analyses tostudy their variability behavior aiming at shedding lighton the leading mechanism that drives the large variabil-ities, and try to understand the typical location and sizeof the obscuring clouds (all the following sources are inthe CDF-S) . XID 63
This source has 848 counts available and its data arebinned to three epochs. It is a moderately luminous source ( ¯ L X = 5 . × erg s − ) at z = 0 . . +0 . − . ) plus an additional soft excess compo-nent. Γ and norm soft do not show variations and arelinked during all the epochs. The unfolded spectra, lightcurves, as well as the 1 σ and 2 σ confidence contours of N H and normalization are displayed in the left columnof Figure 20. It has χ N H = 65 . χ f, obs = 34 . χ L = 7 .
2, indicating significant variabilities. Althoughthe photon index does not vary, the observed spectralshape changes prominently due to the large variationin column density. The absorption is weak in epoch1with N H = 5 . × cm − . Then it transforms intoa highly obscured state with N H = 1 . × cm − inepoch2, and its N H continues to rise to 1 . × cm − in epoch3. The intrinsic flux increases about a factorof 1.4 from epoch1 to epoch2, and decreases about afactor of 1.6 from epoch2 to epoch3. Due to the largevariability amplitude in N H , the observed flux remainsconstant in the first two epochs, but significantly de-clines in epoch3 by a factor of 1.9. The N H variabilitymay be caused by the clumpy cloud moving into thel.o.s. Since the X-ray obscuration type transition hap-pens between epoch1 and epoch2, which corresponds toa rest-frame timescale t < . ≈ r and14 - 195 keV luminosity from Koshida et al. (2014) toestimate the distance. The relation can be describedas log r = − .
04 + 0 . L −
195 keV / ), where r and L − are in units of pc and erg s − , respec-tively. The L −
195 keV = 2 . × erg s − at epoch1 isobtained by extrapolating the power-law to the higherenergy band. Thus the estimated r from the inner torusto the central black hole is ≈ .
14 pc. By assum-ing Keplerian motion, the orbital period of the cloud t orbit = 2 π r ( M BH G ) ≈
519 ( M BH M (cid:12) ) − year. The angu-lar size (viewed from the BH) of the cloud is estimatedas t tran t orbit × ◦ < . ◦ ( M BH M (cid:12) ) . By multiplying r , weestimate the cloud size to be < .
01 ( M BH M (cid:12) ) pc. XID 249
This source can be well fitted by the absorbed power-law model with Γ = 1 . +0 . − . and has the most promi-nent N H variability in our sample ( χ N H = 130.8) asshown in Figure 20. The N H declines from 5 . × cm − to 2 . × cm − , and finally 1 . × cm − during the three epochs. The intrinsic flux of thissource remains roughly constant ( χ flux, in = 2 . ighly Obscured and Compton-thick AGNs in CDFs k e v ( p h o t o n s c m s k e V ) ID : 6312
Epoch1 Epoch2
Energy (keV)12
Epoch3 r e s i d u a l k e v ( p h o t o n s c m s k e V ) ID : 24912
Epoch1 Epoch2
Energy (keV)12
Epoch3 r e s i d u a l k e v ( p h o t o n s c m s k e V ) ID : 4912
Epoch1 Epoch2 Epoch3
Energy (keV)12
Epoch4 r e s i d u a l o b s f l u x ID : 63 p w l f l u x Epoch1 Epoch2 Epoch3 N H o b s f l u x ID : 249 p w l f l u x Epoch1 Epoch2 Epoch3 N H o b s f l u x ID : 49 p w l f l u x Epoch1 Epoch2 Epoch3 Epoch41234 N H N H (10 cm ) n o r m ( p o w e r - l a w ) ID : 63 N H (10 cm ) n o r m ( p o w e r - l a w ) ID : 249 N H (10 cm ) n o r m ( p o w e r - l a w )
12 34
ID : 49
Figure 20. first row : The spectra for three highly variable sources XID 63, XID 249 and XID 49 in the CDF-S. The top panelsshow the unfolded spectra in each epoch and the bottom panels show the data-to-model ratios. second row : The observed0.5–7 keV flux, intrinsic 0.5–7 keV flux (in units of 10 − erg cm − s − ) and N H (in units of 10 cm − ) variability curvesfor the three sources, respectively. third row : Confidence contours of the normalizations of the intrinsic power law and N H forthe three sources. The solid and dashed curves indicate 1 σ and 2 σ confidence contours, respectively. The numbers annotatedrepresent different epochs. Li, Xue, Sun et al. the strong variation in N H causes a significant increasein the observed flux ( χ flux, obs = 31 . N H varia-tion. Using the same method as for XID 63, we esti-mate the distance and the cloud size to be < .
09 pcand < .
006 ( M BH M (cid:12) ) pc for the occultation event fromepoch2 to epoch3. XID 328
This source has the highest column density in ourvariability sample. The epoch-mean L X and N H are3 . × erg s − and 8 . × cm − , respectively,which are consistent with the stacked spectral analy-ses. The spectra can be well described by a power law(Γ = 2 . +0 . − . ) plus Compton reflection continuum andstrong Fe K lines. N H and the reflection flux remainroughly constant in the timespan of 5.8 years in therest-frame, and the observed flux variations result fromthe intrinsic X-ray power variability. According to thestacked spectral analyses, the reflection flux accounts for27.5% of the observed 0.5–7 keV flux, but only accountsfor 3.8% of the intrinsic 0.5–7 keV flux. This source hasalso been reported in the literature as CT candidates(Comastri et al. 2011; Y16). Our results confirm itshighly obscured nature but with higher N H , L X and Γ. XID 49
This source was observed 98 times during the total7 Ms campaign and has 2345 counts available. Thespectra can be explained by a simple absorbed powerlaw (Γ = 1 . +0 . − . ) with different column densitiesand normalizations (see Figure 20). The intrinsic fluxdecreases from 1 . × − erg cm − s − to 7 . × − erg cm − s − during the first two epochs and fi-nally rises again to 1 . × − erg cm − s − in thefourth epoch. N H has a wave-shape variability behaviorand does not follow the intrinsic flux variability pattern,indicating that the varying absorption is not caused bythe change of ionization state. XID 458
This source has a redshift of 2.291 and can be de-scribed by a power law (Γ= 1 . +0 . − . ) with reflectionhump and iron emission lines. Though the N H is not sohigh during all the epochs ( ¯ N H ≈ . × cm − ), thereflection component is required to fit the 7 Ms stackedspectrum. The intrinsic flux increases from epoch1 toepoch3, and declines at epoch4. But the reflection fluxremains roughly constant from epoch1 to epoch2, de-creases at epoch3, and does not show any variability from epoch3 to epoch4. The different variability pat-terns indicate that there may be a time lag betweenthe intrinsic continuum variability and the large-scalereflection flux variability. The decline of the reflectionflux during epoch2 to epoch3 (from observed-frame 2007September to 2010 July, corresponding to rest-frame0.76 years) might result from significant flux declinebefore epoch1, and this decline has just propagated tothe reflection medium (possibly the clumpy torus) andcauses significant diminishment of the reflection flux.This means that the variable continuum signal needsat least 2.4 years (rest-frame timespan from 1999 Octo-ber to 2007 September) to spread to the torus from thecentral emission region, which provides a rough lower-limit estimate, ≈ REMAINING SOURCES
By applying similar analyses to the remaining sources,we find that 17 sources in our sample that show observedflux variabilities can be classified into three types: 29%(5/17) are caused by the change of N H , 53% (9/17) arecaused by the AGN intrinsic luminosity variability, and6% (1/17) are driven by both effects. Note that there aretwo sources (CDF-S XIDs 186 and 876) identified to beflux-variable but showing neither L X - nor N H -variabilityaccording to the χ tests. However, their best-fit in-trinsic luminosities both show obvious variations, butdue to large errors, the corresponding χ L values aresmaller than the critical value. The high N H -variablefraction among flux-variable sources confirms our previ-ous thoughts in Section 7.3 that the N H variability is akey ingredient for investigating the variability in highlyobscured AGNs. CONCLUSIONSIn this study, we present systematic X-ray spectralanalyses of 436 highly obscured AGNs ( N H > cm − )with L X > erg s − identified in the Chandra
DeepFields, which make up the largest dedicated highly ob-scured AGN sample to date, to explore their physicalproperties and evolution. We also carry out detailedlong-term variability analyses for a subsample of 31sources with largest counts available to investigate themain driver of their spectral variability and the typicalvariability amplitude. Below we summarize our mainresults.1. We perform detailed X-ray spectral fitting for 1152AGNs in the
Chandra
Deep Fields with observed-frame 0.5–7 keV net counts >
20 using physicallyappropriate MYTorus-based models, in order to iden-tify heavily obscured ones. By limiting our analyses ighly Obscured and Compton-thick AGNs in CDFs L X >
42 erg s − in order to avoidpossible contamination from star-forming galaxies,436 AGNs are confirmed to be highly obscured, with z = 1 . L X = 43 . − .2. We select 102 Compton-thick (CT) candidates withbest-fit N H > cm − and 1 σ lower limit > × cm − , accounting for ∼
23% of the highly ob-scured sample. The observed log N -log S for CTAGNs prefers the moderate CT number counts aspredicted by Akylas et al. (2012) and Ueda et al.(2014) cosmic X-ray background models, while othermodels (e.g., Gilli et al. 2007; Ballantyne et al. 2011)more or less overestimate or underestimate the num-ber counts.3. We present a new hardness-ratio measure of the ob-scuration level as a function of redshift (the HRcurve), which can be used to select heavily obscuredAGN candidates without resorting to detailed spec-tral fitting. The completeness and accuracy by ap-plying the HR curve on the CDFs AGN populationto identify highly obscured ones are 88% and 80%,respectively.4. We find a strong negative correlation between the softexcess fraction f exs and N H with a Spearman rankcorrelation coefficient ρ = − .
66. By assuming thatthe soft excess originates from the scattered-back con-tinuum and treating the small f ref as an indicator ofhigh covering factor of the obscuring materials, thisresult indicates that a portion of the most heavilyabsorbed AGNs reside in an extremely geometricallyburied environment.5. Among the 31 CDF-S highly obscured sources whichhave optical classification results, 19% (6/31) of themare labeled as BLAGN, indicating that at least forpart of the AGN population, the high-level X-ray ob-scuration is largely a line-of-sight effect, i.e., somehigh-column-density clouds on various scales (butnot necessarily a dust-enshrouded torus) may obscurethe compact X-ray emitter without blocking the en-tire BLR; alternatively, the heavy X-ray obscurationmaybe produced by the BLR itself.6. After considering the errors on the best-fit N H andcorrecting for the sky coverage effect, the Eddingtonbias and the N H -dependent Malmquist bias, we de-rive the intrinsic N H distribution representative of thehighly obscured AGN population as well as its evolu-tion across cosmic time. The intrinsic CT-to-highly-obscured fraction is roughly 52% and is consistentwith no evident redshift evolution. 7. We select 31 sources with 0.5–7 keV net counts > >
900 in the CDF-S and CDF-N, respectively,to perform long-term ( ≈
17 years in the observedframe) spectral variability analyses. We find thatthe flux-variable, L X -variable and N H -variable sourcefractions are 55 ±
13% (17/31), 29 ±
10% (9/31) and19 ±
8% (6/31), respectively. The typical flux, L X and N H variability amplitudes for those variable sourcesare 17%, 19%, and 16%, respectively.8. We calculate the normalized excess variance ( σ )of N H , observed 0.5–7 keV flux and L X to investigatethe intrinsic variability amplitude. No correlation be-tween σ and L X is detected, possibly due to theobserved flux variability for a significant fraction ofhighly obscured AGNs being caused by the change of N H rather than the variation of L X alone, as well asthe limited sample size and the broad redshift span.9. We discuss detailed variability behaviors for 5 sourcesthat show significant N H , L X or observed flux vari-ability. The typical locations and sizes of the occul-tation and reflection clouds are estimated (see Sec-tion 7.4). The main driver for the variability of the17 flux-variable sources can be classified into threetypes: 29% (5/17) are caused by the change of N H ,53% (9/17) are caused by the AGN luminosity vari-ability and 6% (1/17) are driven by both effects. Twosources are not classified due to large measurementerrors.Benefiting from the deepest X-ray surveys to date, ourwork provides meaningful constraints on the propertiesand evolution of the AGN obscuring materials over abroad redshift range, and quantifies the detailed vari-ability behaviors of these hidden sources, which are cru-cial for us to better understand the role that highlyobscured AGNs played in galaxy evolution. However,X-ray data alone are incapable to show the overall per-spective. In a subsequent paper of this series (Li et al.,in preparation), we will further explore the properties ofhighly obscured AGNs and their host galaxies by com-bining the wealth of multi-wavelength data in the CDFs,aiming at testing the merger-triggered AGN-galaxy co-evolution scenario (Sanders et al. 1988).We thank the anonymous referee for valuable com-ments and suggestions. We thank David Ballantyne,Murray Brightman, Giorgio Lanzuisi and YoshihiroUeda for providing relevant data, and Tahir Yaqoob andMislav Balokovi´c for helpful discussions. J.Y.L., Y.Q.X.,M.Y.S., and X.C.Z. acknowledge support from the 973Program (2015CB857004), the National Natural Sci-ence Foundation of China (NSFC-11890693, 11473026,0 Li, Xue, Sun et al.
Aird, J., Coil, A. L., Georgakakis, A., Nandra, K., & Barro,G. 2015, MNRAS, 1892, 451Akylas, A., Georgakakis, A., Georgantopoulos, I.,Brightman, M., & Nandra, K. 2012, A&A, 546, A98Akylas, A., Georgantopoulos, I., Ranalli, P., et al. 2016,A&A, 594, A73Alexander, D., & Hickox, R. 2012, NewAR, 56, 93Alexander, D. M., Bauer, F. E., Brandt, W. N., et al. 2003,AJ, 126, 539Alexander, D. M., Chary, R., Pope, A., et al. 2008, ApJ,687, 835Alexander, D. M., Bauer, F. E., Brandt, W. N., et al. 2011,ApJ, 738, 44Allevato, V., Paolillo, M., Papadakis, I., & Pinto, C. 2013,ApJ, 771, 9Ballantyne, D. R., Draper, A. R., Madsen, K. K., Rigby,J. R., & Treister, E. 2011, ApJ, 736, 56Balokovi´c, M., Brightman, M., Harrison, F. A., et al. 2018,ApJ, 854, 42Baumgartner, W. H., Tueller, J., Markwardt, C. B., et al.2013, ApJS, 207, 19Bianchi, S., Guainazzi, M., & Chiaberge, M. 2006, A&A,448, 499Boissay, R., Ricci, C., & Paltani, S. 2016, A&A, 588, A70Brandt, W. N., & Alexander, D. M. 2015, A&A Rv, 23, 1Brightman, M., & Nandra, K. 2011, MNRAS, 413, 1206Brightman, M., Nandra, K., Salvato, M., et al. 2014,MNRAS, 443, 1999Brightman, M., & Ueda, Y. 2012, MNRAS, 423, 702Broos, P. S., Townsley, L. K., Feigelson, E. D., et al. 2010,ApJ, 714, 1582Buchner, J., Georgakakis, A., Nandra, K., et al. 2015, ApJ,802, 89Burlon, D., Ajello, M., Greiner, J., et al. 2011, ApJ, 728, 58Cash, W. 1979, ApJ, 228, 939 Comastri, A. 2004, in Astrophysics and Space ScienceLibrary, Vol. 308, Supermassive Black Holes in theDistant Universe, ed. A. J. Barger, 245Comastri, A., Ranalli, P., Iwasawa, K., et al. 2011, A&A,526, L9Corral, A., Della Ceca, R., Caccianiga, A., et al. 2011,A&A, 530, A42Corral, A., Georgantopoulos, I., Watson, M. G., et al. 2014,A&A, 569, A71Daddi, E. 2007, ApJ, 670, 173Del Moro, A., Alexander, D. M., Bauer, F. E., et al. 2016,MNRAS, 456, 2105Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature,433, 604Fanali, R., Caccianiga, A., Severgnini, P., et al. 2013,Monthly Notices of the Royal Astronomical Society, 433,648Feruglio, C., Daddi, E., Fiore, F., et al. 2011, ApJL, 729, L4Fiore, F., Brusa, M., Cocchia, F., et al. 2003, A&A, 409, 79Fiore, F., Puccetti, S., Brusa, M., et al. 2009, ApJ, 693, 447Fiore, F., Puccetti, S., Grazian, A., et al. 2012, A&A, 537,A16Georgakakis, A., Aird, J., Buchner, J., et al. 2015, MNRAS,453, 1946Georgantopoulos, I., Akylas, A., Georgakakis, A., &Rowan-Robinson, M. 2009, A&A, 507, 747Georgantopoulos, I., Rovilos, E., Xilouris, E. M., Comastri,A., & Akylas, A. 2011, 526, A86Georgantopoulos, I., Comastri, A., Vignali, C., et al. 2013,A&A, 555, A43Gilli, R., Comastri, A., & Hasinger, G. 2007, A&A, 463, 79Gilli, R., Su, J., Norman, C., et al. 2011, ApJL, 730, L28Gonz´alez-Mart´ın, O. 2018, ApJ, 858, 2Gonz´alez-Mart´ın, O., & Vaughan, S. 2012, A&A, 544, A80 ighly Obscured and Compton-thick AGNs in CDFs Goulding, A. D., Alexander, D. M., Mullaney, J. R., et al.2011, MNRAS, 411, 1231Guainazzi, M., Matt, G., & Perola, G. C. 2005, A&A, 444,119Harrison, C. M., Alexander, D. M., Mullaney, J. R., et al.2012, ApJL, 760, L15Hasinger, G. 2008, A&A, 490, 905Hickox, R. C., & Alexander, D. M. 2018, ARA&A, 56, 625Hopkins, P., & Hernquist, L. 2006, ApJS, 166, 1Iwasawa, K., Gilli, R., Vignali, C., et al. 2012, A&A, 546,A84Kawamuro, T., Ueda, Y., Tazaki, F., Ricci, C., &Terashima, Y. 2016, ApJS, 225, 14Kocevski, D. D., Brightman, M., Nandra, K., et al. 2015,ApJ, 814, 104Kocevski, D. D., Hasinger, G., Brightman, M., et al. 2018,ApJS, 236, 48Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511Koshida, S., Minezaki, T., Yoshii, Y., et al. 2014, ApJ, 788,159Koss, M. J., Assef, R., Balokovi´c, M., et al. 2016, ApJ, 825,85Lansbury, G. B., Alexander, D. M., Aird, J., et al. 2017,ApJ, 846, 20Lanzuisi, G., Civano, F., Elvis, M., et al. 2013, MNRAS,431, 978Lanzuisi, G., Ranalli, P., Georgantopoulos, I., et al. 2015,A&A, 573, 137Lanzuisi, G., Civano, F., Marchesi, S., et al. 2018, MNRAS,480, 2578Liu, T., Tozzi, P., Wang, J. X., et al. 2017, ApJS, 451, 457Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2011, ApJ, 740,37Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2017, ApJS, 228,2Lusso, E., Hennawi, J. F., Comastri, A., et al. 2013, ApJ,777, 86Marchesi, S., Lanzuisi, G., Civano, F., et al. 2016, ApJ,830, 100Marchesi, S., Ajello, M., Zhao, X., et al. 2019, ApJ, 872, 8McHardy, I. M., Koerding, E., Knigge, C., Uttley, P., &Fender, R. P. 2006, Nature, 444, 730Miyaji, T., Hasinger, G., Salvato, M., et al. 2015,Astrophys. J., 804, 104 Morganti, R. 2017, Nature Astronomy, 1, 596Murphy, K. D., & Yaqoob, T. 2009, MNRAS, 397, 1549Paolillo, M., Papadakis, I., Brandt, W. N., et al. 2017,MNRAS, 471, 4398Ranalli, P., Comastri, A., Vignali, C., et al. 2013, A&A,555, A42Ricci, C., Ueda, Y., Koss, M. J., et al. 2015, ApJL, 815, L13Ricci, C., Bauer, F. E., Treister, E., et al. 2017a, MNRAS,468, 1273Risaliti, G., Maiolino, R., & Salvati, M. 1999, ApJ, 522, 157Sanders, D. B., Soifer, B. T., Elias, J. H., et al. 1988, ApJ,325, 74Silverman, J. D., Mainieri, V., Salvato, M., et al. 2010,ApJS, 191, 124Soldi, S., Beckmann, V., Baumgartner, W. H., et al. 2014,A&A, 563, A57Stark, A. A., Gammie, C. F., & Wilson, R. 1992, ApJS, 77,79Tanimoto, A., Ueda, Y., Kawamuro, T., et al. 2018, ApJ,853, 146Tozzi, P., Gilli, R., Mainieri, V., et al. 2006, A&A, 451, 457Treister, E., Urry, C. M., & Virani, S. 2009, ApJ, 696, 110Ueda, Y., Akiyama, M., Hasinger, G., Miyaji, T., &Watson, M. G. 2014, 786, 104Ueda, Y., Hashimoto, Y., Ichikawa, K., et al. 2015, ApJ,815, 1Vaughan, S., Edelson, R., Warwick, R. S., & Uttley, P.2003, MNRAS, 345, 1271Vito, F., Brandt, W. N., Yang, G., et al. 2018, MNRAS,473, 2378Wang, J.-M., Watarai, K.-Y., & Mineshige, S. 2004, ApJL,607, L107Wang, J. X., Malhotra, S., Rhoads, J. E., & Norman, C. A.2004, ApJ, 612, L109Xue, Y. Q. 2017, NewAR, 79, 59Xue, Y. Q., Luo, B., Brandt, W. N., et al. 2016, ApJS, 224,15Xue, Y. Q., Luo, B., Bauer, F. E., et al. 2011, ApJS, 195, 10Yang, G., Brandt, W. N., Luo, B., et al. 2016, ApJ, 831, 145Zhang, F., Yu, Q., & Lu, Y. 2017, ApJ, 845, 88Zhao, X., Marchesi, S., Ajello, M., et al. 2019, ApJ, 870, 60Zheng, X. C., Xue, Y. Q., Brandt, W. N., et al. 2017, ApJ,849, 127 Li, Xue, Sun et al.
41 42 43 44 45 46 log L X, new paras. (erg s ) l o g L X ( e r g s )
23 24 25 log N H, new paras. (cm ) l o g N H ( c m ) Figure 21.
Comparison between the spectral fitting results obtained in Section 3 and those obtained in Appendix A using thenew model parameters.
APPENDIX A. JUSTIFICATION OF THE ADOPTED MODELS AND PARAMETERSIn the MYTorus configuration, since the reproccessed components strongly depend on the inclination angle θ , theequatorial N H and the assumed f ref , the best-fit results will also depend on the assumed geometries and parameters.However, due to the low S/N for most of the spectra, we have to make prior assumptions about the input parametersand simply fix them to values that their representativeness has not yet been physically validated. Here we justify thatour spectral fitting results are not significantly influenced by the assumed parameters.We first fit the spectra for the high count sources (counts > θ and f ref set as free parameters. The best-fitted θ peaks at θ < ◦ , therefore this time we choose to fix it at 65 ◦ for all sources. Then we re-fit the spectra bykeeping f ref free and obtain the best-fit results. For the low count sources, we fix f refref