Survival of the Obscuring Torus in the Most Powerful Active Galactic Nuclei

Dedicated searches for active galactic nuclei (AGN) generally find that the fraction of AGN classified either as optical type-2 (obscured) or X-ray absorbed decreases substantially with increasing luminosity (Lawrence & Elvis 1982; Hasinger et al. 2005; Simpson 2005; Della Ceca et al. 2008; Burlon et al. 2011; Merloni et al. 2014; Ueda et al. 2014; Buchner et al. 2015).

To explain these findings, receding torus models have often been adopted (Lawrence 1991; Simpson 2005; Hönig & Beckert 2007). They postulate that the geometry of the material obscuring the AGN nuclear region, the dusty torus (Antonucci 1993; García-Burillo et al. 2016), changes with AGN luminosity. The torus geometrical covering factor (hereafter f₂) defines the fraction of the sky around the AGN central engine that is obscured. If f₂ decreases with increasing AGN luminosity, then this considerably reduces the probability of finding luminous type-2 AGN (Elitzur 2012).

The observed decrease of the ratio of the torus infrared luminosity and the AGN bolometric luminosity (${L}_{\mathrm{torus}}/{L}_{\mathrm{bol}}$) with ${L}_{\mathrm{bol}}$ has often been interpreted as direct evidence of a receding torus (Maiolino et al. 2007; Treister et al. 2008; Lusso et al. 2013; but see Netzer et al. 2016). These results should be treated with caution since a one-to-one correspondence is not expected between ${L}_{\mathrm{torus}}/{L}_{\mathrm{bol}}$ and f₂ (Stalevski et al. 2016).

We can determine f₂ using radiative transfer models that self-consistently reproduce the emission from dust in the torus heated by the AGN (e.g., Fritz et al. 2006; Hönig & Kishimoto 2010; Nenkova et al. 2008; Schartmann et al. 2008; Stalevski et al. 2016). Using torus models with a clumpy distribution of dust from Nenkova et al. (2008, hereafter N08), we determined, for the first time, the probability density distributions of f₂ for individual objects (Mateos et al. 2016, hereafter M16) for a large, uniformly selected, and complete flux-limited sample of X-ray selected AGN drawn from the Bright Ultra-hard XMM-Newton Survey (BUXS; Mateos et al. 2012, 2013).

Using the distributions of f₂ for the AGN in BUXS, we derive here the intrinsic fraction of optical type-2 AGN at redshifts $0.05\leqslant z\leqslant 1$ as a function of intrinsic (absorption-corrected) rest-frame 2–10 keV luminosity from ${10}^{42}$ to ${10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ (hereafter ${L}_{{\rm{X}}}$). We also investigate whether the decrease of f₂ with ${L}_{{\rm{X}}}$, which we observe in the BUXS sample, is a property of the AGN population. Throughout, errors are 1σ (the 16th and 84th percentiles when referring to distributions) unless otherwise stated. We adopt the concordance cosmology, ${{\rm{\Omega }}}_{{\rm{M}}}=0.3$, ${{\rm{\Omega }}}_{\lambda }=0.7$, and H₀ = 70$\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$.

Our AGN sample is drawn from the BUXS survey. BUXS includes 255 X-ray bright (${f}_{4.5-10\mathrm{keV}}\gt 6\,\times {10}^{-14}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$) AGN detected with XMM-Newton in the 4.5–10 keV band over 44.43 deg². Out of these, 252 have robust redshift z and optical spectroscopic classifications. Objects with detected rest-frame UV/optical broad emission lines (FWHM $\geqslant \,1500$ km ${{\rm{s}}}^{-1}$) are classified as type-1 and those with narrower emission lines as type-2.

Here, we only consider AGN with ${L}_{{\rm{X}}}\geqslant {10}^{42}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ to minimize host galaxy contamination, by increasing the AGN to galaxy contrast ratio, and with $z\leqslant 1$ to avoid strong evolutionary effects. This restricts our sample to 199 objects, with ${10}^{42}\,\leqslant {L}_{{\rm{X}}}\leqslant {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ and $0.05\leqslant z\leqslant 1$.

BUXS is a unique survey to conduct this study. It is the only AGN sample for which we know f₂ for almost all (∼99%) objects. It is sufficiently large to accurately constrain the intrinsic type-2 AGN fraction. For all sources we have good-quality X-ray spectroscopy (∼few hundred counts) that guarantees robust estimates of ${L}_{{\rm{X}}}$. Assuming the worst case, that all three unidentified sources are in the $z,{L}_{{\rm{X}}}$ interval under study, they represent at most ∼1.5% of our sample. Clearly, our results are not affected by identification incompleteness effects that would bias against obscured AGN.

To compute the luminosity dependence of the type-2 AGN fraction, we further divided our sample into three luminosity bins of equal logarithmic width: ${10}^{42}\mbox{--}{10}^{43},{10}^{43}\mbox{--}{10}^{44}$, and ${10}^{44}\mbox{--}{10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ (see Table 1).

Table 1.  Summary of the Properties of Our AGN Samples and the Results of Our Analysis

log(${L}_{{\rm{X}}}$)	n1	n2	$\langle {L}_{{\rm{X}}}{\rangle }_{1}$	$\langle z{\rangle }_{1}$	$\langle {L}_{{\rm{X}}}{\rangle }_{2}$	$\langle z{\rangle }_{2}$	Observed	N2	Intrinsic	Type-2 Fraction	Compton-thick
							Type-2 Fraction		Type-2 Fraction	Missed	Fraction
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
42–43	16	21	42.75	0.10	42.80	0.11	${56.8}_{-8.4}^{+7.5}$	${5}_{-4}^{+6}$	${64.4}_{-9.3}^{+7.1}$	${25.3}_{-14.1}^{+16.7}$	≤${14.5}_{-7.7}^{+14.6}$
43–44	53	33	43.59	0.29	43.48	0.25	${38.4}_{-5.0}^{+5.4}$	${38}_{-9}^{+18}$	${58.7}_{-6.8}^{+5.6}$	${56.9}_{-12.2}^{+7.5}$	≤${32.6}_{-8.3}^{+8.1}$
44–45	55	21	44.50	0.76	44.41	0.62	${27.6}_{-4.5}^{+5.6}$	${40}_{-14}^{+19}$	${54.6}_{-8.7}^{+6.8}$	${69.7}_{-14.4}^{+6.1}$	≤${37.0}_{-10.5}^{+8.9}$

Note. Column 1: X-ray luminosity range in $\mathrm{erg}\,{{\rm{s}}}^{-1}$ in logarithmic units; columns 2 and 3: number of type-1 and type-2 AGN in the bin, respectively; columns 4 to 7: median X-ray luminosity and redshift of type-1 and type-2 AGN, respectively; column 8: observed type-2 AGN fraction; column 9: number of type-2 AGN missed; column 10: intrinsic type-2 AGN fraction; column 11: fraction of type-2 AGN that have escaped X-ray detection; column 12: fraction of type-2 AGN that have escaped X-ray detection over the total population. Fractions are given in percentage units.

Download table as:  ASCIITypeset image

In Mateos et al. (2015), we built the rest-frame UV-to-infrared spectral energy distributions (SEDs) of our objects using data from the Sloan Digital Sky Survey (Abazajian et al. 2009), the Two Micron All Sky Survey (Cutri et al. 2003), the UKIRT Infrared Deep Sky Survey (Lawrence et al. 2007), and the Wide Field Infrared Survey Explorer (WISE; Wright et al. 2010). With an SED decomposition analysis we isolated the emission associated with dust in the torus heated by the AGN at rest-frame wavelengths from ∼1 to ∼22 μm.

In M16, we fitted the torus SEDs with the N08 models using the Bayesian inference tool BayesCLUMPY that provides posterior distributions for all the free parameters of the models (Asensio Ramos & Ramos Almeida 2009). We used truncated uniform prior distributions for all the torus model parameters (six in total) in the ranges listed in Table 1 from M16.

In the N08 models, the dust is distributed in optically thick clouds (${\tau }_{V}\gt 1$ at 5500 Å). The torus inner radius is set by the sublimation temperature of the dust grains ($\approx 1500$ K). The radial distribution of clouds declines as a power law. The vertical angular distribution of clouds has no sharp boundary and it is parameterized with a Gaussian. The geometrical covering factor of the torus f₂ is defined as

$$\begin{eqnarray}&&{f}_{2}=1-{\int }_{0}^{\pi /2}{P}_{\mathrm{esc}}(\beta )\cos (\beta ){\text{}}d\beta ,\end{eqnarray} \tag{ 1 }$$

where ${P}_{\mathrm{esc}}$ is the probability that light from the AGN will escape without being absorbed at an angle β from the torus equatorial plane:

$$\begin{eqnarray}&&{P}_{\mathrm{esc}}(\beta )={e}^{-{N}_{0}\times {e}^{(-{\beta }^{2}/{\sigma }^{2})}}.\end{eqnarray} \tag{ 2 }$$

f₂ depends on the angular width of the torus (σ) and the mean number of clouds along the equatorial direction (N₀). Using the posterior distributions of σ and N₀ we calculated the probability density distribution of f₂ for each source.

Since f₂ represents a geometrical covering factor, in any AGN subpopulation having a given dust covering factor f₂, the fraction of type-2 objects should intrinsically be f₂ (Elitzur 2012).

We started by computing the observed type-2 AGN fraction in BUXS as a function of f₂, fully taking into account the uncertainties in f₂. We divided the range of f₂ into five bins of width Δf₂ = 0.2, the 1σ average error in our f₂ estimates from the individual distributions. For each source, we obtained the fraction of f₂ in each bin by integrating its f₂ probability distribution. The observed type-2 AGN fraction in bin i, F_obsⁱ, is defined as

$$\begin{eqnarray}&&{F}_{\mathrm{obs}}^{i}=\displaystyle \frac{{\sum }_{j=1}^{{n}_{2}}{F}_{2,j}^{i}}{{\sum }_{j=1}^{{n}_{2}}{F}_{2,j}^{i}+{\sum }_{k=1}^{{n}_{1}}{F}_{1,k}^{i}},\end{eqnarray} \tag{ 3 }$$

where n₁ and n₂ are the number of type-1 and type-2 AGN, respectively (see Table 1) and ${F}_{1,k}^{i}$ and ${F}_{2,j}^{i}$ are the fractions of the probability distributions of f₂ in bin i.

To compute the uncertainties in F_obsⁱ, fully taking into account both source Poisson counting noise and the uncertainties in f₂ for each source (the full f₂ distributions), we used a bootstrap error estimate. We generated ${10}^{6}$ mock samples by randomly selecting type-1 and type-2 objects, with their corresponding f₂ distributions, from the original samples with replacement. Each mock sample contained a number of type-1 and type-2 AGN that was calculated from binomial distributions keeping constant the observed number of sources, i.e., (n₁ + n₂), and assuming that the true type-2 AGN fraction is the observed one, n₂/(n₁ + n₂). We computed F_obsⁱ for each simulated data set and then determined the median and percentiles.

Since we will use a chi-squared (${\chi }^{2}$) fit to derive the intrinsic type-2 AGN fraction in the next section, we have corrected for the small asymmetry in the F_obsⁱ errors. For each set of values (${F}_{\mathrm{obs}},{F}_{\mathrm{obs}}^{-},{F}_{\mathrm{obs}}^{+}$) we computed the Gaussian function that has mean F_obs and dispersion ${\rm{\Delta }}{F}_{\mathrm{obs}}$ such that its integral from ${F}_{\mathrm{obs}}-{F}_{\mathrm{obs}}^{-}$ to F_obs + ${F}_{\mathrm{obs}}^{+}$ is 68.3%. We used these ${\rm{\Delta }}{F}_{\mathrm{obs}}$ as σ errors in our ${\chi }^{2}$ fits.

Our results are illustrated in Figure 1. The y-axis represents the observed fraction of AGN classified as optical type-2 in BUXS. The x-axis represents the covering factor of the torus inferred from SED-fitting with N08 models. If BUXS did not miss any highly covered AGN (f₂ ∼ 1), our points should follow the 1:1 relation between the type-2 AGN fraction and f₂. Clearly, this is not the case, especially at ${L}_{{\rm{X}}}\gt {10}^{43}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. There are not enough luminous type-2 AGN with high f₂; therefore, some must have escaped X-ray detection.

Zoom In Zoom Out Reset image size

We have derived the global intrinsic type-2 AGN fraction by requiring that the type-2 AGN fraction and f₂ follow a 1:1 relation, i.e., for each f₂ the intrinsic fraction of type-2 AGN must be equal to f₂. To do so we made the following assumptions.

(i) The AGN missed are all type-2. BUXS is not biased against type-1 AGN.¹⁰ Moreover, the AGN missed cannot be "X-ray weak" AGN, not only because these are a rare population (e.g., Brandt et al. 2000; Risaliti et al. 2003), but also because there is no physical reason why such objects should have tori with the highest f₂ among all AGN. Flux-limited surveys below 10 keV are incomplete for low-z AGN whose line of sight X-ray absorption approaches the Compton-thick limit (${N}_{{\rm{H}}}=1.5\times {10}^{24}\,{\mathrm{cm}}^{-2}$). Thus, the type-2 AGN missed are either Compton-thick or heavily absorbed (${N}_{{\rm{H}}}\sim {\rm{a}}\,\mathrm{few}\,\times {10}^{23}\,{\mathrm{cm}}^{-2}$).

(ii) The stacked probability distribution of f₂ for all the type-2 AGN in BUXS with ${N}_{{\rm{H}}}\gt 4\times {10}^{23}\,{\mathrm{cm}}^{-2}$ (insets in Figure 1) represents well f₂ in the objects missed. This is well justified since heavily absorbed type-2 AGN have f₂ ∼ 1 (Ramos Almeida et al. 2011; M16) and their distributions of f₂ are all very similar, regardless of ${L}_{{\rm{X}}}$ (Alonso-Herrero et al. 2011; Ichikawa et al. 2015; M16). If we assume instead a distribution of f₂ peaking at smaller values, our main results remain unchanged, although only poorer fits are possible.

To compute the number of type-2 AGN missed, N₂, we increased N₂ until we found the best ${\chi }^{2}$ fit to the 1:1 relations in Figure 1. The best-fit type-2 AGN fractions, F_bfⁱ (open squares), are

$$\begin{eqnarray}&&{F}_{{bf}}^{i}=\displaystyle \frac{{\sum }_{j=1}^{{n}_{2}}{F}_{2,j}^{i}+{F}_{\mathrm{CT}}^{i}\times {N}_{2}}{{\sum }_{j=1}^{{n}_{2}}{F}_{2,j}^{i}+{\sum }_{k=1}^{{n}_{1}}{F}_{1,k}^{i}+{F}_{\mathrm{CT}}^{i}\times {N}_{2}},\end{eqnarray} \tag{ 4 }$$

where F_CTⁱ is the fraction of the f₂ distribution used to represent the missing AGN in the bin i. We obtained the uncertainties on N₂ from our fits as $P({N}_{2})\propto {e}^{-{\rm{\Delta }}{\chi }^{2}/2}$. Next, we computed the intrinsic type-2 AGN fraction, F_intr = (n₂ + N₂)/(n₁ + n₂ + N₂), with uncertainties using a Bayesian approach (Wall & Jenkins 2003). We have assumed a binomial distribution similar to that in Section 4 weighted with $P({N}_{2})$. Since $P({N}_{2})$ are only defined for ${N}_{2}\geqslant 0$, our estimates of F_intr are higher than what would be directly obtained from the numbers listed in columns 2, 3, and 9 of Table 1. We have followed this same Bayesian approach to compute the fractions listed in columns 8, 11, and 12 in Table 1.

Since the errors of F_obsⁱ are not independent, we tested the robustness of our ${\chi }^{2}$ fitting with Monte Carlo simulations. We used our simulated values of F_obsⁱ (see Section 4) to determine N₂ with a least squares fitting and found that the results were indistinguishable from those obtained with the ${\chi }^{2}$ fitting.

We now discuss some issues that could affect our results.

In our computations of F_obsⁱ, objects with f₂ ∼ 0.5 contribute to both the upper and lower f₂ neighboring bins, while objects with extreme f₂ values contribute only toward the central bins. This effect could flatten our F_obsⁱ estimates in Figure 1 with respect to the 1:1 relation removing the need of adding any missing sources. To address this issue, we have used Monte Carlo simulations of source samples of the same size as ours but following the 1:1 relation. The simulated objects are assigned f₂ probability distributions drawn from the real type-1/type-2 objects, proportionally to their weight at the needed f₂ value. This results only in a difference of a few percent in our F_intr estimates, not changing our conclusions.

For two type-1 AGN with ${L}_{{\rm{X}}}\gt {10}^{43}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ we could not determine f₂ because they were not detected by WISE. We included both sources in our analysis by assigning them f₂ distributions drawn at random with replacement from the type-1 AGN sample in the corresponding luminosity bin. Since all type-1 AGN have similar distributions of f₂ (M16), our approach is well justified. Therefore, our results are free from systematic uncertainties associated with the lack of detection of X-ray sources at infrared wavelengths.

Type-2 AGN typically have optical extinctions associated with dust in their hosts of ${A}_{V}\leqslant 5$ mag (e.g., Alonso-Herrero et al. 2003, 2011). These are too small to have any noticeable effect on our torus SEDs and thus on our f₂ estimates. Host galaxy dilution may cause a type-1 AGN to be misclassified as type-2, especially at low luminosities (Caccianiga et al. 2007). If there are cases such as these in BUXS, the number is too small to affect our results (M16).

By restricting our analysis to AGN with ${L}_{{\rm{X}}},z$ where BUXS is complete, we found that the uncertainties were larger, but the results remained unchanged. Finally, we have verified that if we force type-1 and type-2 AGN to have the same z distributions on each ${L}_{{\rm{X}}}$ bin (with a bootstrap re-sampling), the results remain the same. Thus, if type-1 and type-2 AGN evolve differently at the z under study (e.g., Reyes et al. 2008), this has no noticeable impact on our analysis.

A non-negligible fraction of luminous type-2 AGN with high f₂ have escaped X-ray detection. The fraction of AGN missed over the total population (last column in Table 1) gives us a strict upper limit to the Compton-thick fraction, since none of the AGN in BUXS are Compton-thick. At ${L}_{{\rm{X}}}\gt {10}^{43}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, Compton-thick AGN cannot contribute more than ${37.0}_{-10.5}^{+8.9}$% to the total AGN population, in agreement with recent estimates (Buchner et al. 2015; Ricci et al. 2015). At ${L}_{{\rm{X}}}\lt {10}^{43}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ ($z\sim 0.1$), our results are also consistent, within the uncertainties, with Compton-thick fractions reported for low-z AGN samples (Bassani et al. 2006; Burlon et al. 2011; Ricci et al. 2015; Akylas et al. 2016).

Figure 2 shows the luminosity dependence of F_intr. For comparison, we also show the observed type-2 AGN fraction in BUXS which, as typically found in flux-limited X-ray surveys, decreases substantially with ${L}_{{\rm{X}}}$. We have parameterized the luminosity dependence of F_intr with a power law of the form F_intr ∝ ${L}_{{\rm{X}}}^{-\alpha }$. A simple ${\chi }^{2}$ fit to a straight line in log–log space using the values and errors of F_intr in Table 1 yields the relation

$$\begin{eqnarray}&&\mathrm{log}\,{F}_{\mathrm{intr}}=-({0.043}_{-0.051}^{+0.051})\times \mathrm{log}\,{L}_{{\rm{X}}}+({1.645}_{-2.218}^{+2.235}).\end{eqnarray} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

The inset in Figure 2 shows the probability density distribution of α, which we find to be relatively narrow around a peak close to zero (3σ upper limit $\alpha \lt 0.19$). When the missing objects are included, the dependence of F_intr on ${L}_{X}$ is no longer significant (just 80.1% of our simulations gave $\alpha \gt 0$). A fit with constant F_intr yields instead ${F}_{\mathrm{intr}}={58.5}_{-4.4}^{+4.1}$% (see Figure 2).

Since most of the radiation absorbed by the dust in the torus is emitted by the accretion disk at UV/optical wavelengths, to test the validity of receding torus models we should investigate whether F_intr varies with AGN bolometric power, ${L}_{\mathrm{bol}}$. If we assume a constant X-ray bolometric correction, our results remain unchanged. If we adopt instead that the conversion from ${L}_{{\rm{X}}}$ to ${L}_{\mathrm{bol}}$ increases with luminosity (Marconi et al. 2004; Hopkins et al. 2007), using one such bolometric correction (Marconi et al. 2004), our AGN span almost four orders of magnitude in ${L}_{\mathrm{bol}}$ (from ${10}^{43}$ to $7\times {10}^{46}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$) yielding an even tighter limit on the slope, $\alpha ={0.033}_{-0.039}^{+0.040}$ (3σ upper limit $\alpha \lt 0.14$). We find evidence that F_intr does not decrease with increasing luminosity, clearly contradicting the expectations of receding torus models.

In our analysis, we have not considered that the uncertainties in ${L}_{{\rm{X}}}$ might move AGN between luminosity bins. Because the bins are 1 dex wide and the uncertainties in ${L}_{{\rm{X}}}$ small (medians of 5% and 35% for type-1 and type-2 AGN, respectively), the effect should be negligible. A simple calculation of the scatter of sources in the bins by counting sources weighted by the distributions of ${L}_{{\rm{X}}}$ revealed that at ${L}_{{\rm{X}}}\gt {10}^{43}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ our results do not change. At ${L}_{X}\lt {10}^{43}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,,\,{N}_{2}$ increases to ∼9, which translates into an F_intr of ∼67%. This has a negligible impact on our estimate of α (${\rm{\Delta }}\alpha \sim 0.01$).

By filling the gap between studies in the local universe and at cosmic epochs when supermassive black hole (SMBH) mass growth peaked, z ∼ 1–2 (Martinez-Sansigre et al. 2005; Ueda et al. 2014; Aird et al. 2015; Assef et al. 2015; Buchner et al. 2015; Del Moro et al. 2016), we show here that luminous highly obscured AGN dominate the population of fast growing SMBHs up to $z\sim 2$ (see also Reyes et al. 2008). Surveys with the Nuclear Spectroscopic Telescope Array have started to reveal this elusive AGN population. At the z of our sample and luminosities $2\times {10}^{43}\lt {L}_{10-40\mathrm{keV}}\lt 2\times {10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, the observed type-2 AGN fraction is consistent with our findings (${53}_{-15}^{+14} \%$, assuming ${L}_{10-40\mathrm{keV}}$/${L}_{2\mbox{--}10\mathrm{keV}}\approx 1;$ Lansbury et al. 2017).

We have determined the intrinsic type-2 AGN fraction at $0.05\leqslant z\leqslant 1$ and at ${10}^{42}\leqslant {L}_{2\mbox{--}10\mathrm{keV}}\leqslant {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. To do so, we used a complete flux-limited sample of 199 X-ray selected AGN drawn from the BUXS survey. For this sample we have robust estimates of the geometrical covering factor of the torus in the framework of N08 clumpy torus models.

Since the distribution of covering factors needs to match the fraction of optical type-2 AGN, we reveal the existence of a substantial population of X-ray undetected objects with high-covering factor tori, which are increasingly numerous at higher AGN luminosities. When these "missing" objects are included, Compton-thick AGN account at most for ${37.0}_{-10.5}^{+8.9}$% of the total AGN population, in agreement with previous estimates. We find that the intrinsic type-2 AGN fraction is ${58.5}_{-4.4}^{+4.1}$% and has a weak and non-significant (less than 2σ) luminosity dependence. This is in clear contradiction with the results generally reported by AGN surveys and the expectations from receding torus models.

Our findings imply that the majority of luminous, rapidly accreting SMBHs reside in highly obscured nuclear environments, many so deeply buried that they remain undetected in X-rays at the depths of <10 keV wide area surveys.

We thank the anonymous referee for the complete and deep revision of our manuscript that helped us to substantially improve our work. S.M. acknowledges financial support through grant AYA2016-76730-P. (MINECO/FEDER). F.J.C., X.B., and A.A.-H. acknowledge financial support through grant AYA2015-64346-C2-1-P (MINECO/FEDER). A.H.-C. acknowledges financial support through grants AYA2015-70815-ERC and AYA2012-31277. C.R.A. acknowledges financial support through grant AYA2016-76682-C3-2-P and the Ramón y Cajal Program through project RYC-2014-15779 (MINECO). T.M. is supported by CONACyT Grants 179662, 252531 and UNAM-DGAPA PAPIIT IN104216. Based on observations collected at the European Organization for Astronomical Research in the Southern hemisphere, Chile. Based on observations made with the William Herschel Telescope-operated by the Isaac Newton Group, the Telescopio Nazionale Galileo-operated by the Centro Galileo Galilei, and the Gran Telescopio Canarias installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias.