Do current X-ray observations capture most of the black-hole accretion at high redshifts?

Abstract

The cosmic black hole accretion density (BHAD) is critical for our understanding of the formation and evolution of supermassive black holes (BHs). However, at high redshifts ($z>3$), X-ray observations report BHADs significantly ($\sim 10$ times) lower than those predicted by cosmological simulations. It is therefore paramount to constrain the high-$z$ BHAD using independent methods other than direct X-ray detections. The recently established relation between star formation rate and BH accretion rate among bulge-dominated galaxies provides such a chance, as it enables an estimate of the BHAD from the star-formation histories (SFHs) of lower-redshift objects. Using the CANDELS Lyman-$\alpha$ Emission At Reionization (CLEAR) survey, we model the SFHs for a sample of 108 bulge-dominated galaxies at $z=0.7$–1.5, and further estimate the BHAD contributed by their high-$z$ progenitors. The predicted BHAD at $z\approx 4$–5 is consistent with the simulation-predicted values, but higher than the X-ray measurements (by $\approx$3–10 times at $z=$4–5). Our result suggests that the current X-ray surveys could be missing many heavily obscured Compton-thick active galactic nuclei (AGNs) at high redshifts. However, this BHAD estimation assumes that the high-$z$ progenitors of our $z=0.7$–1.5 sample remain bulge-dominated where star formation is correlated with BH cold-gas accretion. Alternatively, our prediction could signify a stark decline in the fraction of bulges in high-$z$ galaxies (with an associated drop in BH accretion). JWST and Origins will resolve the discrepancy between our predicted BHAD and the X-ray results by constraining Compton-thick AGN and bulge evolution at high redshifts.

I Introduction

Understanding the relation between supermassive black holes (BHs) and their host galaxies is one of the most important tasks in extragalactic astronomy. The observations of nearby galaxies reveal a tight correlation with an intrinsic dispersion of $\approx 0.3$ dex between bulge stellar mass ($M_{\star}$) and black hole mass ($M_{\rm BH}$; Kormendy & Ho 2013; Saglia et al. 2016).

Tremendous observational effort has been spent to unveil the origin of this bulge-BH mass relation. One interesting clue comes from the recent observations of Yang et al. (2019), which indicate that the star formation rate (SFR) is linearly correlated with the sample-averaged black hole accretion rate (BHAR) among bulge-dominated galaxies at $z\approx 0.5$–2.5 (see also Kocevski et al. 2017; Ni et al. 2019; Ni et al. 2021 for similar conclusions). Their BHAR/SFR ratio ($\approx 1/300$) is similar to the observed BH-bulge mass ratio in the local universe (Kormendy & Ho 2013). This similarity suggests that the observed BHAR-SFR relation is strongly related to the BH-bulge connection. Also, Yang et al. (2019) found that the BHAR-SFR relation does not hold for galaxies that are not bulge-dominated. This result indicates that BHs only coevolve with bulges rather than the disks, consistent with the observations of the local galaxies (Kormendy & Ho 2013). Ni et al. (2021) found the BHAR-SFR relation also holds for a bulge-dominated sample at lower redshift of $z\lesssim 1.2$. Numerical simulations show that the BHAR-SFR relation is driven by fundamental accretion physics in a bulge-dominated morphological structure (Yao et al. in prep.). Therefore, the BHAR-SFR relation is likely a universal correlation that does not depend on redshift.

The BHAR in Yang et al. (2019) was derived by averaging the X-ray detections or stacked fluxes over samples of sources (hundreds of objects per sample). This averaging process is designed to overcome AGN short-term ($\lesssim 10^{7}$ years) variability and approximate long-term average BH accretion rate (Hickox et al. 2014; Yang et al. 2017; Yuan et al. 2018, e.g.,). We note that the Yang et al. (2019) BHAR is dominated by cold (radiative efficient) accretion rather than hot (radiative inefficient) accretion. This is because the average BHAR is mainly ($\gtrsim 80\%$) contributed by X-ray detected sources rather than stacking, and hot-accretion sources are below the sensitivity of the X-ray data in Yang et al. (2019). The low BH-growth contribution from hot accretion is also expected from simulations (Croton et al. 2006; Yuan et al. 2018, e.g.,). Although it is observationally challenging to constrain the star-formation physical scales among the bulge-dominated galaxies, the simulations of bulge-dominated galaxies suggest that star formation mainly occurs on a nuclear scale of $\lesssim 1$ kpc (Yao et al. in prep.).

One feasible application of the BHAR-bulge SFR relation is to infer the black hole accretion density (BHAD; i.e., total BH accretion rate per comoving volume in units of $M_{\odot}$ yr^-1 Mpc^-3). The BHAD, especially at high redshifts, is an important quantity for the studies of BH formation and evolution (Bonoli et al. 2014; Volonteri et al. 2016, e.g.,). However, the high-redshift BHADs from X-ray observations (Vito et al. 2016; Vito et al. 2018, e.g.,) are significantly lower by a factor of $\sim 10$ than the theoretical results, including predictions from EAGLE (Crain et al. 2015), IllustrisTNG (Weinberger et al. 2017), and Horizon-AGN (Volonteri et al. 2016). If the simulated results are correct, then there is significantly more BH accreted mass at high redshifts and the X-ray surveys are highly incomplete, for example, because of heavy obscuration (i.e., current constraints on BH growth are missing a large fraction of highly obscured AGN at high redshifts). Alternatively, if the X-ray results for BH accretion are correct, then there may be significant flaws in the simulation recipes that lead to the systematic overestimation of BH growth in the early universe.

Therefore, it is paramount to infer the BHAD from another independent method. One possibility is to use constraints on the star-formation histories (SFHs) of the bulge-dominated galaxies combined with the observed BHAR-bulge SFR relation (Yang et al. 2019). Recent work has shown that galaxy SFHs can be robustly constrained by modeling observed spectroscopy and photometric data with stellar population synthesis models (Schreiber et al. 2018; Akhshik et al. 2020; Estrada-Carpenter et al. 2020, e.g.,). Here, we use HST/WFC3 grism and broadband photometric data from the CANDELS Lyman-$\alpha$ Emission At Reionization (CLEAR) survey to constrain the SFHs of bulge-dominated galaxies for this purpose (see, Estrada-Carpenter et al. 2019; Estrada-Carpenter et al. 2020; Simons et al., in prep). The CLEAR fields fall in the GOODS-N and GOODS-S fields, which have deep HST $H$-band imaging, allowing robust selection of bulge-dominated galaxies (Huertas-Company et al. 2015a; Huertas-Company et al. 2015b; Yang et al. 2019).

In this work, we take the SFH results for a sample of bulge-dominated galaxies at $z=0.7$–1.5 in CLEAR. We then estimate the BH accretion histories (BHAHs) of these individual galaxies using the observed BHAR–SFR relation from Yang et al. (2019). We sum the BHAHs and divide it by the comoving volume to estimate the redshift evolution of the cosmic black hole accretion density (BHAD) contributed by the progenitors of our bulge-dominated galaxies. As this BHAD history is based on only the bulge-dominated galaxies at $0.7<z<1.5$, it represents a lower bound on the total BHAD at higher redshift, as it ignores any BH growth in disk/irregular galaxies. We compare our predicted BHAD from the SFHs of bulge-dominated galaxies with the results from direct X-ray observations and cosmological simulations.

In this paper we use the existing constraints on galaxy SFHs and the BHAR-SFR relation for bulge-dominated galaxies to predict the BHAD at high redshifts and we use it to address the discrepancy between simulations (which predict higher BHAR at high redshifts) and X-ray surveys (that measure lower BHAR at high redshifts). The organization of this paper is as follows. We present the CLEAR data analyses and our sample selection in §II. In §III, we describe our procedures of BHAD estimation and compare our BHAD with the simulated and observed BHADs in the literature. In §IV, we discuss the possible uncertainties of different types of BHAD, present physical arguments for our results, and perform sanity checks on our results. We summarize our work and discuss future prospects in §V.

Throughout this paper, we assume a cosmology with $H_{0}=70$ km s^-1 Mpc^-1, $\Omega_{M}=0.3$, and $\Omega_{\Lambda}=0.7$. We adopt a Chabrier initial mass function (IMF; Chabrier 2003). Quoted uncertainties are at the $1\sigma$ (68%) confidence level, unless otherwise stated.

II Data and Sample

The analyses in this work are based on the CLEAR survey (a Cycle 23 HST program, PI: C. Papovich), which has 12 pointings of deep (12 orbit) WFC3 G102 slitless grism spectroscopy in the GOODS-South/North fields. The CLEAR fields also have G141 grism and UV-to-8 $\mu$m photometric data. The detailed data reduction and analyses of the CLEAR data are presented in Estrada-Carpenter et al. (2019), Estrada-Carpenter et al. (2020), and Simons et al. (in prep.). In §II.1, we briefly describe the modeling of CLEAR data which yields galaxy properties such as SFH, redshift, and $M_{\star}$. We then define the sample for this work in §II.2.

II.1 CLEAR data modeling

Figure 1: Example spectral fits, SFH, and BHAH for GS-39170 (left) and GN-12078 (right). The top panels show the best-fit spectra, the G102 grism data, the G141 grism data, and the photometry in black, blue, red, and green, respectively. The middle panels show the the derived SFHs based on the spectral modelling, with uncertainties (inner 68th percentile) indicated by the shaded region. The bottom panels show the BHAHs from the SFHs. The BHAH uncertainties (shaded region) are propagated from both of the SFH and $R$ uncertainties (see §III.1).

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

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 where as 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 towards 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 of we can derive our SFH (as the 50th percentile) and errors on the SFHs (inner 68th percentile). Fig. 1 displays example spectra and SFH fits for two sources in our sample (§II.2). The details of the CLEAR catalog will be presented in Simons et al. (in prep.).

II.2 Sample selection

Among the CLEAR objects, we select bulge-dominated galaxies based on the CANDELS machine-learning $H_{160}$ 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 inspections (see Fig. 2 in Yang et al. 2019) and tend to be compact upon profile fitting (see Fig. 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).

The CLEAR catalog provides redshifts, $M_{\star}$, and SFHs (§II.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 $\mu$m) covers important age and metallicity indicators of, e.g., H$\beta$, Mg$b$, and H$\alpha$. We select all bulge-dominated galaxies with stellar mass above $M_{\star}=10^{10}\ M_{\odot}$, which is the CLEAR mass limit at $z\approx 1.5$. Our volume-limited sample has 108 bulge-dominated objects.

III Estimation of Black Hole Accretion Density

III.1 From SFHs to BHAD

Figure 2: The schematic plot describing our procedures to estimate the BHAD. The three major steps are marked and explained at bottom left. These are discussed in more detail in § III.

Figure 3: Black hole accretion density as a function of Universe age. The black curve indicates our estimated BHAD contributed by bulge-dominated galaxies. The grey shaded region indicates 2$\sigma$ uncertainties. The red and blue curves represent observational (X-ray) and theoretical results from the literature, respectively. All of the X-ray BHADs are derived assuming the same bolometric correction and radiation efficiency as adopted by Yang et al. (2019). The red error bars represent the $1\sigma$ bootstrap uncertainties from Vito et al. (2018). We expect similar uncertainties on the other X-ray BHADs. At high redshifts ($z\approx 4$–5), our BHAD is similar to the simulation predictions, but higher than the X-ray measurements.

Fig. 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 broad-band photometric data (§II). This is “step 1”, labeled (1) in Fig. 2.

In “step 2”, labeled (2) in Fig. 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.,

$${\rm BHAH}_{i}={\rm SFH}_{i}\times R,$$ (1)

where $R=10^{-2.48}$ is the BHAR/SFR ratio derived by Yang et al. (2019, see their Eq. 6) and the subscript ($i$) represents the source index in our bulge-dominated sample. We remind the reader that the BHAH estimated here represent the long-term average accretion rate dominated by cold accretion over the cosmic history (see §I). Fig. 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 Eq. 1, i.e.,

$$\delta{\rm BHAH}_{i}={\rm BHAH}_{i}\sqrt{\left(\frac{\delta{\rm SFH}_{i}}{{\rm SFH}_{i}}\right)^{2}+\left(\frac{\delta R}{R}\right)^{2}},$$ (2)

where $\delta R/R=0.05\times\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 §III.2. The SFH uncertainties are from our modeling of CLEAR data (see §II.1). We apply Eq. 4 to upper and lower uncertainties, respectively, because the upper and lower SFH uncertainties are not always equal (e.g., Fig. 1).

In “step 3”, labeled (3) in Fig. 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.,

$${\rm BHAD}=\frac{\sum_{i=1}^{108}{\rm BHAH}_{i}}{V_{c}}\\ =\frac{\sum_{i=1}^{108}{\rm SFH}_{i}}{V_{c}}R,$$ (3)

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

$$\delta{\rm BHAD}={\rm BHAD}\sqrt{\frac{\sum_{i=1}^{108}(\delta{\rm SFH}_{i})^{2}}{({\sum_{i=1}^{108}\rm SFH}_{i})^{2}}+\left(\frac{\delta R}{R}\right)^{2}}.$$ (4)

The resulting bulge BHAD and its 2$\sigma$ error as a function of redshift/Universe age is displayed in Fig. 3. In this work, we do not extend beyond $z=5$, because there is only $\approx 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 $\approx 1$ Gyr (§II.1, see Estrada-Carpenter et al. 2019; Estrada-Carpenter et al. 2020). We will discuss the SFH/BHAD evolution at $z>5$ in a future dedicated work.

Fig. 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_{\star}>10^{10}\ M_{\odot}$. Since active galactic nuclei (AGNs) can also be found in less massive bulges as well as in non-bulge-dominated galaxies (Yang et al. 2019, e.g.,), our estimated BHAD is an lower bound for the total BHAD, i.e., any BHAD measurement similar 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 Fig. 3).

III.2 BHADs based on literature X-ray luminosity functions

In this section, we derive BHADs based on the measured X-ray luminosity functions (XLFs) from the literature (Ueda et al. 2014; Aird et al. 2015; Vito et al. 2018; Ananna et al. 2019, i.e.,), 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.,

$$\frac{dn}{d\log L_{\rm bol}}=\frac{dn}{d\log L_{\rm X}}\frac{d\log L_{\rm X}}{d\log L_{\rm bol}},$$ (5)

where $dn/d\log L_{\rm bol}$ and $dn/d\log L_{\rm X}$ are the BLF and XLF, respectively, and $d\log L_{\rm X}/d\log L_{\rm 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

$${\rm BHAD}=\frac{1-\epsilon}{\epsilon c^{2}}\int_{43}^{49}L_{\rm bol}\frac{dn}{d\log L_{\rm bol}}d\log L_{\rm bol},$$ (6)

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

In Fig. 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_{\rm H}$), and X-ray fluxes. From Fig. 3, these uncertainties are small compared to the difference between Vito et al. (2018) BHAD and our BHAD at $z\approx 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_{\rm H}\gtrsim 10^{24}$ cm^-2), which are largely missed in X-ray surveys especially at high redshift. We discuss this issue in §IV.

IV Discussion

IV.1 Possible causes of the BHAD discrepancy

Figure 4: Same format as Fig. 4 but showing different cases of the high-redshift ($z>2.5$) evolution of the bulge-dominated galaxy fraction ($f_{\rm bulge}$). The solid curve assumes that $f_{\rm bulge}$ does not evolve at $z>2.5$, following the trend at lower redshifts. The dashed and dotted black curves assume that $f_{\rm bulge}$ decreases at $z>2.5$, following $(1+z)^{-3}$ and $(1+z)^{-5}$, respectively. Compared to the original BHAD (solid black curve), the modulated BHADs (dashed and dotted black curves) are more consistent with the results from X-ray observations at $z\approx 4$–5. Therefore, if the X-ray BHADs are physical, then $f_{\rm bulge}$ must evolve strongly in the early universe.

Fig. 3 shows that our BHAD at $z\approx 4$–5 is similar to theoretical predictions, but is $\approx$3–10 times higher than the X-ray results. This high-$z$ BHAD difference is beyond the expectations from the uncertainties of different BHADs (§III). 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 §I). 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 (Brandt & Alexander 2015, e.g.,). 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 (Hasinger 2008; Liu et al. 2017, e.g.,). 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_{\rm CTK}$) is often assumed to be $\approx$ constant due to the lack of observational constraints, especially at $z\gtrsim 3$ (Ueda et al. 2014; Aird et al. 2015; Buchner et al. 2015; Ananna et al. 2019, e.g.,). The BHAD curves of Ueda et al. (2014), Aird et al. (2015), and Ananna et al. (2019) in Figs. 3 and 4 include the contribution from Compton-thick AGNs based on the assumption that $f_{\rm CTK}$ is constant. These BHAD estimations are still lower at $z\approx 4$–5 than the BHAD we infer from the galaxy SFHs, suggesting that the strong assumption of a constant $f_{\rm 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\approx 4$–5 (Fig. 3). We present some physical arguments for a higher intrinsic BHAD than X-ray observed and make practical predictions for future observations in §IV.2.

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 (Kocevski et al. 2017; Ni et al. 2021, e.g.,), 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_{\rm bulge}$) among the $M_{\star}\approx 10^{11.2}~M_{\odot}$ ($z=0$) galaxies’ progenitors is roughly a constant ($\approx 20\%$–30%) at $z\approx 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\gtrsim 2.5$. Probing this possibility is beyond the capability of current facilities due to the lack of $\gtrsim 1.6\ \mu$m high-resolution imaging, but will be testable with JWST imaging.

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

IV.2 Physical arguments for a higher intrinsic BHAD than X-ray observed

Figure 5: Predicted AGN IR luminosity ($L_{\rm 6\mu m}$) function based on our BHAD at $z=4$ (top) and $z=5$ (bottom). The vertical lines represent $5\sigma$ sensitivity from a typical exposure of 1000 seconds for some IR telescopes as labeled. JWST and Origins will be able to sample around or below the break luminosity ($\sim 10^{44.5}$ erg s^-1).

Our SFH-based BHAD and the simulated BHADs are both higher than the X-ray results at $z\approx 4$–5 (Fig. 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 sub-millimeter surveys begin to reveal a large populations of heavily obscured star-forming galaxies at $z\gtrsim 3$ that are faint/undetected in shorter-wavelengths surveys (González-López et al. 2020; Smail et al. 2021, e.g.,). This new population could contribute a significant (or even dominant) fraction of the cosmic star-formation rate density (SFRD; 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 (§IV.1). From X-ray spectral analyses (Vito et al. 2018; Li et al. 2019, e.g.,), most ($\approx 80\%$–90%) of the X-ray selected $z\gtrsim 3$ AGNs are Compton-thin or unobscured, and none of the detected Compton-thick AGNs have $N_{\rm H}>10^{25}$ cm^-2. This absence of $N_{\rm H}>10^{25}$ cm^-2 AGNs is likely a selection effect due to their X-ray faintness (Hickox & Alexander 2018, e.g.,). These results indicates that current X-ray selections could indeed miss many high-$z$ Compton-thick AGNs (especially at $N_{\rm H}>10^{25}$ cm^-2). Our SFH-based and the simulated BHADs (both dominated by cold accretion; §I) are consistent with this interpretation: if these results are correct then it implies a large population of heavily obscured AGN at $z\gtrsim 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_{\rm X}\sim 10^{42.5}$ erg s^-1 at $z\approx 4$–5; e.g., Vito et al. 2018). Although these missed high-$z$ Compton-thick AGNs have weak or none X-ray signals due to heavy obscuration (Vito et al. 2016, e.g.,), 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 $\log L_{X}=41$–47) equals to unity (see §III.2 for the detailed process). We multiply the normalized XLF by our BHAD (Fig. 3). The above procedure yields a XLF that can produce our BHAD. We then convert this XLF to the AGN IRLF using a relation between $L_{X}$ and $L_{\rm 6\mu m}$ (AGN 6 $\mu$m $\nu L_{\nu}$ luminosity; e.g., Stern 2015), i.e.,

$$\begin{split}\frac{d\Phi}{d\log L_{\rm 6\mu m}}=\frac{d\Phi}{d\log L_{X}}\frac{d\log L_{X}}{d\log L_{\rm 6\mu m}}\\ =\frac{d\Phi}{d\log L_{X}}\times(1.024-0.094\log(L_{\rm 6\mu m}/10^{41}\ {\rm erg\ s^{-1}})).\end{split}$$ (7)

We display the resulting IRLF at $z=4$ and $z=5$ on Fig. 5. We mark sensitivities of some current and future IR missions on Fig. 5. These $L_{\rm 6\mu m}$ limits are converted from the flux-density sensitivities for a typical 1000-second 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; Stalevski et al. 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.

From Fig. 5, Spitzer and Herschel can only sample $L_{\rm 6\mu m}$ more than $\approx 10$ times above the IRLF break luminosity ($L^{*}\sim 10^{44.5}$ erg s^-1). 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 ($\sim 1000$ arcmin²), Spitzer (Herschel) can only detect $\approx 2$ (0) objects according to the IRLF in Fig. 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$\mu$m and the MIPS 24$\mu$m filters (Yang et al. 2021, e.g.,). The future missions of JWST and Origins can sample $\lesssim L^{*}$ objects (see Fig. 5) thanks to their unprecedented sensitivities. Origins perform better than JWST because the AGN SED peak (rest-frame $\approx 5$–20 $\mu$m) is out of the JWST coverage at $z\approx 4$–5. For a CANDELS-like deep survey ($\sim 1000$ arcmin²), JWST (Origins) can detect $\approx 20$ ($\approx 60$) objects according to the IRLF in Fig. 5. Thanks to the continuous wavelength coverage of JWST and Origins, the AGN identification will be practically feasible (Yang et al. 2021, e.g.,).

We caution that the IRLF in Fig. 5 assumes that Compton-thick AGNs follow the same intrinsic (obscuration-corrected) $L_{\rm 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_{\rm X}$ distributions can be inferred from the measured IR luminosities by JWST and Origins.

V Summary and Future Prospects

The quantity of BHAD, despite its importance, is debatable as X-ray measurements are often significantly lower than theoretical predictions, particularly at high redshift ($z\gtrsim 3$) where detections are less complete. In this work, we constrain BHAD at $z=1.5$–5 using a novel method, from a sample of $z=0.7$–1.5 bulge-dominated galaxies (§II). Our BHAD estimation (dominated by cold accretion; §I) is based on the BHAR-SFR correlation among bulge-dominated galaxies (Yang et al. 2019) and the galaxy SFHs derived from their broad-band photometry and grism spectroscopy from HST/WFC3 observations in the CLEAR survey.

Our estimated BHAD is consistent with both of the theoretical predictions and the X-ray measurements at $z\lesssim 3$ (Fig. 3). At $z\approx 4$–5, our BHAD agrees with simulations, but it is higher than the X-ray results (by $\approx$3–10 times at $z=4$–5). After considering several causes of this discrepancy (§IV.1), we argue that it stems from two possibilities. Either (1) there exists a large population of heavily obscured Compton-thick AGN at $z\gtrsim 4$ current missed in X-ray surveys (see §IV.2), or (2) the BHAR-SFR relation begins to break down at $z\gtrsim 2.5$, which could result from a significant drop in the frequency of bulge-dominated galaxies. Either scenario can be tested with future observations.

The high-$z$ Compton-thick AGNs likely have strong IR emission. In the future, JWST and Origins will be able to detect dozens of high-$z$ Compton-thick AGNs in their deep broad-band imaging surveys, if this heavily obscured population is mainly responsible for the discrepancy between our BHAD and the X-ray results (see §IV.2). Another way is to identify high-$z$ AGNs with narrow emission lines (but this requires spectroscopy). For example, the high-$z$ version of BPT diagram (Baldwin et al. 1981), which is often used to classify AGNs vs. star-forming galaxies in the local universe, will be available from JWST observations. There are also some AGN-sensitive lines such as [Ne v] 14.32 $\mu$m and [O iv] 25.88 $\mu$m observable by Origins at $z\approx 4$–5 (Satyapal et al. 2020, e.g.,). It will be interesting to further study the Lyman continuum escape fraction ($f_{\rm esc}$) of the JWST/Origins-detected Compton-thick AGNs. We expect that $f_{\rm esc}$ to be low (nearly zero) considering the strong obscuration in X-ray. But if this is not the case, then the high-$z$ Compton-thick AGNs could be an important source of cosmic reionization (Fan et al. 2006; Robertson et al. 2015; Finkelstein et al. 2019, e.g.,).

JWST will also be able to test the frequency of bulge-dominated galaxies at high redshifts. Our BHAD estimation assumes that the high-$z$ progenitors of our sample remain bulge-dominated (§IV.1). This assumption could be impacted by morphological transformation, although observations find such transformation is unlikely prevalent at $z\lesssim 2.5$. The currently available HST $H$-band imaging is shifted into rest-frame UV wavelengths at $z\gtrsim 2.5$, preventing reliable morphological classifications (Conselice 2014; Huertas-Company et al. 2015a, e.g.,). JWST will overcome this issue by providing high-resolution imaging of wavelengths up to $\approx 5\ \mu$m (NIRCam). If morphological transformation is responsible for our reported BHAD difference, JWST will find that bulge-dominated galaxies are very rare, less than a few percent among massive galaxies at $z\approx 5$ (Fig. 4).

In reality, we expect that the Universe will surprise us. For example, it is reasonable to hypothesize that multiple effects may be at play (including others not considered here). We may discover both a higher abundance of obscured AGN, evolution in the bulge-fraction of galaxies, and/or something entirely unexpected. Regardless, the predictions are important as they provide a baseline, and then future studies will lead to an improved understanding of the history of BH accretion and the joint evolution of BH accretion and star-formation in galaxies.