Piercing Through Highly Obscured and Compton-thick AGNs in the Chandra Deep Fields: I. X-ray Spectral and Long-term Variability Analyses

TL;DR

This study analyzes X-ray spectra of over a thousand AGNs in the Chandra Deep Fields to identify highly obscured and Compton-thick AGNs, introduces a new obscuration measure, and examines their variability and evolution.

Contribution

It presents a new hardness-ratio method for selecting highly obscured AGNs with high accuracy and applies it to a large sample, providing insights into their properties and evolution.

Findings

38% of AGNs are highly obscured, including 102 Compton-thick candidates.

The new hardness-ratio measure achieves 88% completeness and 80% accuracy.

The CT-to-highly-obscured fraction is about 52% with no clear redshift evolution.

Abstract

We 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_{\rm H} > 10^{23}\ \rm cm^{-2}$). By fitting spectra with physical models, 436 (38%) sources with $L_{\rm X} > 10^{42}\ \rm erg\ s^{-1}$ are confirmed to be highly obscured, including 102 Compton-thick (CT) candidates. We propose a new hardness-ratio measure of the obscuration level which can be used to select highly obscured AGN candidates. The completeness and accuracy of applying this method to our AGNs are 88% and 80%, respectively. The observed logN-logS relation favors cosmic X-ray background models that predict moderate (i.e., between optimistic and pessimistic) CT number counts. 19% (6/31) of our highly obscured AGNs that have optical classifications are labeled as broad-line AGNs, suggesting that, at least for part of the AGN population, the heavy X-ray obscuration is largely a line-of-sight effect, i.e., some high-column-density clouds on various scales (but not necessarily a dust-enshrouded torus) along our sightline may obscure the compact X-ray emitter. After correcting for several observational biases, we obtain the intrinsic NH distribution and its evolution. The CT-to-highly-obscured fraction is roughly 52% and is consistent with 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, 17 sources show flux variabilities: 31% (5/17) are caused by the change of NH, 53% (9/17) are caused by the intrinsic luminosity variability, 6% (1/17) are driven by both effects, and 2 are not classified due to large spectral fitting errors.

Figures

Tables

Table 1. Table 1: Model definitions used in this work

model	$f_{\mathrm{ref}}$	$f_{\mathrm{exs}}$
A	${10}^{-8}$	0
B	${10}^{-8}$
C	$1.0$	0
D	$1.0$

Table 2. Table 2: X-ray spectral fitting results for highly obscured AGNs

XID	field	RA	DEC	$z$	ztype	HR	$\Gamma$	$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	A
10	CDF-N	188.846392	62.292835	1.173	zphot	0.07	1.80f	$23.88$	43.15	41	C
11	CDF-N	188.847062	62.217590	1.498	zphot	-0.23	1.80f	$24.63$	44.90	51	D
14	CDF-N	188.853852	62.256847	1.652	zphot	0.59	1.80f	$23.81$	44.22	105	A
16	CDF-N	188.869802	62.240976	1.732	zphot	0.25	1.80f	$23.37$	43.97	197	A

Table 3. Table 3: The best-fit models for highly obscured sources with 0.5–7 keV net counts < \textless  100, 100 ≤ counts < 300 100 counts 300 \rm 100\leq counts<300 , counts ≥ \geq 300, counts ≥ \geq 1000, 10 23 ​ cm − 2 superscript 10 23 superscript cm 2 \rm 10^{23}\ cm^{-2}   < \textless   N H < 10 24 ​ cm − 2 subscript 𝑁 H superscript 10 24 superscript cm 2 \textit{$N_{\rm H}$}<\textit{$\rm 10^{24}\ cm^{-2}$} , and N H subscript 𝑁 H N_{\rm H}   ≥ \geq 10 24 ​ cm − 2 superscript 10 24 superscript cm 2 \rm 10^{24}\ cm^{-2} , respectively

cases	$N$	A	B	C	D	$\frac{C+D}{N}$	$N_s$	As	Bs	Cs	Ds	$\frac{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%
$\ge 300$ counts	80	14	5	52	9	76.2%	45	4	2	31	8	86.7%
$\ge 1000$ counts	18	3	2	9	4	72.2%	14	3	1	6	4	71.4%
${10}^{23}{\mathrm{cm}}^{-2}$ $<$ $N_H$ $<$ ${10}^{24}{\mathrm{cm}}^{-2}$	328	120	17	166	24	57.9%	136	27	10	80	19	72.8%
$N_H$ $\ge$ ${10}^{24}{\mathrm{cm}}^{-2}$	108	25	16	44	23	62.0%	55	8	8	18	21	70.9%
$N_H$ $\ge$ ${10}^{24}{\mathrm{cm}}^{-2}$ and  count $\ge$ 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%

Table 4. Table 4: Information of sources with extreme photon indexes

field	ID	$\Gamma$	$N_H$	counts	ztype	zqual
CDFS	98	1.40	0.24${}_{-0.01}{}^{}{}_{}{}^{+0.02}$	3689	zphot	insecure
	135	2.53	1.59${}_{-0.24}{}^{}{}_{}{}^{+0.21}$	617	zspec	insecure
	172	1.40	0.33${}_{-0.08}{}^{}{}_{}{}^{+0.07}$	534	zphot	…
	243	1.40	0.12${}_{-0.01}{}^{}{}_{}{}^{+0.01}$	372	zspec	secure
	249	1.40	0.18${}_{-0.02}{}^{}{}_{}{}^{+0.02}$	830	zphot	…
	597	1.40	0.24${}_{-0.04}{}^{}{}_{}{}^{+0.03}$	351	zspec	secure
	760	1.40	0.72${}_{-0.07}{}^{}{}_{}{}^{+0.08}$	351	zspec	insecure
CDF-N	546	1.40	0.11${}_{-0.01}{}^{}{}_{}{}^{+0.01}$	456	zspec	secure

Table 5. Table 5: Information of X-ray highly obscured BLAGNs in the CDF-S

XID	model	$f_{\mathrm{ref}}$	Oclass	$z$	ztype	zqual	counts	$\log N_H$(${\mathrm{cm}}^{-2}$)	$z$-range	$z_x$	$\log N_{H,x}$ (${\mathrm{cm}}^{-2}$)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
882	A	$\ll 1$	BLAGN	3.19	zspec	Secure	80	23.1	–	–	–
944	C	1	BLAGN	0.96	zphot	Insecure	69	23.9	all	2.00	24.1
16	A	$\ll 1$	BLAGN	3.10	zphot	Insecure	409	23.7	$z>1.4$	1.66	23.2
399	C	1	BLAGN	1.73	zspec	Secure	1077	23.2	–	–	–
968	C	1	BLAGN	2.03	zspec	Secure	500	23.8	–	–	–
977	C	1	BLAGN	4.64	zspec	Insecure	225	23.9	$z>1.2$	3.25	23.9

Table 6. Table 6: The evolution of the observed ( f obs subscript 𝑓 obs f_{\rm obs} ) and the corrected intrinsic ( f int subscript 𝑓 int f_{\rm int} ) CT-to-highly-obscured fraction with redshift. The values listed below are calculated by averaging the results from the three LF models.

$z$	$f_{\mathrm{obs},\mathrm{CH}}$	$f_{\mathrm{int},\mathrm{CH}}$	$\overline{\log L_X}$ ($\mathrm{erg}{s}^{-1}$)
(1)	(2)	(3)	(4)
Total sample
0.0 $-$ 0.8	0.32	$0.57\pm 0.07$	42.9
0.8 $-$ 1.6	0.24	$0.50\pm 0.05$	43.5
1.6 $-$ 2.4	0.24	$0.50\pm 0.05$	43.7
2.4 $-$ 3.0	0.29	$0.54\pm 0.06$	43.9
3.0 $-$ 5.0	0.18	$0.50\pm 0.09$	44.1
0.0 $-$ 5.0	0.25	$0.52\pm 0.06$	43.6
Subsample
1.0 $-$ 2.0	0.25	$0.50\pm 0.05$	43.7
2.0 $-$ 3.0	0.26	$0.51\pm 0.05$	43.8

Table 7. Table 7: Observational epochs, variability sample and binning information

CDF-S	Start time	End time	CDF-N	Start time	End time
I	1999 Oct	2000 Dec	I	1999 Nov	2000 Nov
II	2007 Sep	2007 Nov	II	2001 Feb	2001 Nov
III	2010 Mar	2010 Jul	III	2002 Feb	2002 Feb
IV	2014 Jun	2015 Jan

Table 8. Table 8: Multi-epoch spectral fitting results for the variability sample

XID	$\Gamma$	$N_{\mathrm{H1}}^a$	$N_{\mathrm{H2}}$	$N_{\mathrm{H3}}$	$L_{\mathrm{X1}}^b$	$L_{\mathrm{X2}}$	$L_{\mathrm{X3}}$	flux${}_{}{}^c{}_1{}^{}$	flux2	flux3
CDF-S
63	1.54	0.06${}_{}{}^{+0.01}{}_{-0.01}{}^{}$	0.14${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.19${}_{}{}^{+0.04}{}_{-0.03}{}^{}$	0.47	0.66	0.41	16.0${}_{}{}^{+1.9}{}_{-1.8}{}^{}$	16.6${}_{}{}^{+1.5}{}_{-1.5}{}^{}$	8.8${}_{}{}^{+0.9}{}_{-0.9}{}^{}$
73	1.81	0.64${}_{}{}^{+0.08}{}_{-0.07}{}^{}$	0.56${}_{}{}^{+0.11}{}_{-0.11}{}^{}$	0.67${}_{}{}^{+0.07}{}_{-0.07}{}^{}$	3.37	1.87	3.78	2.7${}_{}{}^{+0.3}{}_{-0.3}{}^{}$	1.7${}_{}{}^{+0.3}{}_{-0.3}{}^{}$	2.9${}_{}{}^{+0.3}{}_{-0.3}{}^{}$
91	2.08	0.16${}_{}{}^{+0.03}{}_{-0.03}{}^{}$	0.12${}_{}{}^{+0.03}{}_{-0.03}{}^{}$	0.14${}_{}{}^{+0.03}{}_{-0.03}{}^{}$	1.29	1.03	1.25	2.3${}_{}{}^{+0.3}{}_{-0.3}{}^{}$	2.1${}_{}{}^{+0.3}{}_{-0.3}{}^{}$	2.4${}_{}{}^{+0.3}{}_{-0.3}{}^{}$
240	1.40	1.80${}_{}{}^{+0.17}{}_{-0.15}{}^{}$	1.41${}_{}{}^{+0.17}{}_{-0.15}{}^{}$	1.34${}_{}{}^{+0.12}{}_{-0.11}{}^{}$	1.76	1.21	1.04	3.7${}_{}{}^{+0.5}{}_{-0.4}{}^{}$	3.4${}_{}{}^{+0.5}{}_{-0.5}{}^{}$	3.3${}_{}{}^{+0.3}{}_{-0.3}{}^{}$
249	1.51	0.50${}_{}{}^{+0.14}{}_{-0.12}{}^{}$	0.28${}_{}{}^{+0.08}{}_{-0.06}{}^{}$	0.14${}_{}{}^{+0.02}{}_{-0.01}{}^{}$	0.25	0.16	0.18	2.1${}_{}{}^{+0.3}{}_{-0.3}{}^{}$	2.4${}_{}{}^{+0.3}{}_{-0.3}{}^{}$	4.2${}_{}{}^{+0.3}{}_{-0.3}{}^{}$
328	2.12	1.91${}_{}{}^{-1.02}{}_{-0.12}{}^{}$	1.31${}_{}{}^{+0.11}{}_{-0.10}{}^{}$	1.57${}_{}{}^{+0.06}{}_{-0.06}{}^{}$	4.37	1.69	4.94	4.2${}_{}{}^{+0.4}{}_{-0.4}{}^{}$	3.4${}_{}{}^{+0.5}{}_{-0.4}{}^{}$	6.0${}_{}{}^{+0.4}{}_{-0.4}{}^{}$
367	2.07	0.42${}_{}{}^{+0.03}{}_{-0.02}{}^{}$	0.38${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.47${}_{}{}^{+0.06}{}_{-0.04}{}^{}$	0.14	0.13	0.07	3.7${}_{}{}^{+0.2}{}_{-0.3}{}^{}$	3.8${}_{}{}^{+0.3}{}_{-0.3}{}^{}$	2.2${}_{}{}^{+0.2}{}_{-0.2}{}^{}$
399	1.40	0.24${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.23${}_{}{}^{+0.03}{}_{-0.03}{}^{}$	0.32${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.38	0.13	0.32	2.9${}_{}{}^{+0.3}{}_{-0.3}{}^{}$	1.5${}_{}{}^{+0.3}{}_{-0.3}{}^{}$	2.6${}_{}{}^{+0.2}{}_{-0.2}{}^{}$
551	1.84	1.56${}_{}{}^{+0.11}{}_{-0.10}{}^{}$	1.62${}_{}{}^{+0.15}{}_{-0.14}{}^{}$	1.54${}_{}{}^{+0.09}{}_{-0.09}{}^{}$	4.02	2.31	3.89	2.0${}_{}{}^{+0.2}{}_{-0.2}{}^{}$	1.5${}_{}{}^{+0.2}{}_{-0.2}{}^{}$	1.9${}_{}{}^{+0.2}{}_{-0.2}{}^{}$
621	2.11	0.18${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.15${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.17${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.40	0.32	0.30	2.3${}_{}{}^{+0.2}{}_{-0.2}{}^{}$	2.0${}_{}{}^{+0.3}{}_{-0.3}{}^{}$	1.8${}_{}{}^{+0.2}{}_{-0.2}{}^{}$
658	1.67	0.15${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.18${}_{}{}^{+0.03}{}_{-0.03}{}^{}$	0.15${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.54	0.57	0.52	2.0${}_{}{}^{+0.2}{}_{-0.2}{}^{}$	2.0${}_{}{}^{+0.2}{}_{-0.2}{}^{}$	2.0${}_{}{}^{+0.2}{}_{-0.2}{}^{}$
733	1.84	0.21${}_{}{}^{+0.04}{}_{-0.03}{}^{}$	0.17${}_{}{}^{+0.03}{}_{-0.03}{}^{}$	0.12${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.94	0.96	0.92	1.8${}_{}{}^{+0.2}{}_{-0.2}{}^{}$	2.0${}_{}{}^{+0.2}{}_{-0.2}{}^{}$	2.2${}_{}{}^{+0.2}{}_{-0.2}{}^{}$
752	2.24	0.15${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.15${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.11${}_{}{}^{+0.01}{}_{-0.01}{}^{}$	0.28	0.33	0.25	4.6${}_{}{}^{+0.7}{}_{-0.7}{}^{}$	5.4${}_{}{}^{+0.7}{}_{-0.7}{}^{}$	5.0${}_{}{}^{+0.5}{}_{-0.5}{}^{}$
785	1.97	0.54${}_{}{}^{+0.05}{}_{-0.05}{}^{}$	0.56${}_{}{}^{+0.05}{}_{-0.04}{}^{}$	0.53${}_{}{}^{+0.03}{}_{-0.02}{}^{}$	4.73	5.97	9.26	16.2${}_{}{}^{+2.8}{}_{-2.7}{}^{}$	18.8${}_{}{}^{+2.1}{}_{-2.1}{}^{}$	29.7${}_{}{}^{+1.8}{}_{-1.8}{}^{}$
818	1.97	0.14${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.11${}_{}{}^{+0.03}{}_{-0.03}{}^{}$	0.19${}_{}{}^{+0.03}{}_{-0.02}{}^{}$	0.61	0.46	0.93	1.1${}_{}{}^{+0.1}{}_{-0.1}{}^{}$	0.9${}_{}{}^{+0.1}{}_{-0.1}{}^{}$	1.5${}_{}{}^{+0.1}{}_{-0.1}{}^{}$
840	1.60	0.58${}_{}{}^{+0.05}{}_{-0.04}{}^{}$	0.59${}_{}{}^{+0.04}{}_{-0.04}{}^{}$	0.40${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.53	0.49	0.59	3.9${}_{}{}^{+0.3}{}_{-0.3}{}^{}$	3.8${}_{}{}^{+0.3}{}_{-0.3}{}^{}$	5.1${}_{}{}^{+0.3}{}_{-0.3}{}^{}$
846	1.56	0.14${}_{}{}^{+0.03}{}_{-0.02}{}^{}$	0.11${}_{}{}^{+0.03}{}_{-0.02}{}^{}$	0.15${}_{}{}^{+0.03}{}_{-0.03}{}^{}$	0.84	0.86	0.90	2.1${}_{}{}^{+0.2}{}_{-0.2}{}^{}$	2.2${}_{}{}^{+0.2}{}_{-0.2}{}^{}$	2.2${}_{}{}^{+0.2}{}_{-0.2}{}^{}$
CDF-N
66	1.92	0.34${}_{}{}^{+0.09}{}_{-0.06}{}^{}$	0.24${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.29${}_{}{}^{+0.03}{}_{-0.03}{}^{}$	1.19	1.49	1.10	7.6${}_{}{}^{+1.0}{}_{-1.0}{}^{}$	12.7${}_{}{}^{+0.7}{}_{-0.7}{}^{}$	8.1${}_{}{}^{+0.7}{}_{-0.7}{}^{}$
143	1.46	0.15${}_{}{}^{+0.02}{}_{-0.02}{}^{}$	0.15${}_{}{}^{+0.01}{}_{-0.01}{}^{}$	0.11${}_{}{}^{+0.01}{}_{-0.01}{}^{}$	2.07	2.79	2.42	10.2${}_{}{}^{+0.9}{}_{-0.9}{}^{}$	13.6${}_{}{}^{+0.6}{}_{-0.7}{}^{}$	13.1${}_{}{}^{+0.7}{}_{-0.7}{}^{}$

Table 9. Table 9: The Chi-squared value χ 2 superscript 𝜒 2 \chi^{2} , the Δ ​ AIC Δ AIC \rm\Delta AIC  value Δ Δ \Delta , and the normalized excess variance σ 2 superscript 𝜎 2 \sigma^{2} of the observed 0.5–7 keV flux, intrinsic rest-frame L 0.5 − 7 ​ keV subscript 𝐿 0.5 7 keV L\rm_{0.5-7\,keV} and N H subscript 𝑁 H N_{\rm H} .

XID	${\chi }_{flux}^2$	${\chi }_{L_X}^2$	${\Delta }_{L_X}$	${\chi }_{N_H}^2$	${\Delta }_{N_H}$	${\sigma }_{flux}^2$	${\sigma }_{L_X}^2$	${\sigma }_{N_H}^2$
CDF-S 4 Epochs
49	53.3	21.9	11.9	17.3	11.3	0.074	0.056	0.073
81	69.1	23.7	13.9	3.1	$-$3.2	0.058	0.037	0
98	3.6	2.9	$-$3.1	0.4	$-$5.6	0	0	0
186	13.4	9.6	2.6	2.5	$-$3.9	0.007	0.009	0
214	2.9	0.4	$-$6.0	3.0	$-$3.3	0	0	0
236	7.0	2.7	$-$3.6	4.4	$-$2.1	0.008	0	0
308	86.7	44.9	32.4	1.4	$-$4.7	0.066	0.054	0
458	21.7	28.4	13.5	7.5	2.8	0.024	0.056	0.013
746	6.8	7.7	2.0	7.8	1.9	0.004	0.010	0.003
876	12.0	5.5	$-$0.9	0.4	$-$5.6	0.005	0.002	0
933	10.3	10.3	2.2	4.9	$-$1.4	0.002	0.004	0.005
760	2.6	3.0	$-$3.1	5.6	$-$1.7	0	0	0
CDF-S 3 Epochs
63	34.2	7.2	2.8	65.2	29.45	0.054	0.026	0.147
73	7.7	5.7	0.7	0.9	$-$3.27	0.029	0.032	0
240	0.6	1.6	$-$0.9	1.6	$-$0.52	0	0	0
249	31.0	2.2	$-$2.9	130.8	14.38	0.086	0	0.150
328	20.8	18.3	8.5	2.4	1.14	0.050	0.095	0
367	33.8	9.6	1.4	1.4	$-$0.38	0.044	0.046	0
399	13.5	14.6	8.3	2.6	$-$1.31	0.051	0.118	0
551	3.6	2.4	$-$2.2	0.1	$-$3.95	0.004	0	0
658	0.1	0.2	$-$3.8	0.8	$-$3.14	0	0	0
752	0.7	1.9	$-$2.3	6.8	1.16	0	0	0.002
785	26.0	10.6	4.3	0.2	$-$3.87	0.062	0.057	0
840	12.0	1.7	1.1	16.1	9.59	0.013	0	0.012
846	0.3	0.1	$-$3.9	1.1	$-$2.96	0	0	0
621	2.4	2.0	$-$2.0	1.3	$-$2.83	0	0	0
91	0.7	1.0	$-$3.1	1.1	$-$2.92	0	0	0
818	13.9	12.8	8.1	3.2	$-$0.7	0.035	0.067	0.004
733	2.0	0.04	-3.9	4.9	0.2	0	0	0.012
CDF-N 3 Epochs
66	25.8	4.2	$-$0.4	8.6	0.29	0.050	0	0
143	11.8	6.8	2.4	8.8	4.28	0.012	0.008	0.011

Equations

$${\rm HR}=\left\{\begin{array}[]{lr}-0.240z^{3}+0.230z^{2}+0.680z-0.200,\\ \ (z\ \leq\ 0.8);\\ -0.004z^{3}+0.077z^{2}-0.530z+0.740,\\ \ (z\ >\ 0.8).\end{array}\right.$$

$$F(N_{\mathrm{H}})dN{}_{\mathrm{H}}=f(N_{\mathrm{H}})d{N}_{\mathrm{H}}{\int }_0^{z_{\max }}\frac{dV}{dz}dz$$

$$\times {\int }_{L_{\mathrm{cut}(N_{\mathrm{H}},z)}}^{L_{\max }}N(L_X,z)d{L}_X,$$

$$f=\frac{{\int }_{L_{\mathrm{cut}(N_{\mathrm{H}},z)}}^{L_{\max }}N(L_{\mathrm{X}},z)d{L}_{\mathrm{X}}}{{\int }_{L_{\mathrm{cut}(23,\mathrm{z})}}^{L_{\max }}N(L_{\mathrm{X}},z)d{L}_{\mathrm{X}}}.$$

$$f=\left\{\begin{array}[]{lr}1,&{\rm if\ 23.0\ \rm cm^{-2}\leq log\,}N_{\rm H}<23.25\ \rm cm^{-2},\\ <1,&{\rm if\ log\,}N_{\rm H}\geq 23.25\ \rm cm^{-2}.\\ \end{array}\right.$$

$$f_{\mathrm{CH}}=\frac{N(\log N{}_{\mathrm{H}}\ge 24)}{N(\log N{}_{\mathrm{H}}\ge 23)},$$

$${\chi }^2=\sum _{i=1}^N\frac{(x_i-\overline{x}{)}^2}{{\sigma }_{x,i}^2}$$

$$\text{σnxv2}=\frac{1}{N⟨x{⟩}^2}\sum _{i=1}^N[(x_i-⟨x⟩{)}^2-(\delta x_i{)}^2],$$

$$s_{\mathrm{D}}^2=\frac{1}{N-1}\sum _{i=1}^N[(x_i-⟨x⟩{)}^2-(\delta x_i{)}^2-\text{σnxv2}⟨x{⟩}^2{]}^2.$$

Full text

galaxies: active — galaxies: nuclei — quasars: general — X-rays: galaxies

1 introduction

Highly obscured active galactic nuclei (AGNs), which are defined as AGNs with hydrogen column density ($N_{\rm H}$) larger than $\rm 10^{23}\ cm^{-2}$, are believed to represent a crucial phase of active galaxies. According to our knowledge of co-evolution of supermassive black holes (SMBHs) and their host galaxies (for reviews, see, e.g., Alexander & Hickox, 2012; Kormendy & Ho, 2013) and the hierarchical galaxy formation model, AGN activity may be triggered in a dust enshrouded environment, into which gas inflows due to either internal (e.g., disk instabilities; e.g., Hopkins & Hernquist, 2006) or external (e.g., major mergers; e.g., Di Matteo et al., 2005) processes both fuel and obscure the SMBH accretion, resulting in short-lived heavily obscured AGNs (e.g., Fiore et al., 2012; Morganti, 2017). The subsequent AGN feedback process may blow out the obscuring material and leave out an unobscured optically bright quasar. Compared with unobscured AGNs, AGNs in the highly obscured phase tend to have smaller BH masses, higher Eddington ratios ($\rm\lambda_{Edd}$; e.g., Lanzuisi et al., 2015) and larger merger fractions (e.g., Kocevski et al., 2015; Ricci et al., 2017a), which may indicate a fast growth state of central SMBHs (e.g., Goulding et al., 2011). Moreover, the cosmic X-ray background (CXB) synthesis models also require a sizable population of highly obscured AGNs, or even Compton-thick (CT) AGNs ($N_{\rm H}$ $\rm\gtrsim 10^{24}\ cm^{-2}$; see, e.g., Comastri 2004; Xue 2017; Hickox & Alexander 2018 for reviews), to reproduce the peak of CXB at 20$-$30 keV (e.g., Gilli et al., 2007, but see Treister et al. 2009). Therefore, the study of highly obscured AGNs across cosmic epochs is vital for our understanding of the AGN triggering mechanism, SMBH growth, AGN environment and the origin of CXB.

Thanks to the powerful penetrability of high energy X-ray photons, X-ray observations provide a great window to uncover the mysterious veil of these heavily obscured sources that are likely missed in optical surveys. In the past twenty years or so, the deep X-ray surveys 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., Baumgartner et al., 2013) and NuSTAR (e.g., Lansbury et al., 2017) have provided relatively unbiased AGN samples thanks to their unprecedented depths and sensitivities, which allow us to identify a significant population of heavily obscured AGNs (e.g., Risaliti et al., 1999; Brightman & Ueda, 2012; Ricci et al., 2015) using either X-ray color, spectral analysis and/or stacking technique (e.g., Alexander et al., 2011; Iwasawa et al., 2012; Georgantopoulos 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 heavily obscured systems, such as MIR excess (e.g., Daddi, 2007; Alexander et al., 2008; Luo et al., 2011), high 24 $\mu m$ to optical flux ratio (e.g., Fiore et al., 2009) and high X-ray to optical flux ratio (e.g., Fiore et al., 2003).

Among various methods, X-ray spectroscopy provides the most direct and unambiguous way to measure the column density of the obscuring materials. Several previous studies have focused on deriving the intrinsic $N_{\rm H}$ distribution corrected for the survey biases. Tozzi et al. (2006) presented a $N_{\rm H}$ distribution that has an approximately log-normal shape peaking at \sim$$\rm 10^{23}\ cm^{-2} and with an excess at \sim$$\rm 10^{24}\ cm^{-2}. Liu et al. (2017) (hereafter L17) reported a similar $N_{\rm H}$ distribution that peaks in a higher $N_{\rm H}$ range, due to the inclusion of more sources with low X-ray luminosities ($L_{\rm X}$) and high redshifts than that of Tozzi et al. (2006), which are expected to have relatively high $N_{\rm H}$ values. However, both works only focused on bright AGNs, and neither was dedicated to or optimized for investigating highly obscured sources. In particular, L17 excluded CT AGNs in their work and only focused on the Compton-thin population. Hence the absorption distribution and evolution of the most deeply buried AGNs are still unclear, especially at high redshifts (Vito et al., 2018). Therefore, unveiling the apparently faint, CT regime using the deepest X-ray survey data is indispensable for fully understanding the entire AGN population.

There have also been several attempts to constrain the obscured AGN fraction and CT fraction on the basis of modeling 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 studies support the picture of evolved absorption with redshift and 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 $\sim$50%, largely due to the limited sample sizes, the use of different X-ray spectral models, the unknown contribution from Compton reflection (Treister et al., 2009; Akylas et al., 2016) and the relatively poor quality of X-ray spectra at high $N_{\rm H}$ circumstances. Therefore, deeper X-ray observations as well as the physically appropriate spectral models for highly obscured AGNs are needed to robustly characterize the obscuration properties.

Among the X-ray surveys, the Chandra Deep Fields (CDFs) surveys (see, e.g., Brandt & Alexander 2015 and Xue 2017 for reviews), which consist of the $\rm 7\ Ms$ Chandra Deep Field-South survey (CDF-S; Luo et al., 2017) and the $\rm 2\ Ms$ Chandra Deep Field-North survey (CDF-N; Xue et al., 2016) along with the $\rm 250\ ks$ Extended Chandra Deep Field-South survey (E-CDF-S; Xue et al., 2016), provide us the most promising data to study highly obscured AGNs. In particular, the $\rm 7\ Ms$ CDF-S, which is the deepest X-ray survey to date, significantly improves the count statistics that allow us to extract high-quality X-ray spectra, detect more faint, highly obscured sources and perform more robust spectral analyses compared with previous 4 Ms analyses (e.g., Brightman et al., 2014). Furthermore, recent works suggest that the power spectral density (PSD) break frequency of AGN light curves might be related to $N_{\rm H}$ variability (Zhang et al., 2017; González-Martín, 2018), which makes obscuration a very important factor in investigating the driving mechanism of AGN variability. Benefiting from the very long timespan of the $\rm 7\ Ms$ CDF-S data (16.4 years in the observed frame), we are able to, for the first time, quantify the detailed variability behavior for a large, dedicated sample of highly obscured AGNs, in order to better understand the location of the obscuring materials and their contribution to AGN variability (also see Yang et al., 2016, Y16 hereafter).

In this study, we construct the largest dedicated highly obscured AGN sample in the deepest Chandra surveys which enables us to extend the studies of deeply buried sources to lower luminosities and higher redshifts in great details, and present systematic X-ray spectral and long-term variability analyses to study their evolution and physical properties. This paper is organized as follows. We describe our data reduction procedure and sample selection in Section 2. In Sections 3 and 4, we present detailed X-ray spectral analyses of our sample, focusing on the column density and luminosity distributions of highly obscured AGNs as well as their relations; the number counts of CT AGNs and the constraint on CXB models. The reprocessed components, the covering factor of the obscuring materials and their indications to AGN structures are discussed in Section 5. In Section 6, by correcting for several observational biases, we constrain the intrinsic $N_{\rm H}$ distribution representative for the highly obscured AGN population and study its evolution across cosmic time. In Section 7, we select a subsample of highly obscured AGNs that have largest counts available and perform detailed long-term X-ray variability analyses, in order to find out the variable fraction as well as the main driven mechanism of variability.

Throughout this paper, we adopt a Galactic column density of $N\rm_{H}=8.8\times 10^{19}\ cm^{-2}$ for the CDF-S, and $N\rm_{H}=1.6\times 10^{20}\ cm^{-2}$ for the CDF-N, respectively (Stark et al., 1992). We adopt cosmological parameters of $\rm H_{0}=70.0\ km\ s^{-1}\ Mpc^{-1}$, $\rm\Omega_{M}=0.30$, and $\rm\Omega_{\Lambda}=0.70$. All given errors are at 1$\sigma$ confidence level unless otherwise stated.

2 data reduction and sample selection

This work is based on the Chandra data in the $\rm 7\ Ms$ CDF-S main-source catalog (Luo et al., 2017) and the $\rm 2\ Ms$ CDF-N main-source catalog (Xue et al., 2016). Since the exposure in the E-CDF-S is $\sim$8–28 times shallower than in the CDF-N and CDF-S, we exclude the E-CDF-S data in this work. The redshift information for each source is adopted from the two catalogs that presented a preferred redshift value given by the ZFINAL keyword. The ZFINAL redshift values were collected and selected based on various published catalogs and followed a general preference order of spectroscopic redshift over photometric redshift (see Section 4.3 of Luo et al. 2017 and Section 2.3.4 of Xue et al. 2016 for details). The redshift (spectroscopic redshift) completenesses in the CDF-S and the CDF-N main source catalogs are 99.4% (64.8%) and 95.2% (52.4%), respectively; and the corresponding mean $1\sigma$ photometric redshift errors are about 0.21 and 0.17, respectively.

The source spectra from individual observations were extracted using the ACIS Extract (AE) software package (Broos et al., 2010). AE generates the point-spread function (PSF) model based on the MARX ray-tracing simulator and constructs a polygonal extraction region that corresponds to an encircled-energy fraction of $\sim$90%. For crowded sources, AE adopts smaller extraction regions to avoid overlapping polygonal regions. The background spectra were extracted using the BETTER_BACKGROUNDS algorithm. The most significant aspect of the above photometry and spectral extraction procedure, compared to the widely used circular-aperture extraction, is that it can obtain photometry and spectra as accurate as possible and remove the contamination from neighboring sources to faint sources to the greatest extent (see Section 3.2 of Xue et al. 2011 for details). This is extremely important for our work since we are dealing with the highly obscured AGNs that are generally fainter and expected to have limited counts due to significant obscuration.

The spectra eventually used in this work are the merged spectra for which all the individual observations were matched to a common $K_{s}$-band astrometric frame (see Section 2.2.1 in Xue et al. 2016 ) and stacked using the MERGE_OBSERVATIONS algorithm in AE. The corresponding response matrix files (rmf) and ancillary response files (arf) were also generated and combined during this stage.

We construct our sample by selecting sources which (1) are classified as AGN (TYPE=AGN; see Section 4.5 in Luo et al. 2017 for details); (2) have 0.5–7 keV net counts $>20$ (FB$\_$COUNTS $>$ 20) to allow basic X-ray spectral fitting; and (3) have redshift measurements (ZFINAL $!=-1$) in the main catalogs. The resulting full sample consists of 1152 sources, with 660 from the CDF-S and 492 from the CDF-N, respectively. The counts and redshift distributions for the full sample are shown 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 have spectroscopic redshifts.

3 X-ray Spectral Analysis

3.1 Spectral Fitting Models

We use MYTorus-based models (Murphy & Yaqoob, 2009) to fit the observed-frame 0.5–7 keV spectra of the full sample in order to identify heavily obscured sources. Due to limited counts, we do not bin the spectra because it may lose some key information of the sources. We use the Cash statistic (Cstat in XSPEC; Cash, 1979) as our spectral fitting statistic. Cstat has a similar probability distribution to $\chi^{2}$ statistics and has been proved to be more appropriate in the low-counts regime. Since Cstat is not appropriate for the background-subtracted spectra, we simultaneously fit the source and background spectra, with the latter (with 1642 median counts) being fit with the $cplinear$ model (Broos et al., 2010) that properly describes the observed Chandra background by a continuous piecewise-linear function in ten energy segments (an example is shown in Figure 3). In this way, we are able to maximize the usage of information relevant to the sources.

We adopt two models to fit the source spectra with different components and degrees of freedom (d.o.f) in order to find a statistically robust best-fit one:

The MYTorus baseline model: $phabs\times{\rm(MYTZ}\times\\ zpow+f_{\rm ref}\times{\rm MYTS}+f_{\rm ref}\times gsmooth{\rm(MYTL)})$. 2. 2.

The soft-excess model: phabs\times{\rm(MYTZ}\times\\ zpow+f_{\rm ref}\times{\rm MYTS}+f_{\rm ref}\times gsmooth{\rm(MYTL)}+\textit{f_{\rm exs}}\times zpow).

These models include all the typical spectral features found in highly obscured AGNs. The $phabs$ component models the Galactic photoelectric absorption. The MYTZ$\times zpow$ term represents the zeroth-order transmitted power-law continuum across the torus that takes into account both photoelectric absorption and Compton scattering processes. MYTS and $gsmooth\,$(MYTL) stand for the reflection component and the broadened Fe $\rm K{\footnotesize\alpha}$, $\rm K{\footnotesize\beta}$ fluorescent emission lines, respectively. A second power law is added to represent the soft excess component often found in AGN spectra, possibly originating from the zeroth-order continuum being scattered by the extended Compton-thin (CN) materials in obscured AGNs (Guainazzi et al., 2005; Bianchi et al., 2006; Corral et al., 2011). During spectral fitting, the inclination angle $\theta$ is fixed at $75^{\circ}$, which is the average value within a range of $60^{\circ}$ to $90^{\circ}$ where the torus intercepts our line of sight (l.o.s). The nominal normalizations of MYTS, MYTL and the second power law are set to be the same as that of the first power law, and we use the constants $f_{\rm ref}$ and $f_{\rm exs}$ to represent the real normalizations. $f_{\rm exs}$ is allowed to vary between $0-0.1$, and we assume a 200 keV high-energy cutoff throughout. One thing to mention is that the best-fit $N_{\rm H}$ values in MYTorus represent the equatorial values, and we always use the l.o.s values $N_{{\rm H,l.o.s}}=(1-4\,\cos^{2}\theta)^{\frac{1}{2}}N\rm_{H}$ in this work.

3.2 Model Selection

Several previous works suggested that the inclusion of the soft excess and reflection components in the spectral models is crucial for us to correctly estimate $N_{\rm H}$ of highly obscured sources (e.g., Brightman et al., 2014; Lanzuisi et al., 2015). However, the low counts of many sources do not allow us to apply complex models with free parameters since it could lead to large degeneracies. Therefore, we choose the model components and the parameter spaces according to the following criteria:

We fix $\Gamma$ 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, we set $\Gamma$ free to find a best-fit value. 2. 2.

For all sources, we first fit the spectra with a free $f_{\rm ref}$. If $f_{\rm ref}$ is less than $10^{-5}$, we then fix it to $10^{-8}$ that is an arbitrary value set as the lower limit for $f_{\rm ref}$ (i.e., indicating a negligible reflection component); otherwise, we fix $f_{\rm ref}$ at 1 that is the default value adopted in MYTorus. 3. 3.

To determine whether we should add a soft excess component, we compare the Cstat between models with and without considering soft excess emission. The best-fit model is chosen to be the one that has a statistically robust low Cstat value. More specifically, adding a new soft excess component should improve the Cstat at least for $\Delta C=3.84$ with 1 more d.o.f. This criterion is based on the fact that $\Delta C$ approximately follows the $\chi^{2}$ distribution, thus $\Delta C=3.84$ is roughly consistent with a $\textgreater$ 95% confidence level (Tozzi et al., 2006; Liu et al., 2017). 4. 4.

In order to avoid extremely untypical $\Gamma$ caused by the degeneracy between $\Gamma$ and $N_{\rm H}$, we re-fit the spectra of those sources with $\Gamma$ pegged at 1.4 (i.e., the lowest value permitted by MYTorus) or $\Gamma$ $\textgreater$ 2.4 (i.e., the typical maximum photon index for high-$\rm\lambda_{Edd}$ AGNs; e.g., Wang et al., 2004; Fanali et al., 2013) by fixing their $\Gamma$ at 1.8. If $\Delta C=C_{\rm fix}-C_{\rm free}$ $\textgreater$ 3.84, we adopt the free $\Gamma$ value; otherwise, we adopt $\Gamma$=1.8.

From now on, for convenience, we refer model 1 (2) with negligible $f_{\rm ref}$ as model A (B) and model 1 (2) with $f_{\rm ref}$ fixed at 1.0 as model C (D), as summarized in Table 1. In order to ensure the consistency of the models used for subsequent bias corrections (see Section 6) as much as possible, we have to assume fixed parameters in case of low counts. We justify our models in Appendix A by validating that the usage of fixed parameters does not significantly affect our results. We also compare our fitting results with several previous works (see Section 4.1) and those obtained from the Borus model (Baloković et al. 2018; see Section 5.3) and find good consistency.

We note that the simple absorbed power-law model (e.g., $phabs\times zwabs\times zpow$; hereafter model $z$) has been widely used in the literature to obtain rough estimates of AGN parameters. However, such a model is not appropriate for highly obscured AGNs since $zwabs$ only models photoelectric absorption and does not take into account the Compton scattering process, which is particularly important in the high $N_{\rm H}$ regime. Therefore, we also fit the spectra using model $z$, aiming at directly testing how accurately such a simple model reproduces the main spectral parameters.

4 Spectral Fitting Results

On the basis of the spectral fitting results, we find that 39% (458/1152) of sources in the full sample are identified to be highly obscured AGNs (hereafter the highly obscured sample; HOS). The main fitting results for HOS using MYTorus models A–D are presented in Table 2 and the full version is available on line. The best-fit models for sources in different count bins and $N_{\rm H}$ ranges are summarized in Table 3. In the following analyses, we only consider the 436 sources that have the absorption-corrected, rest-frame 2–10 keV luminosity (calculated from the lumin command after deleting all the additional components and only keeping the $zpow$ component) $L_{\rm X}$ $\textgreater$ $\rm 10^{42}\ erg\ s^{-1}$ to avoid possible contamination from star-forming galaxies.

4.1 Comparisons With Previous Works and the Model of $phabs\times zwabs\times zpow$

Thanks to the increased exposure time, the improved data reduction procedure (see Table 1 in Xue et al. 2016 for a summary), the updated redshift measurements and the usage of physically appropriate spectral fitting models for highly obscured AGNs, we are able to extract higher-quality spectra and obtain more robust parameter constraints than previous works in the same field (e.g., Tozzi et al., 2006), especially for faint sources. Before performing further analyses, we first make direct comparisons with several works which presented spectral fitting parameters in the CDF-S to understand how different methods influence the spectral fitting results. The works used for comparisons are:

Liu et al. (2017) (L17): L17 performed an extensive X-ray spectral analysis for the bright sources (hard-band counts $>$ 80) in the 7 Ms CDF-S using the $wabs\times(plcabs\times power-law+zgauss+power-law+plcabs\times pexrav\times constant)$ model. The background spectrum was also fitted with the cplinear model. The redshift information used in the two works is the same. 2. 2.

Brightman et al. (2014) (B14): B14 presented X-ray spectral fitting results in the 4 Ms CDF-S using the BNtorus model (Brightman & Nandra, 2011). 112 sources are found common between the two works; 24 of them have redshifts different by more than 0.2 ($\Delta z>0.2$) compared with the values adopted in the updated 7 Ms catalog, and are neglected in the following comparison. 3. 3.

Buchner et al. (2015) (B15): The 4 Ms CDF-S spectra were analyzed using a physically motivated torus model and a Bayesian methodology to estimate spectral parameters. 114 sources are found common between the two works and 30 of them with $\Delta z>0.2$ are excluded.

As shown in Figure 4, our measured $N_{\rm H}$ and $L_{\rm X}$ values are in general agreement with previous works, despite that for individual sources, the usage of different model configurations and data may result in large discrepancies. The largest distinction happens between B14 and our work that ten highly obscured AGNs in our sample are reported to be unobscured in B14. To understand the discrepancy, we re-fit their spectra using the same model and method as described in B14 and the results again confirm their heavily obscured nature. Therefore, the large discrepancy could be due to the adopted data of different depths as well as the different data reduction and spectral extraction methods used in the two works.

The comparisons between the results obtained through MYTorus and model $z$ are shown in Figure 5. Note that since MYTorus does not allow $N_{\rm H}$ to vary below 10^{22}\ \textit{\rm cm^{-2}}, a number of unobscured sources are thus have best-fit \textit{N_{\rm H}}=10^{22}\ \textit{\rm cm^{-2}}. Therefore, we only consider sources with MYTorus \textit{{\rm log},N_{\rm H}}>22.2\ \textit{\rm cm^{-2}} in the comparison.

The $N_{\rm H}$ and the observed 0.5–7 keV flux derived by the two models are generally consistent. But at a given $N_{\rm H}$, the absorption-corrected intrinsic $\rm 0.5-7\ keV$ flux ratio and the intrinsic $L_{\rm X}$ ratio between the two models dramatically increase with $N_{\rm H}$. It is obvious that model $z$ underestimates the intrinsic luminosity due to neglect of the Compton scattering process, and the discrepancy is already evident at \textit{N_{\rm H}}\sim\textit{\rm 10^{23}\ cm^{-2}} (see also Burlon et al., 2011). Therefore such models with only photoelectric absorption taken into consideration (e.g., $zwabs$, $zphabs$, $ztbabs$) must be used with caution in the highly obscured regime.

4.2 Photon Index Distribution

The photon index distribution for the 62 highly obscured sources with free-$\Gamma$ during the fitting process is shown in Figure 6. The mean value of the distribution is $1.82\pm 0.04$, in agreement with previous X-ray spectral analyses in the CDF-S (e.g., Tozzi et al., 2006; Liu et al., 2017) and COSMOS (e.g., Lanzuisi et al., 2013). Note that our $\Gamma$ distribution peaks at $\Gamma$ = 1.4 instead of the mean value 1.8. This is because the photon index is restricted within $1.4-2.6$ in MYTorus. If we use MYTZ alone which does not limit the $\Gamma$ range to fit the spectra, most sources with $\Gamma$ pegged at 1.4 will have smaller $\Gamma$, thus the distribution will appear more symmetric with a larger dispersion.

We also try to search the potential correlations among $\Gamma$, $N_{\rm H}$ and $L_{\rm X}$ using the Spearman rank correlation test. No correlation between $\Gamma$ and redshift (Spearman’s $\rho$ = 0.00, $p$-value = 0.99) is detected, indicating that the inner disk and corona structures have little evolution across cosmic time. The Spearman tests also suggest no significant correlation between $\Gamma$ and $N_{\rm H}$ ($L_{\rm X}$).

It is noteworthy that there are 8 sources with $\Gamma$ = 1.4 or $\Gamma$ $\textgreater$ 2.5 (see Table 4). Unlike some of the high-count sources that we artificially fix $\Gamma$ at 1.8 (see Section 3.2), these sources have significantly worse fits if we do not allow $\Gamma$ to vary freely (the average improvement of $\Delta C$ is 9.0 if we set $\Gamma$ free). Sources with very flat photon indexes are often considered to be reflection-dominated CT candidates and their extremely obscured nature may not be revealed by our best-fit $N_{\rm H}$ (Georgantopoulos et al., 2009, 2011). Moreover, since we are dealing with the stacked spectra here, sources which possessed large spectral variability may also exhibit untypical spectral shape, as might be the case for XID 249 (see Section 7.4). Additionally, we cannot rule out the possibility that the extreme photon indexes of some sources are wrongly measured due to insecure photometric redshifts.

4.3 Hardness Ratio versus Redshift

For highly obscured AGNs, the significant absorption and scattering of soft X-ray photons may lead to large hardness ratio (HR), thus the HR may be used as an indicator of obscuration level (e.g., Wang et al., 2004). In Figure 7 we show the observed HR (here we adopt the definition of HR as (H$-$S)/(H+S), where H stands for the observed-frame 2–7 keV count rate and S stands for the 0.5–2 keV count rate) as a function of redshift for highly obscured sources and CT candidates confirmed by the spectral fitting as well as less obscured AGNs (best-fit $N_{\rm H}$ $\textless$ $\rm 10^{23}\ cm^{-2}$) that contaminate our sample. Heavily obscured sources have significantly larger HR than less obscured sources as expected. We also calculate the effective photon index $\Gamma_{\rm eff}$ (obtained by fitting the spectra using XSPEC model $phabs\times zpow$ with $N_{\rm H}$ fixed at the Galactic value) for the HOS. 90% of them have $\Gamma_{\rm eff}$ $\textless$ 1.0, and the median $\Gamma_{\rm eff}$ is only $-0.38$ for CT AGNs, which again verify their heavily obscured nature.

The anti-correlation between HR and $z$ is due to that for high redshift sources, the observed soft band actually corresponds to a much harder band in the rest-frame and the photons may not be absorbed significantly. Note that the HRs of the most heavily obscured CT AGNs are not always the largest at a given redshift. This can be ascribed to the additional soft excess component softens the CT AGN spectra, and Compton-thin (CN) AGNs with intrinsically flat photon indexes are also plausible to have larger HRs than CT AGNs.

To test whether a simple HR value can be used to select heavily obscured AGNs, we calculate the evolution of HR as a function of redshift for typical X-ray spectral parameters. Since our purpose is to find the critical HR for highly obscured AGNs as a function of redshift, we simply consider the threshold condition: a source with \textit{N_{\rm H}}=\textit{\rm 10^{23}\ cm^{-2}}. Additionally, we set $f_{\rm ref}$ at 1.0, the high-energy cutoff at 200 keV and fix the photon index $\Gamma$ of the two power-laws at 2.0. The constant $f_{\rm exs}$ is set to 5% that is typical for AGNs with \textit{N_{\rm H}}\sim\textit{\rm 10^{23}\ cm^{-2}} (e.g., Brightman et al., 2014). The reason why we adopt \textit{\Gamma}=2.0 instead of the mean value 1.8 is that $\sim 96\%$ of the sources with a well-constrained photon index in L17 have \textit{\Gamma}\lesssim 2.0. Therefore, with this value we expect that our derived HR curve will be able to successfully select out most of highly obscured AGNs. The model is shown in Figure 2a and the derived model HRs at different reshifts are shown in green in Figure 7 (curve A).

It is worth noting that not all CT AGNs have larger HRs than CN AGNs. We also show the HR curve derived for a CT AGN with \textit{N_{\rm H}}=\textit{\rm 10^{24}\ cm^{-2}}. This time $f_{\rm exs}$ is chosen to be 1% due to the fact that $f_{\rm exs}$ decreases with increasing $N_{\rm H}$ (see Section 5.2). The model is shown in Figure 2b and the simulated HR curve is also shown in Figure 7 (curve B). It is clear that at low redshifts, even a very small soft excess fraction can easily dominate the soft X-ray spectra that could make CT AGNs even softer than CN AGNs. Therefore, to avoid missing a large number of high-$N_{\rm H}$ sources at low redshifts, we use curve B at $z\leq 0.8$ and curve A at $z>0.8$ as our selection curve in the following analysis (i.e., the combined curve).

For a source with given HR and redshift values, we consider it as a heavily obscured candidate if the observed HR lies above the combined HR curve. By applying this curve to our full sample, 480 sources will be selected as highly obscured candidates (hereafter the candidate sample). 80% (382/480) of sources in the candidate sample are indeed highly obscured as confirmed from spectral fitting, i.e., the accuracy is 80%. The parameter distributions for the selected candidates and confirmed sources are shown in Figure 8 (top panel). Only 20% of the candidate sources are contaminated by less obscured sources with most of them having $N_{\rm H}$ only slightly smaller than $\rm 10^{23}\ cm^{-2}$ and lying close to the boundary curve, as expected. The redshift, luminosity and soft excess fraction distributions for the selected candidates are very similar to sources which are confirmed to have \textit{N_{\rm H}}>\textit{\rm 10^{23}\ cm^{-2}}.

In addition, we define the completeness fraction as the ratio of the amount of highly obscured sources (confirmed) lie above the selection curve to the total highly obscured source amount (i.e., all red and blue points in Figure 7). The completeness fractions as a function of several parameters are shown in the second row in Figure 8. The average completeness fraction for the selection curve is 88% (382/436) and remains $>90\%$ for most high-redshift and high-luminosity bins. While at lower redshifts and larger $f_{\rm exs}$, the completeness fraction significantly drops to $<80\%$, since the soft excess dominates the observed spectrum which makes highly obscured sources lie below the selection curve.

We also show the result after changing the soft excess fraction for the model CT AGN to 5% (see curve C in Figure 7). This time the completeness fraction at low redshift is improved to 93%, but the accuracy dramatically drops to 68% since a large amount of less obscured AGNs are mistakenly included. As a trade-off, we choose to adopt the combined curve A and B as our selection curve. The functional form of this curve can be written as:

[TABLE]

We note that the conclusions in the following sections remain largely unchanged if we use the HR-selected sample instead, suggesting that a simple HR value can be used to identify highly obscured AGNs and select a representative highly obscured sample as long as the redshift and soft excess effects are properly taken into account while deriving the boundary curve.

4.4 Redshift and Luminosity Distributions

In Figure 9 we show the redshift distribution of the HOS in the right-bottom panel. 191 sources have spectroscopic redshifts and 245 sources have photometric redshifts. The peak appears at $z\sim 1-2.5$, which is known as the peak epoch of both star-formation and black hole accretion activities (e.g., Aird et al., 2015). The source amount slightly decreases down to $z\sim 0.5$, and significantly drops at $z<0.5$. Since Chandra only pointed at small patches of sky in these two surveys, the small number of sources detected at $z<0.5$ is mainly due to the small volume probed and the insufficient penetrability of the corresponding rest-frame energies at low redshifts; thus, the detection probability is severely limited in finding nearby highly obscured AGNs. The downward trend toward high redshifts is primarily due to the flux limits and may be partly due to the fact that photoelectric absorption in the soft band of mildly obscured (\textit{N_{\rm H}}\sim\textit{\rm 10^{23}\ cm^{-2}}) sources shifts out of the Chandra soft bandpass at $z\sim 2.5$, making it more difficult to constrain $N_{\rm H}$ (Vito et al., 2018).

The $L_{\rm X}$ distribution for the HOS is displayed in the top-left panel of Figure 9. The ${\rm log}\,L_{\rm X}$ ($L_{\rm X}$ in units of $\rm erg\ s^{-1}$) peaks in the range of 43.5–44.0 with the mean value of $43.65\pm 0.03$. The luminous sources with ${\rm log}\,L_{\rm X}$ $\textgreater$ 44.0 (44.5) account for 32% (11%) of the whole sample. CT AGNs, in particular, have even higher luminosities with the mean value of $44.10\pm 0.06$, because of the minimum detectable $L_{\rm X}$ for a flux-limited survey significantly increases with increasing $N_{\rm H}$ (see Figure 16). The correlations among $L_{\rm X}$, $N_{\rm H}$ and redshift are also plotted in Figure 9. The well-known Malmquist bias can be clearly seen in the bottom-left $L_{\rm X}$-$z$ corner and we will discuss the influence of this bias to our results in Section 6.4.

4.5 Observed $N_{\rm H}$ Distribution

The observed distributions of the best-fit $N_{\rm H}$ in six luminosity bins of the HOS are shown in Figure 10. The small amount of high $N_{\rm H}$ sources observed in the smaller $L_{\rm X}$ bins will be clearly illustrated in Section 6.2 (see Figure 16) that the minimum detectable luminosity significantly increases with $N_{\rm H}$. We show the normalized cumulative $N_{\rm H}$ distributions in Figure 11a. Sources with best-fit $N_{\rm H}$ $\textgreater$ $\rm 1.0\times 10^{24}\ cm^{-2}$ account for $\sim$25% (108/436) of the HOS. The distribution obtained in Brightman et al. (2014), which utilized the 4 Ms CDF-S data, is shown in the same plot for comparison. The K-S tests imply that the $N_{\rm H}$ distributions from the CDF-S and the CDF-N are consistent with $p$-value = 0.96, but are different from that in Brightman et al. (2014) with $p$-value $\ll 0.001$, due to that we identify more sources between ${\rm log}\,N_{\rm H}$ = 23.5 – 24.0 $\rm cm^{-2}$ while the distribution in Brightman et al. (2014) peaks at ${\rm log}\,N_{\rm H}$ = 23.0 – 23.5 $\rm cm^{-2}$.

Although we have significantly improved the photon statistics for each source using the deepest data, we should caution that most of CT AGNs discovered still have counts $\textless$ 200 and their $N_{\rm H}$ might be poorly constrained. Therefore we adopt a more conservative criterion that a source will be regarded as CT candidate only if it has the best-fit \textit{N_{\rm H}}>\textit{\rm 10^{24}\ cm^{-2}} and the 1$\sigma$ lower limit on $N_{\rm H}$ is constrained to be greater than $\rm 5\times 10^{23}\ cm^{-2}$. This criterion selects 102 sources, with 66 from CDF-S and 36 from CDF-N, accounting for $\sim$23% of the HOS, to be CT candidates. Note that there are some well-studied CT candidates in the literature that are not classified as CT AGNs in our sample due to their best-fit $N_{\rm H}$ $\textless$ $\rm 10^{24}\ cm^{-2}$, such as XID 328, XID 551 (XIDs 153, 202 in Comastri et al., 2011, respectively), XID 539 (XID 403 in Gilli et al., 2011) and XID 375 ($BzK$8608 in Feruglio et al., 2011) in the CDF-S. However, all these sources have best-fit $N_{\rm H}$ $\textgreater$ $8\times$ $\rm 10^{23}\ cm^{-2}$, which shows good consistency with previous studies, thus we confirm their heavily obscured nature using deeper Chandra data.

We note that at f_{\rm 2-10\,keV}>10^{-15}\ \textit{\rm erg\ cm^{-2}\ s^{-1}}, the number of our CT AGNs in the CDF-N (22 sources) is higher than Georgantopoulos et al. (2009) which presented 10 CT candidates in the same field. This may because of Georgantopoulos et al. (2009) mainly focused on searching reflection-dominated CT AGNs and possibly missed some transmission-dominated sources. We also fit the 22 CT candidates using the BNtorus model (Brightman & Nandra, 2011) and find that 20 sources have best-fit \textit{N_{\rm H}}>\textit{\rm 10^{24}\ cm^{-2}} and $f_{\rm 2-10\ keV}$ $\textgreater$ $10^{-15}$ $\rm erg\ cm^{-2}\ s^{-1}$, 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. The observed log$\,N-$ log$\,S$ is defined as the observed source amount divided by survey area. However, as we will discuss in Section 6.2, the sky coverage is not a constant value across the whole flux and $N_{\rm H}$ ranges, hence we assign a weighting factor $1/\omega_{\rm area}$ (see Section 6.2 for details) for each source while calculating their cumulative number counts. The corrected results are shown in Figure 11b, and several number counts measurements presented in previous works in the 4 Ms CDF-S (Brightman & Ueda, 2012), X-UDS (Kocevski et al., 2018) and COSMOS (Lanzuisi et al., 2018) fields are also displayed for comparison. Note that Brightman & Ueda (2012) reported the result in the 0.5–8 keV band and we convert the 0.5–8 keV flux into 2–10 keV by assuming a \textit{\Gamma}=-0.4 power-law which is the median effective photon index for our CT AGNs.

The corrected number counts for the CDF-S and CDF-N are generally consistent within the error bars at higher fluxes, despite that we obtain six CT AGNs in the CDF-N but only three in the CDF-S at $f_{\rm 2-10\ keV}$ $\textgreater$ 4\times 10^{-15}\ \textit{\rm erg\ cm^{-2}\ s^{-1}}, possibly caused by the low number statistics and the cosmic variance due to the small sky coverage of CDFs (e.g., Harrison et al., 2012). While at the faint end, both CDF-S and CDF-N number counts flatten due to the limited detection ability, and the CDF-N number counts are significantly smaller than the CDF-S, because the three times shallower exposure limits its detectability of faint CT sources (see Figure 11c). Our results are also in agreement with the 4 Ms CDF-S result at the faint end, the COSMOS result at the bright end, and the X-UDS result at moderate fluxes.

We also compare the observed log$\,N-$ log$\,S$ with several CXB model predictions in $z=0-5$ (Gilli et al., 2007; Ballantyne et al., 2011; Akylas et al., 2012; Ueda et al., 2014). For the Gilli et al. (2007)111http://www.bo.astro.it/~gilli/count.html/. We assume a high-$z$ declined LF (Vito et al., 2018)., Akylas et al. (2012)222http://indra.astro.noa.gr/xrb.html. We assume a 40% CT fraction and the default 4.5% reflection fraction. and Ueda et al. (2014)333http://www.kusastro.kyoto-u.ac.jp/ỹueda/xrb2014.html CXB models, the luminosity range used to derive the predicted number counts is \textit{{\rm log},L_{\rm X}}=42-45.5\ \textit{\rm erg\ s^{-1}}, similar to our sample; while for the Ballantyne et al. (2011) model, the predicted result is presented in \textit{{\rm log},L_{\rm X}}=41.5-48\ \textit{\rm erg\ s^{-1}}. As can be seen from Figure 11b, our observed log$\,N-$ log$\,S$ prefers the moderate CT number counts as predicted by Akylas et al. (2012) and Ueda et al. (2014), while other models more or less overestimate or underestimate the number counts.

In summary, the observed parameter distributions and relationships presented in this section provide a basic description of the highly obscured AGN population and are crucial for distinguishing various CXB models. However, these observed distributions, in particular, the observed $N_{\rm H}$ distribution, are influenced by several biases. We will discuss more details about these biases and reconstruct the intrinsic $N_{\rm H}$ distribution representative for the highly obscured AGN population in Section 6.

5 The Origin of Heavy X-ray Obscuration

5.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/radio counterpart positions presented in the X-ray source catalogs (CP_RA and CP_DEC) with the Silverman et al. (2010) E-CDF-S optical spectroscopic catalog using a 0.5*′′* 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 summarized in Table 5.

Most of the six sources have sufficient counts and reliable redshift measurements to constrain $N_{\rm H}$. Three sources have insecure spectroscopic redshifts although they are labeled as BLAGN, and the 7 Ms main catalog chooses to adopt the photometric redshift as ZFINAL for two of them. This could be due to that they have low S/N optical spectra or only one emission line is available for sources within particular redshift ranges, which makes it difficult to conclusively determine the line nature, thus giving an insecure redshift. For three sources with insecure redshifts, we set redshift as a free parameter in the spectral fitting and obtain the best-fit X-ray redshift and corresponding $N_{\rm H}$. All of them are still best-fitted by \textit{N_{\rm H}}>\textit{\rm 10^{23}\ cm^{-2}}. The redshift range in which the source will remain X-ray highly obscured is also listed. In general, the $N_{\rm H}$ estimates for the six sources are robust. This result suggests that the heavy X-ray absorption in a fraction of sources is largely a l.o.s effect caused by some compact clumpy clouds obscuring the central X-ray emitting region, but the global covering factor (CF) for the high-$N_{\rm H}$ materials is limited thus the BLR is not blocked; or maybe the heavy X-ray obscuration is produced by the BLR itself.

5.2 Soft Excess Fraction Dependences

There are 80 sources in our sample that require a second power law to fit the soft excess component. The origin of this component is still a puzzle. Many different models, e.g., warm Comptonization or blurred ionized reflection from the disk, partially ionized absorption in a wind, or power-law continuum in other directions scattered into the line of sight from large-scale Compton-thin matters, are proposed to explain the diverse situations found in different sources (e.g., Boissay et al., 2016). However, for highly obscured AGNs, the scattered mechanism is preferred since either blurred reflection or warm Comptonization from the accretion disk will be significantly attenuated by the obscuring materials.

We show the $f_{\rm exs}$, which represents the relative normalization of the soft component with respect to the intrinsic power law, as a function of $N_{\rm H}$ in Figure 12. We find a significant anti-correlation between the two parameters: Spearman’s $\rho=-0.66$, $p$-value = $\ll 0.001$. The mean $f_{\rm exs}$ is 3.7% for the total sample. The scattered fraction decreases from 6.1% to 4.3% and 1.3% for $N_{\rm H}$ $\textless$ 5 $\times$ $\rm 10^{23}\ cm^{-2}$, 5 $\times$ $\rm 10^{23}\ cm^{-2}$ $\textless$ $N_{\rm H}$ $\textless$  1.5 $\times$ $\rm 10^{24}\ cm^{-2}$ and $N_{\rm H}$ $\textgreater$ 1.5 $\times$ $\rm 10^{24}\ cm^{-2}$, respectively.

To verify that this anti-correlation is not caused by the parameter degeneracy that sources with a large soft excess fraction and a high $N_{\rm H}$ might be misclassified as low $N_{\rm H}$, low $f_{\rm exs}$ and low intrinsic luminosity sources owing to the low S/N, we assume that CT AGNs have exactly the same $f_{\rm exs}$ as CN AGNs (here we simply choose $f_{\rm exs}$ = 5%). Then we generate fake spectra using model D for a sample of simulated sources, which have redshifts and $N_{\rm H}$ uniformly distributed between $z=0-4$ with an interval of $\Delta z=0.5$ and \rm 23\ cm^{-2}<\textit{{\rm log},N_{\rm H}}<25\ cm^{-2} with an interval of \Delta\textit{{\rm log},N_{\rm H}}=0.25. For each source we simulate 100 spectra with 0.5–7 keV counts $\sim$ 50, 100 and 150, respectively. We then fit the fake spectra and perform Spearman rank correlation test for the output $N_{\rm H}$ and $f_{\rm exs}$. The Spearman’s $\rho\approx 0.1$ (insensitive to the simulated spectral counts), indicating no apparent correlation, thus the observed anti-correlation is intrinsic.

By assuming that the second power law originates from the scattered-back continuum, Brightman & Ueda (2012) found that the scattered fraction depends on the opening angle, i.e., sources with small opening angles have lower $f_{\rm exs}$. Moreover, Ueda et al. (2015) and Kawamuro et al. (2016) found that the [O III] and [O IV] to hard X-ray luminosity ratios are lower in low scattering fraction sources, inferring that the low $f_{\rm exs}$ AGNs are buried in small opening angle torus. Therefore, the anti-correlation between $f_{\rm exs}$ and $N_{\rm H}$ indicates that high $N_{\rm H}$ sources might preferentially reside in high CF toruses. This is also consistent with the results from Brightman & Ueda (2012) and Lanzuisi et al. (2015), who found a similar anti-correlation between $f_{\rm exs}$ and $N_{\rm H}$ and explained it as highly obscured AGNs being also heavily geometrically obscured (but see the discussion in the next section).

5.3 Reprocessed Components and Covering Factor

Unlike local CT AGNs in which the prominent neutral Fe lines and reflection hump (hereafter the reprocessed components) are prevalently detected in the high-quality hard X-ray spectra and can be served as unambiguous signatures of being Compton-thick (e.g., Tanimoto et al., 2018; Zhao et al., 2019; Marchesi et al., 2019), there are 41% of our highly obscured sources and 38% of the CT candidates having negligible $f_{\rm ref}$(see Table 3), respectively. It should be noted that the low detection rate of reprocessed components is not in conflict with their highly obscured nature. The poor S/N for low-count sources and the narrow Chandra spectral coverage make it challenging to detect narrow iron lines and to distinguish the transmission and reflection components, which may cause a significant overestimation of the fractions.

But we note that it is possible that highly obscured AGNs may have weak reprocessed components if the global CFs of high-$N_{\rm H}$ materials are low or the central AGN is geometrically fully buried by CT materials such that even the reflected photons cannot escape (e.g., Brightman et al., 2014). To distinguish different scenarios, we use the Borus model (Baloković et al., 2018)444In the Borus model, $\rm CF_{\rm tor}=1.0$ corresponds to a fully covered torus while $\rm CF_{\rm tor}=0.1$ represents a typical disk-like covering (see http://www.astro.caltech.edu/~mislavb/download/)., which allows CF to vary freely, to fit the spectra for all sources with counts $\textgreater$ 200. This time we treat photon index as a free parameter in order to make a direct comparison between the two models. The distribution of the derived $\rm CF_{\rm tor}$, which is defined as the cosine value of the opening angle $\theta_{\rm tor}$ measured from the symmetry axis towards the equatorial plane, is shown in Figure 13a. We also show the comparisons of the main spectral fitting parameters with those obtained from the MYTorus model and find good consistency.

For sources with negligible reprocessed components in the MYTorus model ($f_{\rm ref}\ll 1$), the Borus results are consistent with the weak torus scenario (\textit{\rm CF_{tor}}\approx 0.1), which might suggest that their heavy obscuration is simply a l.o.s effect (i.e., some high-column-density clouds on various scales along our sightline may obscure the compact X-ray emitter), without the necessity to invoke a strongly-buried nuclear environment. In contrast, sources with $f_{\rm ref}$ = 1.0 are mostly best fitted by a highly-covered model (0.8<\textit{\rm CF_{tor}}<1.0) that have torus covering angles ($90^{\circ}-\theta_{\rm tor}$) between $65^{\circ}-90^{\circ}$. We emphasize that though the l.o.s and the global torus $N_{\rm H}$ are linked in the spectral fitting, those fully-buried sources do not need to be covered by very high $N_{\rm H}$ materials in all directions, since the best-fit $N_{\rm H}$ will more likely converge to the values determined by the l.o.s component because the photoelectric absorption is a much stronger spectral feature. Moreover, we find that the average $\rm CF_{tor}$ for \textit{f_{\rm exs}}<0.05 sources is 0.60; while for \textit{f_{\rm exs}}>0.05 sources, the \overline{\textit{\rm CF_{tor}}} is 0.32. This might provide evidence for that the lower soft-excess fraction in some sources is caused by a geometrically more buried structure (see Section 5.2).

However, we caution here that the current S/N and spectral coverage do not allow us to constrain CF from X-ray spectral analysis. If we fix $\rm CF_{tor}$ at 0.1 for sources that are best fitted by a fully-buried model, we may still obtain acceptable fits with Cstat values only slightly increasing (an example is shown in Figure 13b). Conclusively disentangling several scenarios is beyond the scope of this work, and multi-wavelength data are needed to further shed light on the nature of the obscuring materials (Li et al., in preparation).

6 Intrinsic $N_{\rm H}$ distribution for highly obscured AGNs

In Section 4.5, we present the observed $N_{\rm H}$ distribution which is affected by several sample incompleteness. In order to derive the intrinsic $N_{\rm H}$ distribution, we should take into account the errors on best-fit $N_{\rm H}$ and correct for several survey biases (Liu et al., 2017). To perform such corrections, we restrict our analysis to the 394 sources with off-axis angle $\textless$ 10.0*′*, for the sake of avoiding large background contamination, extremely small sky coverage (see Section 6.2) and limited detectable fraction (see Section 6.4).

6.1 Errors on Best-fit $N_{\rm H}$

To consider the uncertainty of the spectral fitting, we perform a resampling procedure to the best-fit $N_{\rm H}$. Given the asymmetric errors, we assume that the errors on $N_{\rm H}$ obey the “half-gaussian” distribution. We first generate 1000 $N_{\rm H}$ for each source with the $\sigma$ of the half-gaussian distribution equals to the lower 1$\sigma$ error, and then generate another 1000 $N_{\rm H}$ with $\sigma$ equals to the upper 1$\sigma$ error. The mean value $\mu$ is set to the best-fit $N_{\rm H}$. The resampled $N_{\rm H}$ distribution (calculated by averaging the 2000 resampled distributions and also corrected for the sky coverage effect; see Section 6.2) is shown in the shaded region of Figure 18. It has an extended tail down to \textit{N_{\rm H}}<10^{23}\ \textit{\rm cm^{-2}} which contains only 4.0% of the resampled data, suggesting that even if considering the spectral fitting errors, most of our sources are still consistent with being highly obscured.

6.2 Sky Coverage Effect

Due to Chandra’s instrument features, the point-source PSF size increases and the effective exposure time dramatically decreases toward large off-axis angles. The sensitivities of detecting faint sources reduce prominently at the outskirt of the field, leading to small sky coverage while the observed flux is low. Therefore, a correction must be made.

We calculate the energy flux to count rate conversion factor (ECF) by assuming the soft-excess model (model D) with $\Gamma$, $f_{\rm ref}$ and $f_{\rm exs}$ fixed at 1.8, 1.0 and 1%, respectively. Combining ECF and the exposure map, we build a sensitivity map that represents the flux limits corresponding to the 20-count cut. Using this map, we calculate the sky coverage as a function of observed 0.5–7 keV flux, $N_{\rm H}$ and redshift as shown in Figure 14. Then we measure the sky coverage for each source based on their observed flux, $N_{\rm H}$ and redshift obtained from spectral fitting. We define $\omega_{\rm area}$ as the ratio between the source sky coverage and the maximum sky coverage of the two Chandra surveys (484.2 arcmin2 and 447.5 arcmin2 in the CDF-S and the CDF-N, respectively). To correct for the sky coverage effect, we simply weigh each source by $1/\omega_{\rm area}$ while resampling the observed $N_{\rm H}$ distribution. To avoid extremely large weights, we apply another cut for the 22 sources with sky coverage less than 100 $\rm arcmin^{2}$, setting their weighting factors to the median factor 1.3, since these sources lie close to the count cut and we cannot rule out the possibility that their extremely small sky coverage may be a result from inappropriate spectral modeling.

6.3 Eddington Bias

Considering the measurement error of net counts, sources with intrinsically low counts may exceed our count cut (20 photons), and sources with high counts may be missed in our sample. Since the number of faint and bright sources are not equal, this error leads to the Eddington bias. To correct for this bias, we convolve the cumulative count distribution for AGNs in the two Chandra survey catalogs with the count errors using Poisson sampling in order to obtain the intrinsic count distribution. We follow the ‘pseudo-deconvolved’ method proposed in L17 (see Section 5.1.3 of L17 for details), by shifting the observed curve leftward with the value equal to the displacement between the observed and convolved curves. We then obtain the deconvolved curve that represents the intrinsic cumulative count distribution.

The difference between the deconvolved and the observed curves can be treated as the number of sources missed or mis-included (depending on the adopted count cut) in our sample. As shown in Figure 15, at a 20-photon count cut, we miss 11.6 and 6.5 sources in the CDF-S and CDF-N, respectively. This is due to the fact that though sources with counts $\textless$ 20 may outnumber the bright sources, such faint sources are not detected in the source catalogs. By controlling the sample size (i.e., the cumulative source number in the $y$-axes) to be the same, we obtain the effective net counts of 22.4 and 21.1 for the CDF-S and CDF-N, respectively. Thus the 20-photon count cut will be replaced by the new effective count cut while performing other corrections.

6.4 Malmquist Bias

So far we have only obtained the corrected observed $N_{\rm H}$ distribution for our sample. In order to derive the intrinsic $N_{\rm H}$ distribution for the highly obscured AGN population, we should also consider the undetectable part of the Chandra survey, which is know as the Malmquist bias that sources with lower luminosities will be missed by a flux-limited survey. More specifically, given the count rate limit of the survey (the effective count cut in our situation), whether a source is detectable or not depends on its redshift, spectral shape, column density, intrinsic luminosity as well as its position in the field, since large off-axis angles lead to significant lower sensitivities.

To quantify this bias, we generate fake X-ray spectra to find out the minimum $L_{\rm X}$ ($L_{\rm cut}$, corresponding to our effective count cut) in different redshift and $N_{\rm H}$ bins, which can be detected at a specific position in the field. We divide the two CDFs into 10 annuli, with each corresponding to an interval of 1*′*. While running XSPEC simulations, we use the averaged exposure time, rmf file and arf file calculated from all sources in each annulus. Model D is used in the simulation with the parameters set as the same as in Section 6.2.

We show some of the detectable boundary curves in Figure 16. We define the detectable region as the area between the boundary curve and the maximum log $L_{\rm X}$ = 46 $\rm erg\ s^{-1}$. As we can clearly see, sources with high $N_{\rm H}$ and large off-axis angles have significantly small detectable regions. We follow Equation 3 in T06,

[TABLE]

where $f$($N_{\rm H}$) represents the intrinsic $N_{\rm H}$ distribution, $F$($N_{\rm H}$) represents the observed $N_{\rm H}$ distribution, and $N$($L_{\rm X}$, $z$) is the X-ray LF. The detectable source fraction for each $N_{\rm H}$ bin relative to the ${\rm log}\,N_{\rm H}$ = 23 $\rm cm^{-2}$ bin in a given redshift interval can be calculated as:

[TABLE]

Thus for the bin we choose,

[TABLE]

Note that T06 derived the above Equation 2 by assuming that $f(N_{\rm H},L_{\rm X},z)$, which represents the possibility of detecting a source with column density in the range of $N_{\rm H}$ to $N_{\rm H}$ + d$N_{\rm H}$ for a given $L_{X}$ and $z$, varying slowly as $L_{\rm X}$ and $z$ change. There is evidence that the obscured fraction of AGNs evolves with both luminosity and redshift. However, we do not consider such a complicated evolution and simply extract $f(N_{\rm H},L_{\rm X},z)$ from the integral. Therefore, by multiplying $1/f$ to the observed $N_{\rm H}$ distribution $F(N_{\rm H})\,dN_{\rm H}$ in each bin, sources with different $N_{\rm H}$ are corrected to cover the same observable space with respect to the ${\rm log}\,N_{\rm H}$ =$\,23-23.25\ \rm cm^{-2}$ bin and the Malmquist bias is corrected. To avoid the case of the corrected distribution being dominated by sources with extremely small detectable fractions, we apply a cut that while $f$ $\textless$ 0.25, the weighting factor $1/f$ will be simply set to 1/0.25. Adopting a lower (0.2) or higher (0.3) detectable fraction cut will lead to $\sim 2\%$ systematic difference of the final intrinsic CT-to-highly-obscured fraction.

Apparently, the correction depends on both spectral models and the detailed shape of X-ray LFs. Therefore we compare the results using different X-ray LFs. We combine the $z$ $\textless$ 3 LF from Miyaji et al. (2015) and $z$ $\textgreater$ 3 LF from Georgakakis et al. (2015) as our first LF model, which is the same as in L17. As for the second model, we adopt Aird et al. (2015) LFs for obscured AGNs (22 $\rm cm^{-2}$ $\textless$ $N_{\rm H}$ $\textless$ 24 $\rm cm^{-2}$). We do not use their LFs for CT AGNs, since they obtained a systematically lower CT fraction than other works and our result in Section 4.6 prefers higher CT number counts (but see the discussion in Section 6.4 of Aird et al. (2015) that this low CT fraction should not be ruled out) and the spectral model they used may not be appropriate for CT AGNs (see Section 4.1). The third LF we used is from Ueda et al. (2014). The three LF models are shown in Figure 17. We note that the adopted obscured AGN LFs may not be suitable for CT AGNs, but an extensive exploration of this issue is beyond the scope of this work.

6.5 Intrinsic $N_{\rm H}$ Distribution

Figure 18 shows the intrinsic $N_{\rm H}$ distributions in five redshift bins after correcting for all the aforementioned biases. We define the intrinsic CT-to-highly-obscured source fraction $f_{\rm CH}$ as

[TABLE]

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_{\rm H}$ distributions are quite different from the observed ones, especially in the CT regime. This highlights the importance of carefully correcting for the observational biases in order to uncover the underlying distributions. Though the intrinsic $N_{\rm H}$ distributions are different between high-redshift and low-redshift bins, the CT-to-highly-obscured fraction is in accordance with no evident redshift evolution given the uncertainties (i.e., $f_{\rm CH}\approx 0.52$ for all redshift bins; see Table 6). Note that Buchner et al. (2015) also claimed that the CT fraction is constant across cosmic time, but the CT fraction in their work is defined as the number of CT AGNs to the total AGN population. To avoid possible biases induced by the fact that different redshift bins are sensitive to different luminosity ranges, we also select two subsamples with 43.0\ \textit{\rm erg\ s^{-1}}<\textit{{\rm log},L_{\rm X}}<44.5\ \textit{\rm erg\ s^{-1}} at $z=1-2$ and $z=2-3$, respectively, which have similar $L_{\rm X}$ distributions (see Figure 9), and perform the same corrections. The results are listed in Table 6 and are still consistent with no redshift evolution.

Vito et al. (2018) presented the $N_{\rm H}$ distribution for the high-redshift ($z$ $\textgreater$ 3) AGNs also in the CDFs. Their analyses were based on modeling the X-ray spectra using the $wabs\times zwabs\times zpow$ model for the $z$ $\textgreater$ 3 sources by carefully taking into account the photometric redshift errors. The comparison of the two works is shown in the $z=3.0-5.0$ bin in Figure 18. Our result shows a higher CT fraction, which might be ascribed to the different redshift values adopted in the two works as well as the large weighting factors while we correcting for the undetectable part (see Section 6.4).

However, we note that all these aforementioned corrections are based on sources we do observe, including the correction of the undetectable part in a given $N_{\rm H}$ bin for which we have to make assumptions regarding the unknown AGN population. The series of corrections can be simply understood as multiplying an evolved weighting factor to the observed distribution, which have several limitations. First of all, as mentioned before, different redshift bins sample different luminosity ranges. Several authors have found that the obscured fraction changes with luminosity (Buchner et al., 2015; Lusso et al., 2013), thus the $N_{\rm H}$ distribution might be luminosity-dependent. However, in Equation 2, we ignore this effect to simplify the calculation. Second, for the extremely buried sources ($N_{\rm H}$ $\textgreater$ $\rm 10^{25}\ cm^{-2}$), or intrinsically very faint sources, the corrections cannot be made since there are hardly any such sources observed due to the limit of the survey and the restrictions imposed by our sample-selection count cut. Therefore, the part of the $N_{\rm H}$ distribution of sources that are under the detection limit is actually still missed in our final result. Furthermore, our identification method of CT AGNs is based on the absorption turnover imprinted in the spectra. For a source with $N_{\rm H}$ $\textgreater$ $\rm 10^{25}\ cm^{-2}$, the turnover occurs at rest-frame $\sim$ 20 keV. This implies that for low-redshift sources, the characteristic feature of CT AGNs may not be observable in the Chandra spectra of limited energy range, which is coupled with the uncertainties induced by the photometric redshifts, leading to great challenges in correctly measuring the obscuration level. It is plausible that $N_{\rm H}$ $\textgreater$ $\rm 10^{25}\ cm^{-2}$ sources may contribute significantly to heavily obscured AGN population (Risaliti et al., 1999), but such sources are missed in our sample and the corrections at this part is beyond our attainment. Therefore, the $N_{\rm H}$ distribution at the highest $N_{\rm H}$ bin displayed in Figure 18 is highly uncertain and incomplete so that the CT-to-highly-obscured fractions in Table 6 should be better treated as lower limits.

7 Spectral variability analyses

X-ray variability is a ubiquitous feature among AGNs that can provide useful information about AGN properties and relevant underlying physical processes (e.g., Vaughan et al., 2003; McHardy et al., 2006; González-Martín & Vaughan, 2012; Soldi et al., 2014; Zheng et al., 2017). In this section, we aim to study the main driving mechanism that causes the variability of highly obscured AGNs by investigating the variability of some main spectral parameters, such as the luminosity, column density and reflection strength.

To perform detailed long-term spectral variability analyses of highly obscured AGNs and expand the sample size, we select our sources from the HOS based on only one additional criterion, i.e., the net 0.5–7 keV counts in the total $\rm 7\ Ms$ exposure in the CDF-S or in the total $\rm 2\ Ms$ exposure in the CDF-N should be larger than 700 or 900, respectively. We bin data from neighboring observations as one epoch to enhance the S/N. For sources in the CDF-S that have counts $>$ 1500, we bin the data into 4 epochs; for sources having 700 $<$ counts $<$ 1500, we bin the data into 3 epochs. The data for all the sources in the CDF-N are binned into 3 epochs. Detailed binning information can be found in Table 7. These binning strategies and the adopted count criteria can roughly ensure that the average counts in each epoch are at least $\sim 200$. Thirty-one sources, which satisfy our selection criteria, are selected for subsequent variability analyses. The redshift distribution of this variability sample is shown in Figure 9 and the count distribution is shown in Figure 1. The mean counts for the variability sample is 1615. The binning process is carried out using the combine_spectra tool in CIAO to generate the source spectra, background spectra, rmf and arf files in each epoch.

7.1 Method

We use the same spectral fitting models as described in Section 3.1 to perform spectral variability analyses, except for removing the constants $f_{\rm ref}$ and $f_{\rm exs}$ in the model and untying the normalizations of the intrinsic power law and other spectral components. We simultaneously fit the spectra in each epoch using the best-fit model for each source obtained in Section 4, but this time we allow $norm_{\rm ref}$ to vary freely to search for possible variation in the Compton-scattered continuum. Considering the relatively low counts in each epoch, the fitting strategies are adopted as follows:

The uncertainties on $\Gamma$ can cause large degeneracies if we set it free in all epochs. Given that some sources may have large variability in $\Gamma$, we first let $\Gamma$ in each epoch vary freely and obtain $C_{\rm free}$. Then for the first two adjacent epochs (i.e., epoch1 and epoch2), we link $\Gamma$ and measure the $\Delta C$ with respect to $C_{\rm free}$. If $\Delta C$ $\textgreater$ 3.84 (for $\Delta\rm d.o.f=1$), $\Gamma$ is considered to be different between the two epochs at $\textgreater$ 95% confidence level, and we will set $\Gamma$ free; otherwise, we will link $\Gamma$. Then for the second epoch pair (i.e., epoch2 and epoch3), we link their $\Gamma$ and compare the Cstat value with the last step. After traversing all the epoch pairs, if no variability is detected, $\Gamma$ is linked together. 2. 2.

The reflection and soft excess components are often considered to be produced in large-scale clouds, e.g., torus, for highly obscured AGNs. Since the timespan in the rest-frame is (1+$z$) times less than in the observed frame and the typical redshift of our sample is relatively high, it is reasonable to assume that in such a short timescale, the large-scale components may not vary dramatically. To better constrain the normalization of the intrinsic power law, we tie the normalizations of the reflection component and the soft excess component in each epoch, respectively, unless setting them free leading to $\Delta C$ $\textgreater$ 3.84. Finally, only one source shows a significant change in the reflection component, and the soft excess fluxes remain constant for all sources. 3. 3.

For those sources with a weak and constant reflection component confirmed in the last step, we delete the MYTS and MYTL components from the model and only use MYTZ to fit the spectra, in order to reduce the number of free parameters. This procedure may break the self-consistency of the MYTorus model to some extent but will not influence the variability analysis, since we only remove a small constant component in each epoch, which has very little influence on the calculated $\chi^{2}$. 4. 4.

The $N_{\rm H}$ and normalization of the intrinsic power law always vary freely.

The simultaneous fitting yields the best-fit model parameters $\Gamma$, $norm_{\rm ref}$ and $norm_{\rm scat}$. Then we fit the spectra in each epoch with $\Gamma_{i}$, $norm_{\rm ref,i}$ and $norm_{\rm scat,i}$ fixed at the best-fit value obtained in the simultaneous fitting to obtain the $N_{\rm H}$, observed flux, intrinsic flux and the corresponding errors in each epoch. The spectral fitting parameters of each epoch are listed in Table 8. No significant photon index variability is detected for all sources.

7.2 $L_{\rm X}$, $N_{\rm H}$ and Observed Flux Variability

We identify $L_{\rm X}$, $N_{\rm H}$ and observed flux variable sources based on the classical $\chi^{2}$ test. For illustration, we explain how we determine $L_{\rm X}$-variable sources. First, we use the $cflux$ model in XSPEC to calculate the observed and absorption-corrected $\rm 0.5-7$ keV fluxes and corresponding errors. Then we derive the intrinsic rest-frame $\rm 0.5-7$ keV luminosity using equation $L_{\rm 0.5-7\,keV}=4\pi d_{L}^{2}f_{\rm 0.5-7\,keV,int}\,(1+z)^{\Gamma_{\rm int}-2}$, where the “int” subscript represents the intrinsic value. Since the photon index and redshift are all fixed during spectral fitting of each single epoch, the error on $L_{\rm 0.5-7\,keV}$ is totally attributed to the error of the intrinsic flux $f_{\rm 0.5-7\,keV,int}$. Thus the $\chi^{2}$ of the rest-frame $L_{\rm 0.5-7\,keV}$ actually equals the $\chi^{2}$ of the observed-frame $f_{\rm 0.5-7\,keV,int}$. If a source does not vary, i.e., the intrinsic luminosity (or $N_{\rm H}$, flux) is constant, assuming that the errors obey the Gaussian distribution, the $\chi^{2}$ value calculated from

[TABLE]

should obey the $\chi^{2}$ distribution with ($N-1$) d.o.f, where $x$ represents the tested parameter ($L_{\rm X}$, $N_{\rm H}$ or flux), $\sigma_{x,i}$ means the 1$\sigma$ error in the $i$th epoch, and $N$ is the number of epochs. The $\chi^{2}$ test results are listed in Table 9. When the source with 3 (4) binning epochs satisfies $\chi^{2}$ $\textgreater$ 7.8 (9.8), we regard it as being variable at $\textgreater$ 98% confidence level.555The choice of this confidence level is because Y16 showed that it roughly corresponds to $\Delta\rm AIC=4$, so that we can directly compare the $\chi^{2}$ results with those obtained using the AIC method. We also check the $L_{\rm X}$- and $N_{\rm H}$-variable identification results using the Akaike information criterion (AIC) method as described in Y16 that does not need to assume the Gaussian errors (see Section 3.2.3 of Y16 for details). The $\rm\Delta AIC$ values for each source are also listed in Table 9. Sources with \textit{\rm\Delta AIC}>4.0 are assigned as being variable. It can be seen that 30/31 and 30/31 of sources have consistent results in identifying $L_{\rm X}$ and $N_{\rm H}$ variability while using the $\chi^{2}$ test and the AIC method, respectively.666We decide to use the $\chi^{2}$ test results in the following analysis since for the two discrepant sources: CDFS XID 933 and CDFN XID 66, the AIC results show neither $L_{\rm X}$ nor $N_{\rm H}$ variability, which is inconsistent with their flux-variable nature.

Based on the $\chi^{2}$ test results, sources that show observed flux variability account for $55\pm 13$% (17/31) of the entire sample. The resulting $L_{\rm X}$ and $N_{\rm H}$ variable source fractions are $29\pm 10\%$ (9/31) and $19\pm 8\%$ (6/31), respectively. We do not find any correlation between $L_{\rm X}$ and $N_{\rm H}$ variability patterns, suggesting that the main reason that causes the variation of $N_{\rm H}$ is likely not the change of ionization parameter induced by the variable $L_{\rm X}$, but is likely the occultation of the clumpy clouds moving in/out of our l.o.s.

7.3 Variability Amplitude Estimation

We use the normalized excess variance $\sigma^{2}_{\rm nxv}$ (nxv) to estimate the variability amplitude. $\sigma^{2}_{\rm nxv}$  is defined as

[TABLE]

where $N$ is the number of binning epochs (3 or 4), $x_{i}$ and $\delta x_{i}$ are the best-fit parameters (observed flux, $L_{\rm X}$, $N_{\rm H}$) and their 1$\sigma$ errors, $\langle x\rangle$ is the unweighted mean of the calculated parameter. The error on $\sigma^{2}_{\rm nxv}$ is estimated as $s_{\rm D}/(\langle x\rangle^{2}\sqrt{N})$, where

[TABLE]

We calculate $\sigma^{2}_{\rm nxv}$ and corresponding errors on the intrinsic rest-frame 0.5–7 keV luminosity ($\sigma^{2}_{\rm nxv}$,L), observed 0.5–7 keV flux ($\sigma^{2}_{\rm nxv}$,F) and $N_{\rm H}$ ($\sigma^{2}_{\rm nxv}$,N) for each source. The results are listed in Table 9. To check whether our calculation of $\sigma^{2}_{\rm nxv}$ is affected by the limited counts (e.g., Zheng et al. 2017; hereafter Z17), we perform Spearman rank correlation tests and find no significant correlation between $\sigma^{2}_{\rm nxv}$,F, $\sigma^{2}_{\rm nxv}$,N, $\sigma^{2}_{\rm nxv}$,L and counts with Spearman’s $\rho=0.07$, $p$-value = 0.70; $\rho=-0.11$, $p$-value = 0.56 and $\rho=0.32$, $p$-value = 0.08, respectively.

We also calculate the fractional root-mean-square (frms) variability amplitude, which is defined as $\sqrt{\sigma^{2}_{\rm nxv}}$ for sources having $\sigma^{2}_{\rm nxv}$ $\textgreater$ 0 and can be treated as the percentage of the variability amplitude. It is obvious that our sample lacks sources that display large variability amplitudes. The mean frms values of the total sample for the observed flux, intrinsic luminosity and $N_{\rm H}$ are 12%, 10% and 6%, respectively. While for sources having positive nxv, the mean (maximum) frms values are 17% (29%), 19% (34%), and 16% (39%), respectively.

In addition, we do not find any correlation between $\sigma^{2}_{\rm nxv}$,F, $\sigma^{2}_{\rm nxv}$,L and $L_{\rm X}$ (Figure 19) whereas other studies suggested a negative trend, i.e., luminous AGNs show relatively small variations (Y16; Z17; Paolillo et al., 2017). We perform a Spearman correlation test for the 25 common sources using the $\sigma^{2}_{\rm nxv}$,F derived in Z17 which is obtained from the light curve analysis, the result is still consistent with no correlation between $L_{\rm X}$ and $\sigma^{2}_{\rm nxv}$,F (Spearman’s $\rho=-0.31$, $p$-value = 0.13). As pointed out by Allevato et al. (2013), the calculation of $\sigma^{2}_{\rm nxv}$ is biased at low counts and uneven cadence and should be restricted to large samples. Therefore, the limited source number and the broad redshift span may prevent us from detecting such a correlation. But beyond that, since Y16 and Z17 mainly focused on the AGNs with largest counts available in the CDF-S rather than highly obscured AGNs, the variability behavior and its underlying physical drivers of our sources and those in Y16 and Z17 may be intrinsically different.

Recently, González-Martín (2018) found a new X-ray variability plane for AGNs that links the characteristic break frequency of the PSD ($f_{\rm break}$) with bolometric luminosity and $N_{\rm H}$. This makes $N_{\rm H}$ play a non-negligible role in understanding AGN variability. As we will show in the next section, although Z17 argued that the $N_{\rm H}$ variation does not contribute significantly to the total variability in their sample, the $N_{\rm H}$ variability does play a substantial role for highly obscured AGNs. Nevertheless, the small $\sigma^{2}_{\rm nxv}$ values indicate that the state of highly obscured AGNs does not vary dramatically while concerning the average source properties on a 17-year timescale in the observed-frame.

7.4 Detailed Variability Analysis

There are some sources in our sample that show significant variability of the observed flux, intrinsic luminosity or $N_{\rm H}$. In this section, we perform detailed analyses to study their variability behavior aiming at shedding light on the leading mechanism that drives the large variabilities, and try to understand the typical location and size of the obscuring clouds (all the following sources are in the CDF-S) .

XID 63

This source has 848 counts available and its data are binned to three epochs. It is a moderately luminous source ($\bar{L}\rm_{X}=\rm 5.1\times 10^{43}\ erg\ s^{-1}$) at $z=0.737$. The best-fit model of this source is the absorbed power-law ($\Gamma=1.54^{+0.21}_{-0.22}$) plus an additional soft excess component. $\Gamma$ and $norm_{\rm soft}$ do not show variations and are linked during all the epochs. The unfolded spectra, light curves, as well as the 1$\sigma$ and 2$\sigma$ confidence contours of $N_{\rm H}$ and normalization are displayed in the left column of Figure 20. It has $\chi_{N_{\rm H}}^{2}=65.2$, $\chi_{f,\rm obs}^{2}=34.2$ and $\chi_{L}^{2}=7.2$, indicating significant variabilities. Although the photon index does not vary, the observed spectral shape changes prominently due to the large variation in column density. The absorption is weak in epoch1 with $N_{\rm H}$ $\rm=5.8\times 10^{22}\ cm^{-2}$. Then it transforms into a highly obscured state with $N_{\rm H}$ $=\rm 1.4\times 10^{23}\ cm^{-2}$ in epoch2, and its $N_{\rm H}$ continues to rise to $\rm 1.9\times 10^{23}\ cm^{-2}$ in epoch3. The intrinsic flux increases about a factor of 1.4 from epoch1 to epoch2, and decreases about a factor of 1.6 from epoch2 to epoch3. Due to the large variability amplitude in $N_{\rm H}$, the observed flux remains constant in the first two epochs, but significantly declines in epoch3 by a factor of 1.9. The $N_{\rm H}$ variability may be caused by the clumpy cloud moving into the l.o.s. Since the X-ray obscuration type transition happens between epoch1 and epoch2, which corresponds to a rest-frame timescale $t<6.0$ years ($\approx$ 10.4 years in the observed frame). By assuming this transition as an eclipse event, we can roughly estimate the distance and angular size of the cloud using the same method in Y16, thus to constrain the cloud size. We use the empirical relation between the inner torus radius $r$ and 14 - 195 keV luminosity from Koshida et al. (2014) to estimate the distance. The relation can be described as ${\rm log}\,r=-1.04+0.5\,{\rm log}\,(L_{\rm 14-195\,keV}/10^{44})$, where $r$ and $L_{\rm 14-195keV}$ are in units of pc and $\rm erg\ s^{-1}$, respectively. The $L\rm_{14-195\,keV}=2.5\times 10^{44}\ erg\ s^{-1}$ at epoch1 is obtained by extrapolating the power-law to the higher energy band. Thus the estimated $r$ from the inner torus to the central black hole is $\approx 0.14$ pc. By assuming Keplerian motion, the orbital period of the cloud $t_{\rm orbit}=2\pi\frac{r^{\frac{3}{2}}}{(M_{\rm BH}G)^{\frac{1}{2}}}\approx 519\,(\frac{M_{\rm BH}}{10^{8}M_{\odot}})^{-\frac{1}{2}}$ year. The angular size (viewed from the BH) of the cloud is estimated as $\frac{t_{\rm tran}}{t_{\rm orbit}}\times 360^{\circ}<4.1^{\circ}(\frac{M_{\rm BH}}{10^{8}M_{\odot}})^{\frac{1}{2}}$. By multiplying $r$, we estimate the cloud size to be $<0.01\,(\frac{M_{\rm BH}}{10^{8}M_{\odot}})^{\frac{1}{2}}\rm\,pc$.

XID 249

This source can be well fitted by the absorbed power-law model with $\Gamma$ = $1.51^{+0.26}_{-0.28}$ and has the most prominent $N_{\rm H}$ variability in our sample ($\chi^{2}_{N\rm_{H}}$ = 130.8) as shown in Figure 20. The $N_{\rm H}$ declines from $\rm 5.0\times 10^{23}\ cm^{-2}$ to $\rm 2.8\times 10^{23}\ cm^{-2}$, and finally $\rm 1.5\times 10^{23}\ cm^{-2}$ during the three epochs. The intrinsic flux of this source remains roughly constant ($\chi_{flux,\rm in}^{2}=2.2$), but the strong variation in $N_{\rm H}$ causes a significant increase in the observed flux ($\chi_{flux,\rm obs}^{2}=31.0$). The variability analyses of each single epoch spectrum provide very different results compared with the $\rm 7\ Ms$ stacked spectrum ($\Gamma$ = 1.20 when using MYTZ alone which does not limit the range of $\Gamma$), also indicating strong $N_{\rm H}$ variation. Using the same method as for XID 63, we estimate the distance and the cloud size to be $<0.09\rm\ pc$ and $<0.006\,(\frac{M_{\rm BH}}{10^{8}M_{\odot}})^{\frac{1}{2}}\rm\,pc$ for the occultation event from epoch2 to epoch3.

XID 328

This source has the highest column density in our variability sample. The epoch-mean $L_{\rm X}$ and $N_{\rm H}$ are $\rm 3.7\times 10^{44}\ erg\ s^{-1}$ and $\rm 8.5\times 10^{23}\ cm^{-2}$, respectively, which are consistent with the stacked spectral analyses. The spectra can be well described by a power law ($\Gamma$ = $2.12^{+0.24}_{-0.24}$) plus Compton reflection continuum and strong Fe K lines. $N_{\rm H}$ and the reflection flux remain roughly constant in the timespan of 5.8 years in the rest-frame, and the observed flux variations result from the intrinsic X-ray power variability. According to the stacked spectral analyses, the reflection flux accounts for 27.5% of the observed 0.5–7 keV flux, but only accounts for 3.8% of the intrinsic 0.5–7 keV flux. This source has also been reported in the literature as CT candidates (Comastri et al. 2011; Y16). Our results confirm its highly obscured nature but with higher $N_{\rm H}$, $L_{\rm X}$ and $\Gamma$.

XID 49

This source was observed 98 times during the total $\rm 7\ Ms$ campaign and has 2345 counts available. The spectra can be explained by a simple absorbed power law ($\Gamma=1.66^{+0.11}_{-0.12}$) with different column densities and normalizations (see Figure 20). The intrinsic flux decreases from $\rm 1.57\times 10^{-14}\ erg\ cm^{-2}\ s^{-1}$ to $\rm 7.45\times 10^{-15}\ erg\ cm^{-2}\ s^{-1}$ during the first two epochs and finally rises again to $\rm 1.39\times 10^{-14}\ erg\ cm^{-2}\ s^{-1}$ in the fourth epoch. $N_{\rm H}$ has a wave-shape variability behavior and does not follow the intrinsic flux variability pattern, indicating that the varying absorption is not caused by the change of ionization state.

XID 458

This source has a redshift of 2.291 and can be described by a power law (\Gamma$$=1.83^{+0.14}_{-0.13}) with reflection hump and iron emission lines. Though the $N_{\rm H}$ is not so high during all the epochs ($\bar{N}\rm_{H}\approx\rm 1.4\times 10^{23}\ cm^{-2}$), the reflection component is required to fit the $\rm 7\ Ms$ stacked spectrum. The intrinsic flux increases from epoch1 to epoch3, and declines at epoch4. But the reflection flux remains roughly constant from epoch1 to epoch2, decreases at epoch3, and does not show any variability from epoch3 to epoch4. The different variability patterns indicate that there may be a time lag between the intrinsic continuum variability and the large-scale reflection flux variability. The decline of the reflection flux during epoch2 to epoch3 (from observed-frame 2007 September to 2010 July, corresponding to rest-frame 0.76 years) might result from significant flux decline before epoch1, and this decline has just propagated to the reflection medium (possibly the clumpy torus) and causes significant diminishment of the reflection flux. This means that the variable continuum signal needs at least 2.4 years (rest-frame timespan from 1999 October to 2007 September) to spread to the torus from the central emission region, which provides a rough lower-limit estimate, $\approx$ 0.7 pc, of the location of the reflecting cloud.

Remaining sources

By applying similar analyses to the remaining sources, we find that 17 sources in our sample that show observed flux variabilities can be classified into three types: 29% (5/17) are caused by the change of $N_{\rm H}$, 53% (9/17) are caused by the AGN intrinsic luminosity variability, and 6% (1/17) are driven by both effects. Note that there are two sources (CDF-S XIDs 186 and 876) identified to be flux-variable but showing neither $L_{\rm X}$- nor $N_{\rm H}$-variability according to the $\chi^{2}$ tests. However, their best-fit intrinsic luminosities both show obvious variations, but due to large errors, the corresponding $\chi^{2}_{L}$ values are smaller than the critical value. The high $N_{\rm H}$-variable fraction among flux-variable sources confirms our previous thoughts in Section 7.3 that the $N_{\rm H}$ variability is a key ingredient for investigating the variability in highly obscured AGNs.

8 Conclusions

In this study, we present systematic X-ray spectral analyses of 436 highly obscured AGNs ($N_{\rm H}$ $\textgreater$ $\rm 10^{23}\ cm^{-2}$) with $L_{\rm X}$ $\textgreater$ $10^{42}$ $\rm erg\ s^{-1}$ identified in the Chandra Deep Fields, which make up the largest dedicated highly obscured AGN sample to date, to explore their physical properties and evolution. We also carry out detailed long-term variability analyses for a subsample of 31 sources with largest counts available to investigate the main driver of their spectral variability and the typical variability amplitude. Below we summarize our main results.

We perform detailed X-ray spectral fitting for 1152 AGNs in the Chandra Deep Fields with observed-frame 0.5–7 keV net counts $\textgreater$ 20 using physically appropriate MYTorus-based models, in order to identify heavily obscured ones. By limiting our analyses to sources with \textit{{\rm log},L_{\rm X}}>42\ \textit{\rm erg\ s^{-1}} in order to avoid possible contamination from star-forming galaxies, 436 AGNs are confirmed to be highly obscured, with $\overline{z}=1.9$ and \overline{\textit{{\rm log},L_{\rm X}}}=\rm 43.6\ \textit{\rm erg\ s^{-1}}. 2. 2.

We select 102 Compton-thick (CT) candidates with best-fit $N_{\rm H}$ $\textgreater$ $\rm 10^{24}\ cm^{-2}$ and 1$\sigma$ lower limit $\textgreater$ $\rm 5\times 10^{23}\ cm^{-2}$, accounting for $\sim$ 23% of the highly obscured sample. The observed log$\,N$-log$\,S$ for CT AGNs prefers the moderate CT number counts as predicted by Akylas et al. (2012) and Ueda et al. (2014) cosmic X-ray background models, while other models (e.g., Gilli et al., 2007; Ballantyne et al., 2011) more or less overestimate or underestimate the number counts. 3. 3.

We present a new hardness-ratio measure of the obscuration level as a function of redshift (the HR curve), which can be used to select heavily obscured AGN candidates without resorting to detailed spectral fitting. The completeness and accuracy by applying the HR curve on the CDFs AGN population to identify highly obscured ones are 88% and 80%, respectively. 4. 4.

We find a strong negative correlation between the soft excess fraction $f_{\rm exs}$ and $N_{\rm H}$ with a Spearman rank correlation coefficient $\rho=-0.66$. By assuming that the soft excess originates from the scattered-back continuum and treating the small $f_{\rm ref}$ as an indicator of high covering factor of the obscuring materials, this result indicates that a portion of the most heavily absorbed AGNs reside in an extremely geometrically buried environment. 5. 5.

Among the 31 CDF-S highly obscured sources which have optical classification results, 19% (6/31) of them are labeled as BLAGN, indicating that at least for part of the AGN population, the high-level X-ray obscuration is largely a line-of-sight effect, i.e., some high-column-density clouds on various scales (but not necessarily a dust-enshrouded torus) may obscure the compact X-ray emitter without blocking the entire BLR; alternatively, the heavy X-ray obscuration maybe produced by the BLR itself. 6. 6.

After considering the errors on the best-fit $N_{\rm H}$ and correcting for the sky coverage effect, the Eddington bias and the $N_{\rm H}$-dependent Malmquist bias, we derive the intrinsic $N_{\rm H}$ distribution representative of the highly obscured AGN population as well as its evolution across cosmic time. The intrinsic CT-to-highly-obscured fraction is roughly 52% and is consistent with no evident redshift evolution. 7. 7.

We select 31 sources with 0.5–7 keV net counts $\textgreater$ 700 and $\textgreater$ 900 in the CDF-S and CDF-N, respectively, to perform long-term ($\approx 17$ years in the observed frame) spectral variability analyses. We find that the flux-variable, $L_{\rm X}$-variable and $N_{\rm H}$-variable source fractions are $55\pm 13\%$ (17/31), $29\pm 10\%$ (9/31) and $19\pm 8\%$ (6/31), respectively. The typical flux, $L_{\rm X}$ and $N_{\rm H}$ variability amplitudes for those variable sources are 17%, 19%, and 16%, respectively. 8. 8.

We calculate the normalized excess variance ($\sigma_{\rm nxv}^{2}$) of $N_{\rm H}$, observed 0.5–7 keV flux and $L_{\rm X}$ to investigate the intrinsic variability amplitude. No correlation between $\sigma_{\rm nxv}^{2}$ and $L_{\rm X}$ is detected, possibly due to the observed flux variability for a significant fraction of highly obscured AGNs being caused by the change of $N_{\rm H}$ rather than the variation of $L_{\rm X}$ alone, as well as the limited sample size and the broad redshift span. 9. 9.

We discuss detailed variability behaviors for 5 sources that show significant $N_{\rm H}$, $L_{\rm X}$ or observed flux variability. The typical locations and sizes of the occultation and reflection clouds are estimated (see Section 7.4). The main driver for the variability of the 17 flux-variable sources can be classified into three types: 29% (5/17) are caused by the change of $N_{\rm H}$, 53% (9/17) are caused by the AGN luminosity variability and 6% (1/17) are driven by both effects. Two sources are not classified due to large measurement errors.

Benefiting from the deepest X-ray surveys to date, our work provides meaningful constraints on the properties and evolution of the AGN obscuring materials over a broad redshift range, and quantifies the detailed variability behaviors of these hidden sources, which are crucial for us to better understand the role that highly obscured AGNs played in galaxy evolution. However, X-ray data alone are incapable to show the overall perspective. In a subsequent paper of this series (Li et al., in preparation), we will further explore the properties of highly obscured AGNs and their host galaxies by combining 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 comments and suggestions. We thank David Ballantyne, Murray Brightman, Giorgio Lanzuisi and Yoshihiro Ueda for providing relevant data, and Tahir Yaqoob and Mislav Baloković for helpful discussions. J.Y.L., Y.Q.X., M.Y.S., and X.C.Z. acknowledge support from the 973 Program (2015CB857004), the National Natural Science Foundation of China (NSFC-11890693, 11473026, 11421303), the CAS Frontier Science Key Research Program (QYZDJ-SSW-SLH006), and the K.C. Wong Education Foundation. M.Y.S and X.C.Z. acknowledge the support from the National Natural Science Foundation of China (NSFC-11603022) and the China Postdoctoral Science Foundation (2016M600485). X.C.Z also acknowledges the financial support from CSC(China Scholarship Council)-Leiden University joint scholarship program. FV acknowledges support from CONICYT and CASSACA through the Fourth call for tenders of the CAS-CONICYT Fund. T.M.H. acknowledges the support from the Chinese Academy of Sciences (CAS) and the National Commission for Scientific and Technological Research of Chile (CONICYT) through a CAS-CONICYT Joint Postdoctoral Fellowship administered by the CAS South America Center for Astronomy (CASSACA) in Santiago, Chile. B.L. acknowledges financial support from the National Key R&D Program of China grant 2016YFA0400702 and National Natural Science Foundation of China grant 11673010. CV acknowledges financial support from the Italian Space Agency (ASI) under the contracts ASI-INAF I/037/12/0 and ASI-INAF n.2017-14-H.0. RG acknowledges support from the agreement ASI-INAF n.2017-14-H.O. X.S. acknowledges support from the NSFC-11573001 and 11822301.

Appendix A Justification of the adopted models and parameters

In the MYTorus configuration, since the reproccessed components strongly depend on the inclination angle $\theta$, the equatorial $N_{\rm H}$ and the assumed $f_{\rm 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 parameters and simply fix them to values that their representativeness has not yet been physically validated. Here we justify that our spectral fitting results are not significantly influenced by the assumed parameters.

We first fit the spectra for the high count sources (counts $\textgreater$ 200) with $\theta$ and $f_{\rm ref}$ set as free parameters. The best-fitted $\theta$ peaks at $\theta<70^{\circ}$, therefore this time we choose to fix it at $65^{\circ}$ for all sources. Then we re-fit the spectra by keeping $f_{\rm ref}$ free and obtain the best-fit results. For the low count sources, we fix $f_{\rm ref}$ at 2.0 for sources with detected reprocessed components, which is roughly the average value for high count sources.

The comparison of the results with those obtained in Section 3.1 is shown in Figure 21. Despite of small scatters, the results are consistent. Therefore our simple choices of the fixed parameters are reasonable. The consistency between the MYTorus results and the Borus results obtained in Section 5.3 also validates our spectral fitting strategy.