Do Current X-Ray Observations Capture Most of the Black-hole Accretion at High Redshifts?

The methodology to derive the stellar populations and SFHs for our galaxies is outlined in Estrada-Carpenter et al. (2019, 2020). For each galaxy we derive posteriors for stellar population properties such as stellar metallicity, SFH span, SFH bins, redshift, A(V), and M_⋆. ⁶

For our SFHs we used the flexible approach outlined in Leja et al. (2019), wherein SFHs are modeled using a set of time bins and the mass generated in each time bin is fit for. The maximum span of our SFHs is determined by redshift and is set to be the age of the universe. We allow the overall span of the time bins to vary, which produces smoother SFHs, whereas keeping the bins set would produce a step-wise SFH. The amount of time bins used depends on the UVJ color–color diagram classification of the galaxy, where quiescent selected galaxies use 10 time bins, and star-forming use 6. This was done because quiescent galaxies form most of their mass early on and using more time bins would allow for a higher temporal resolution at early formation times, while star-forming galaxies will tend to form more of their mass later; therefore, there is less of a need for higher temporal resolution at early times (this choice also reduces the run time of our SED fits). For our sample, we use a continuity prior (Leja et al. 2019). This prior weighs toward a more continuous SFH and against a bursty history.

The resulting SFHs were derived by sampling the posteriors of the SFH span, SFH bins, and stellar mass and generating 5000 iterations of SFHs. From this sample, we can derive our SFH (as the 50th percentile) and errors on the SFHs (inner 68th percentile). Figure 1 displays example spectra and SFH fits for two sources in our sample (Section 2.2). The details of the CLEAR catalog will be presented in R. C. Simons et al. (in preparation).

Zoom In Zoom Out Reset image size

Among the CLEAR objects, we select bulge-dominated galaxies based on the CANDELS machine-learning H₁₆₀ morphological classifications (Huertas-Company et al. 2015a), following the same criterion as in Yang et al. (2019). ⁷ The machine-learning-selected bulge-dominated galaxies have round and smooth shapes upon visual inspection (see Figure 2 in Yang et al. 2019) and tend to be compact upon profile fitting (see Figure C4 of Ni et al. 2021). It is critical to select bulge-dominated galaxies, as BHAR is only correlated with the SFR in bulge-dominated galaxies not the SFR in other types of galaxies (Yang et al. 2019).

Zoom In Zoom Out Reset image size

The CLEAR catalog provides redshifts, M_⋆, and SFHs (Section 2.1). We focus on the redshift range of z = 0.7–1.5. This redshift range guarantees that the grism spectroscopy (0.8–1.7 μm) covers important age and metallicity indicators of, e.g., Hβ, Mgb, and Hα. We select all bulge-dominated galaxies with stellar mass above M_⋆ = 10¹⁰ M_⊙, which is the CLEAR mass limit at z ≈ 1.5. Our volume-limited sample has 108 bulge-dominated objects.

Figure 2 shows the schematic for the steps we applied to convert the SFHs of individual bulge-dominated galaxies to an estimate of the BHAD. We take the SFH and the associated uncertainties for each galaxy derived by modeling the grism spectroscopy and broadband photometric data (Section 2). This is "step 1," labeled (1) in Figure 2.

In "step 2," labeled (2) in Figure 2, we multiply the SFHs by the BHAR–SFR relation from Yang et al. (2019). This procedure yields a BH accretion history (BHAH) for each object in our sample, i.e.,

$$\begin{eqnarray}&&{\mathrm{BHAH}}_{i}={\mathrm{SFH}}_{i}\times R,\end{eqnarray} \tag{ 1 }$$

where R = 10^−2.48 is the BHAR/SFR ratio derived by Yang et al. (2019), see their Equation (6) and the subscript (i) represents the source index in our bulge-dominated sample. We remind the reader that the BHAH estimated here represents the long-term average accretion rate dominated by cold accretion over the cosmic history (see Section 1). Figure 1 shows two example BHAHs and the associated uncertainties. The BHAH uncertainties are propagated from both of the SFH and R uncertainties, using the standard error-propagation formula based on Equation (1), i.e.,

$$\begin{eqnarray}&&\delta {\mathrm{BHAH}}_{i}={\mathrm{BHAH}}_{i}\sqrt{{\left(\displaystyle \frac{\delta {\mathrm{SFH}}_{i}}{{\mathrm{SFH}}_{i}}\right)}^{2}+{\left(\displaystyle \frac{\delta R}{R}\right)}^{2}},\end{eqnarray} \tag{ 2 }$$

where $\delta R/R=0.05\times \mathrm{ln}10=0.12$ is the fitting uncertainty in Yang et al. (2019), under the assumptions of radiation efficiency and bolometric correction. We address the systematic uncertainties arising from these assumptions in Section 3.2. The SFH uncertainties are from our modeling of CLEAR data (see Section 2.1). We apply Equation (4) to upper and lower uncertainties, respectively, because the upper and lower SFH uncertainties are not always equal (e.g., Figure 1).

In "step 3," labeled (3) in Figure 2, we sum the BHAHs for the 108 bulge-dominated galaxies and divide it by the comoving volume (V_c ) of z = 0.7–1.5 covered by the CANDELS/CLEAR area (61 arcmin²), i.e.,

$$\begin{eqnarray}&&\mathrm{BHAD}=\displaystyle \frac{\sum _{i=1}^{108}{\mathrm{BHAH}}_{i}}{{V}_{c}}=\displaystyle \frac{\sum _{i=1}^{108}{\mathrm{SFH}}_{i}}{{V}_{c}}R,\end{eqnarray} \tag{ 3 }$$

where we apply Equation (1). We propagate the SFH and R uncertainties into the BHAD using the standard error-propagation formula, i.e.,

$$\begin{eqnarray}&&\delta \mathrm{BHAD}=\mathrm{BHAD}\sqrt{\displaystyle \frac{{\sum }_{i=1}^{108}{(\delta {\mathrm{SFH}}_{i})}^{2}}{{\left({\sum }_{i=1}^{108}{\mathrm{SFH}}_{i}\right)}^{2}}+{\left(\displaystyle \frac{\delta R}{R}\right)}^{2}}.\end{eqnarray} \tag{ 4 }$$

The resulting bulge BHAD and its 2σ error as a function of redshift/universe age is displayed in Figure 3. In this work, we do not extend beyond z = 5, because there is only ≈1 Gyr cosmic time at z > 5 and our SFH measurements may not have the time resolution sufficiently high to probe the detailed SFH evolution in the first ≈1 Gyr (Section 2.1; see Estrada-Carpenter et al. 2019, 2020). We will discuss the SFH/BHAD evolution at z > 5 in a future dedicated work.

Zoom In Zoom Out Reset image size

Figure 3 also displays theoretical BHADs from cosmological simulations. The Horizon-AGN BHAD curve (Volonteri et al. 2016) is compiled by Vito et al. (2018). The EAGLE (Crain et al. 2015; Schaye et al. 2015) and IllustrisTNG (Weinberger et al. 2017; Pillepich et al. 2018) results are from the corresponding database, where we use the simulation sets of "RefL0100N1504" (EAGLE) and "TNG100-1" (IllustrisTNG). When deriving the simulated BHAD at a given redshift, we add up the BHARs from different galaxies and divide them by the simulated comoving volume.

We remind the reader that our BHAD only accounts for BH growth in bulge-dominated galaxies with M_⋆ > 10¹⁰ M_⊙. Since active galactic nuclei (AGNs) can also be found in less massive bulges as well as in non-bulge-dominated galaxies (e.g., Yang et al. 2019), our estimated BHAD is a lower bound for the total BHAD, i.e., any BHAD measurement similar to or above our values should be considered as consistent with our estimation. Therefore, our BHAD is consistent with the simulated results at z = 1.5–5 in general (see Figure 3).

In this section, we derive BHADs based on the measured X-ray luminosity functions (XLFs) from the literature (i.e., Ueda et al. 2014; Aird et al. 2015; Vito et al. 2018; Ananna et al. 2019), and compare the results with our SFH-based BHAD. We do not use the BHADs from the literature directly, because those are estimated based on different assumptions of radiation efficiencies and bolometric corrections. Below, when deriving the XLF-based BHADs, we adopt the same radiation efficiency and bolometric correction used by Yang et al. (2019) for the BHAR–SFR relation. In this way, we effectively address the systematic uncertainties due to radiation efficiency and bolometric correction.

To derive the BHAD from each XLF, we first convert the XLF to the bolometric luminosity function (BLF), i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{d\mathrm{log}{L}_{\mathrm{bol}}}=\displaystyle \frac{{dn}}{d\mathrm{log}{L}_{{\rm{X}}}}\displaystyle \frac{d\mathrm{log}{L}_{{\rm{X}}}}{d\mathrm{log}{L}_{\mathrm{bol}}},\end{eqnarray} \tag{ 5 }$$

where ${dn}/d\mathrm{log}{L}_{\mathrm{bol}}$ and ${dn}/d\mathrm{log}{L}_{{\rm{X}}}$ are the BLF and XLF, respectively, and $d\mathrm{log}{L}_{{\rm{X}}}/d\mathrm{log}{L}_{\mathrm{bol}}$ is the derivative of the luminosity-dependent bolometric correction from Hopkins et al. (2007), which is also adopted by Yang et al. (2019). ⁸ From the BLF, we can calculate the BHAD by

$$\begin{eqnarray}&&\mathrm{BHAD}=\displaystyle \frac{1-\epsilon }{\epsilon {c}^{2}}{\int }_{43}^{49}{L}_{\mathrm{bol}}\displaystyle \frac{{dn}}{d\mathrm{log}{L}_{\mathrm{bol}}}d\mathrm{log}{L}_{\mathrm{bol}},\end{eqnarray} \tag{ 6 }$$

where c is the speed of light and is the radiation efficiency, and the integral limits 43 and 49 $[\mathrm{log}(\mathrm{erg}\ {{\rm{s}}}^{-1})]$ correspond to $\mathrm{log}{L}_{X}=42$ and 47 $[\mathrm{log}(\mathrm{erg}\ {{\rm{s}}}^{-1})]$ under the Hopkins et al. (2007) bolometric correction. We adopt = 0.1, which is the value used by Yang et al. (2019). Figure 3 displays the resulting XLF-based BHADs. These BHADs are consistent with our BHAD at low redshifts (z ≲ 3), but are lower than our values by a factor of ≈3–10 at z ≈ 4–5.

In Figure 3, we also show the uncertainties of Vito et al. (2018) BHAD. These uncertainties were calculated by Vito et al. (2018) based on binned high-z AGNs, employing a bootstrap technique. This bootstrap technique properly accounts for statistical fluctuations due to limited AGN sample sizes at different luminosities. It also propagates different types of errors, including photometric redshift, column density (N_H), and X-ray fluxes. From Figure 3, these uncertainties are small compared to the difference between Vito et al. (2018) BHAD and our BHAD at z ≈ 4–5. Therefore, the XLF uncertainties cannot explain the discrepancy between our BHAD and the X-ray results. The uncertainties of the other XLF-based BHADs (Ueda et al. 2014; Aird et al. 2015; Ananna et al. 2019), although not publicly available, ⁹ should be comparable to the Vito et al. (2018) uncertainties at high redshifts. This is because the BHAD uncertainties in all of the X-ray works are dominated by the relatively small number of high-z AGNs detected in the existing deep X-ray surveys (e.g., CDF-S and CDF-N). However, we caution that all of the XLF works could suffer from significant systematic uncertainties due to Compton-thick AGNs (N_H ≳ 10²⁴ cm⁻²), which are largely missed in X-ray surveys especially at high redshift. We discuss this issue in Section 4.

Figure 3 shows that our BHAD at z ≈ 4–5 is similar to theoretical predictions, but is ≈3–10 times higher than the X-ray results. This high-z BHAD difference is beyond the expectations from the uncertainties of different BHADs (Section 3). We note that hot radiative-inefficient accretion should not be responsible for this BHAD discrepancy, because both of our BHAD and the X-ray BHADs should be dominated by cold accretion (see Section 1). We discuss some possible causes of the discrepancy below.

X-ray observation is often reliable in AGN selection owing to the nearly universal X-ray emission from AGN and the weak contamination from host galaxies (e.g., Brandt & Alexander 2015). However, AGNs could be missed by X-ray surveys due to strong obscuration. Vito et al. (2016) performed an X-ray stacking analysis for high-z X-ray undetected galaxies using the deepest X-ray survey, CDF-S (Luo et al. 2017). They found the stacked X-ray emission is negligible compared to that from X-ray detected AGNs. Their result suggests that, if a large number of high-z AGNs are missed, they must be heavily obscured, likely at the Compton-thick level (e.g., Hickox & Alexander 2018).

AGN obscuration generally increases toward high redshift (e.g., Hasinger 2008; Liu et al. 2017). The fraction of obscured AGNs (both Compton-thick and Compton-thin AGN) has been modeled as a positive function of redshift, but the ratio of Compton-thick and Compton-thin AGNs (f_CTK) is often assumed to be ≈ constant due to the lack of observational constraints, especially at z ≳ 3 (e.g., Ueda et al. 2014; Aird et al. 2015; Buchner et al. 2015; Ananna et al. 2019). The BHAD curves of Ueda et al. (2014), Aird et al. (2015), and Ananna et al. (2019) in Figures 3 and 4 include the contribution from Compton-thick AGNs based on the assumption that f_CTK is constant. These BHAD estimations are still lower at z ≈ 4–5 than the BHAD we infer from the galaxy SFHs, suggesting that the strong assumption of a constant f_CTK may be incorrect, especially at high redshifts. Therefore, the population of Compton-thick AGN may be dominant over the Compton-thin population at high redshifts. This is also supported by simulations, since the simulated BHADs are also higher than the X-ray BHADs at z ≈ 4–5 (Figure 3). We present some physical arguments for a higher intrinsic BHAD than X-ray observed and make practical predictions for future observations in Section 4.2.

Zoom In Zoom Out Reset image size

In our calculation of the BHAD, we assume that the BHAR–SFR relation still holds for the progenitors of our bulge-dominated galaxies. Our BHAD prediction would be overestimated if many of the galaxy progenitors were non-bulge-dominated at earlier times (e.g., Kocevski et al. 2017; Ni et al. 2021), which could be a result of morphological transformation in star-forming disk galaxies, caused by, e.g., major mergers. As one counter example, Huertas-Company et al. (2015b) found that the bulge-dominated fraction (f_bulge) among the M_⋆ ≈ 10^11.2 M_⊙ (z = 0) galaxies' progenitors is roughly a constant (≈20%–30%) at z ≈ 0–2.5, suggesting that morphological transformation between bulge-dominated and other types is not prevalent at least at these redshifts (see, e.g., Mortlock et al. 2013 and Conselice 2014, who arrive at similar conclusions). However, it is still possible that such transformation happens frequently at z ≳ 2.5. Probing this possibility is beyond the capability of current facilities due to the lack of ≳ 1.6 μm high-resolution imaging, but will be testable with JWST imaging.

We quantify the effects of possible high-redshift f_bulge evolution based on the assumption that our BHAD declines following BHAD ∝ f_bulge ∝ (1 + z)^−γ at z > 2.5. In Figure 4 we plot the BHAD evolution curves in Figure 4 for different values of f_bulge. From Figure 4, to match the BHAR inferred from X-ray surveys would require a very steep power-law index of γ ≈ 3–5. This means that, if the X-ray BHADs are correct, then f_bulge must decline by a large factor of ≈10–50 from z ≈ 2.5 to z ≈ 5. This provides a testable prediction for JWST observations of galaxies at these redshifts.

Our SFH-based BHAD and the simulated BHADs are both higher than the X-ray results at z ≈ 4–5 (Figure 3). It is understandable that simulations predict a relatively strong BH accretion process at high redshifts, because cold gas, which fuels both star formation and AGN, is likely abundant and concentrated in the early universe.

At high redshifts, large amounts of dust associated with the gas can totally obscure the (rest-frame) UV light from intensive star formation. The recent development of submillimeter surveys begin to reveal large populations of heavily obscured star-forming galaxies at z ≳ 3 that are faint/undetected in shorter-wavelength surveys (e.g., González-López et al. 2020; Smail et al. 2021). This new population could contribute a significant (or even dominant) fraction of the cosmic SFR density (e.g., Wang et al. 2019; Gruppioni et al. 2020).

Likewise, the gas/dust-rich environment at high redshifts could also obscure high-z AGN activity. If Compton-thick obscuration is common at high redshifts, then the X-ray selected AGNs will be highly incomplete, leading to an underestimation of BHAD (Section 4.1). From X-ray spectral analyses (e.g., Vito et al. 2018; Li et al. 2019), most (≈80%–90%) of the X-ray selected z ≳ 3 AGNs are Compton-thin or unobscured, and none of the detected Compton-thick AGNs have N_H > 10²⁵ cm⁻². This absence of N_H > 10²⁵ cm⁻² AGNs is likely a selection effect due to their X-ray faintness (e.g., Hickox & Alexander 2018). These results indicate that current X-ray selections could indeed miss many high-z Compton-thick AGNs (especially at N_H > 10²⁵ cm⁻²). Our SFH-based and the simulated BHADs (both dominated by cold accretion; Section 1) are consistent with this interpretation: if these results are correct then it implies a large population of heavily obscured AGN at z ≳ 3 than currently found in X-ray surveys.

Since the missed AGNs are undetected by the currently deepest X-ray surveys (i.e., CDF-S and CDF-N), their apparent (uncorrected for obscuration) luminosities must lie below the survey sensitivity (L_X ∼ 10^42.5 erg s⁻¹ at z ≈ 4–5; e.g., Vito et al. 2018). These missed high-z Compton-thick AGNs have weak or none X-ray signals due to heavy obscuration (e.g., Vito et al. 2016), but they are likely luminous at IR wavelengths due to dust re-emission. Therefore, they can be detected by IR telescopes. It is thereby useful to quantitatively predict the high-z AGN IR luminosity function (IRLF) based on our SFH-based BHAD.

To perform this task, we first take the XLF from Aird et al. (2015). We then normalize the XLF at a given redshift so that the corresponding BHAD (integrated over $\mathrm{log}{L}_{X}=41\mbox{--}47$) equals unity (see Section 3.2 for the detailed process). We multiply the normalized XLF by our BHAD (Figure 3). The above procedure yields an XLF that can produce our BHAD. We then convert this XLF to the AGN IRLF using a relation between L_X and L_6μm (AGN 6 μm ν L_ν luminosity; e.g., Stern 2015), i.e.,

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d{\rm{\Phi }}}{d\mathrm{log}{L}_{6\mu {\rm{m}}}}=\displaystyle \frac{d{\rm{\Phi }}}{d\mathrm{log}{L}_{X}}\displaystyle \frac{d\mathrm{log}{L}_{X}}{d\mathrm{log}{L}_{6\mu {\rm{m}}}}\\ =\ \displaystyle \frac{d{\rm{\Phi }}}{d\mathrm{log}{L}_{X}}\times (1.024-0.094\mathrm{log}({L}_{6\mu {\rm{m}}}/{10}^{41}\ \mathrm{erg}\ {{\rm{s}}}^{-1})).\end{array}\end{eqnarray} \tag{ 7 }$$

We display the resulting IRLF at z = 4 and z = 5 on Figure 5. We mark sensitivities of some current and future IR missions on Figure 5. These L_6μm limits are converted from the flux-density sensitivities for a typical 1000 s exposure, assuming a K correction based on an AGN IR spectral template generated by x-cigale (Boquien et al. 2019; Yang et al. 2020). x-cigale employs a clumpy torus model, skirtor (Stalevski et al. 2012, 2016). We set the viewing angle to 70°, typical for obscured AGNs (Yang et al. 2020), and leave other parameters as the default values. Our conclusion below is not sensitive to the IR model parameters in x-cigale.

Zoom In Zoom Out Reset image size

From Figure 5, Spitzer and Herschel can only sample L_6μm more than ≈10 times above the IRLF break luminosity (L^∗ ∼ 10^44.5 erg s⁻¹). This means that Spitzer and Herschel are not able to effectively detect the predicted Compton-thick AGNs, as the IRLF declines sharply above L^∗. For a CANDELS-like deep survey (∼1000 arcmin²), Spitzer (Herschel) can only detect ≈2 (0) objects according to the IRLF in Figure 5. Also, the task of AGN identification is challenging for Spitzer, as there is a large wavelength "gap" between the coverages of the IRAC 8 μm and the MIPS 24 μm filters (e.g., Yang et al. 2021). The future missions of JWST and Origins can sample ≲ L^∗ objects (see Figure 5) thanks to their unprecedented sensitivities. Origins performs better than JWST because the AGN SED peak (rest-frame ≈5–20 μm) is out of the JWST coverage at z ≈ 4–5. For a CANDELS-like deep survey (∼1000 arcmin²), JWST (Origins) can detect ≈20 (≈60) objects according to the IRLF in Figure 5. Thanks to the continuous wavelength coverage of JWST and Origins, the AGN identification will be practically feasible (e.g., Yang et al. 2021).

We caution that the IRLF in Figure 5 assumes that Compton-thick AGNs follow the same intrinsic (obscuration-corrected) L_X distribution as Compton-thin and unobscured AGNs at high redshifts. This assumption can be tested in the future using the Compton-thick samples detected by JWST and Origins as above. If this assumption turns out to be incorrect, the Compton-thick intrinsic L_X distributions can be inferred from the measured IR luminosities by JWST and Origins.