The Clustering of Undetected High-redshift Black Holes and Their Signatures in Cosmic Backgrounds

TL;DR

This paper develops a new model to analyze the contribution of high-redshift black holes to cosmic background fluctuations, finding that early black hole growth cannot explain the observed signals, and ruling out alternative galactic foreground explanations.

Contribution

The paper introduces a novel framework combining semi-analytic black hole growth models with hydrodynamical spectra to assess their impact on cosmic background fluctuations.

Findings

High-redshift black holes at z>6 cannot account for the fluctuations due to excessive early growth.

Black holes at z~2-3 produce poor fits and overestimate local black hole mass density.

Galactic warm dust foregrounds are unlikely to cause the observed fluctuations.

Abstract

There exist hitherto unexplained fluctuations in the Cosmic Infrared Background (CIB) on arcminute scales and larger. These have been shown to cross-correlate with the Cosmic X-ray Background (CXB), leading several authors to attribute the excess to a high-redshift growing black hole population. In order to investigate potential sources that could explain this excess, in this paper, we develop a new framework to compute the power spectrum of undetected sources that do not have constant flux as a function of halo mass. In this formulation, we combine a semi-analytic model for black hole growth and their simulated spectra from hydrodynamical simulations. Revisiting the possible contribution of a high-redshift black hole population, we find that too much black hole growth is required at early epochs for z>6 accretion to explain these fluctuations. Examining a population of accreting black holes at more moderate redshifts, z\sim 2-3, we find that such models produce a poor fit to the observed fluctuations while simultaneously overproducing the local black hole mass density. Additionally, we rule out the hypothesis of a missing Galactic foreground of warm dust that produces coherent fluctuations in the X-ray via reflection of Galactic X-ray binary emission. Although we firmly rule out accreting massive black holes as the source of these missing fluctuations, additional studies will be required to determine their origin.

Figures

Tables

Table 1. Table 1: Integrated mass densities of central black holes at z = 6 𝑧 6 z=6 predicted by the semi-analytic model. The total mass density (left) includes seed masses, while the accreted mass density (right) does not.

Seed Type	$\log {\rho }_{\mathrm{tot}}\left(6\right)\left[M_{\odot }{\mathrm{Mpc}}^{-3}\right]$	$\log {\rho }_{\mathrm{acc}}\left(6\right)\left[M_{\odot }{\mathrm{Mpc}}^{-3}\right]$
Heavy	$3.9755\pm 0.0098$	$3.429\pm 0.007$
Light	$3.685\pm 0.003$	$3.674\pm 0.003$

Table 2. Table 2: Best fit parameters for an unknown extra-galactic population of AGN with a characteristic host halo mass of 10 13 ​ M ⊙ superscript 10 13 subscript M direct-product 10^{13}\ \mathrm{M}_{\odot} . F tot , 3.6 subscript 𝐹 tot 3.6 F_{\mathrm{tot},3.6} has units of erg ​ s − 1 ​ cm − 2 ​ ster − 1 erg superscript s 1 superscript cm 2 superscript ster 1 \mathrm{erg}\;\mathrm{s}^{-1}\;\mathrm{cm}^{-2}\;\mathrm{ster}^{-1} . Note that the data require an extreme 1:100 flux ratio between the soft X-ray and the 3.6 μ ​ m 𝜇 𝑚 \mu m bands.

	$\log \left(F_{\mathrm{tot},3.6}\right)$	$\log \left(F_{3.6}/F_{4.5}\right)$	$\log \left(F_x/F_{3.6}\right)$	$\log C_{X-\mathrm{IR}}$
AGN_z2	$-{6.317}_{-0.011}^{+0.010}$	$-{0.063}_{-0.013}^{+0.014}$	$-{2.09}_{-0.46}^{+0.04}$	$-{0.6}_{-0.7}^{+0.2}$
AGN_z3	$-{6.438}_{-0.011}^{+0.010}$	$-{0.064}_{-0.013}^{+0.014}$	$-{2.08}_{-0.38}^{+0.05}$	$-{0.6}_{-0.6}^{+0.2}$

Table 3. Table 3: Fitting metrics for various alternate models. Generic extra-galactic populations do not provide good fits due to the gentle slope in the CIB data on arcminute scales. A simple power law is highly preferred over an extra-galactic population on the basis of the Bayesian Information Criterion (BIC).

	${\chi }^2$	${\chi }^2/\mathrm{DOF}$	BIC
Power Law	63.6	1.27	83.6
AGN_z2	1140	22.4	1160
AGN_z3	1090	21.5	1110

Table 4. Table 4: Best fit parameters for an unknown foreground component. Note that while a flux ratio between the CXB and CIB of ∼ 0.01 similar-to absent 0.01 \sim 0.01 is preferred, the large error bars on large scales in this dataset allow for little cross-power between the CIB and CXB. P 1000 subscript 𝑃 1000 P_{1000} has units of erg 2 ​ s − 2 ​ cm − 4 ​ ster − 1 superscript erg 2 superscript s 2 superscript cm 4 superscript ster 1 \mathrm{erg}^{2}\;\mathrm{s}^{-2}\;\mathrm{cm}^{-4}\;\mathrm{ster}^{-1} .

$\log P_{1000}$	$\alpha$	$\log \left(F_{4.5}/F_{3.6}\right)$	$\log \left(F_x/F_{3.6}\right)$	$\log C_{X-\mathrm{IR}}$
$-{20.30}_{-0.058}^{+0.038}$	$-{1.47}_{-0.03}^{+0.04}$	$-{0.105}_{-0.010}^{+0.010}$	$-{1.9778}_{-5.6118}^{+0.0002}$	$-{0.7}_{-3.8}^{+0.3}$

Equations

$$\eta (F_x,F_{IR})=\{\begin{array}{ll}0& \mathrm{if}{F}_{IR}<25,\\ \zeta (F_x)& \mathrm{otherwise}\end{array}$$

$$P_{\mathrm{tot}}=P(q)+P_{\mathrm{SN}},$$

$$P(q)=\int \frac{H(z)}{c{d}_c^2(z)}\frac{d{F}_i}{dz}\frac{d{F}_j}{dz}{P}_{i,j}(q{d}_c^{-1},z)dz,$$

$$\frac{d{F}_i}{dz}=\int \frac{d^{3}N}{d{F}_{i}d\Omega dz}\eta (F_x,F_{IR})d{F}_i.$$

$$\frac{d^{3}N}{d{F}_{i}d\Omega dz}=\frac{d^{2}V}{dzd\Omega }\frac{d^{2}N}{d{F}_{i}dV}.$$

$$\frac{d^{2}V}{dzd\Omega }=c{d}_L^2(1+z{)}^{-1}\frac{dt}{dz}.$$

$$P(k)=P^{1h}(k)+P^{2h}(k)$$

$$F^c=⟨F^c(M_h,z)\eta (F_x,F_{IR})⟩,$$

$$F^s=⟨F^s(M_h,z)\eta (F_x,F_{IR})⟩.$$

$$\bar{F}=\int \frac{dn}{d{M}_h}\left[F^c+F^s\right]d{M}_h.$$

$$B^c=\int \frac{dn}{d{M}_h}\frac{F^c}{\bar{F}}b(M_h,z)d{M}_h$$

$$B^s=\int \frac{dn}{d{M}_h}\frac{F^s}{\bar{F}}b(M_h,z)u(k|M_h,z)d{M}_h$$

$$P_{i,j}^{1h}(k)=\int \frac{dn}{d{M}_h}\times$$

$$\frac{(F_i^s{F}_j^c+F_j^s{F}_i^c)u(k|M_h,z)+F_i^s{F}_j^s{u}^2(k|M_h,z)}{{\bar{F}}_i{\bar{F}}_j}d{M}_h,$$

$$P_{i,j}^{2h}(k)=P^{\mathrm{lin}}(k,z)(B_i^c+B_i^s)(B_j^c+B_j^s),$$

$$b_{\mathrm{eff}}(z)=\int \frac{dn}{d{M}_h}\frac{F^c+F^s}{\bar{F}}b(M_h,z)d{M}_h.$$

$$P_{\mathrm{SN}}=\int \int F_i{F}_j\frac{d^{2}N}{d{F}_{i}d{F}_j}\eta (F_x,F_{IR})d{F}_{i}d{F}_j$$

$$F_{\mathrm{MW}}\sim L_{\mathrm{MW}}\frac{1}{\pi (15\mathrm{kpc}{)}^2}\frac{1}{4\pi }\sim 3\times 10^{-8}\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}{\mathrm{ster}}^{-1}.$$

$$F_{\mathrm{ref}}=F_{\mathrm{MW}}(1-e^{-N_H{\sigma }_T}),$$

$$F_{\mathrm{ref}}\sim 2\times 10^{-12}\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}{\mathrm{ster}}^{-1}.$$

$$\frac{dF}{d\Omega }=\frac{d^{3}E}{dtdAd\Omega },$$

$$\frac{dF}{dV}=\frac{f_{\mathrm{bol}}}{4\pi d_L^2}\frac{dL}{dV},$$

$$\frac{dF}{d\Omega }=\int \frac{dF}{dV}\frac{d^{2}V}{d\Omega dz}dz.$$

$${\rho }_{\bullet }=\frac{1-ϵ}{ϵc^2}\int \frac{dL}{dV}\frac{dt}{dz}dz.$$

$${\rho }_{\mathrm{extra}}=\frac{1-ϵ}{ϵ}\frac{4\pi (1+z_{\mathrm{eff}})}{f_{\mathrm{bol}}{c}^3}\frac{dF}{d\Omega },$$

$$1+z_{\mathrm{eff}}=\int(1+z)\frac{dF}{dV}\frac{d^{2}V}{d\Omega dz}dz\bigg{/}\int\frac{dF}{dV}\frac{d^{2}V}{d\Omega dz}dz\,.$$

$${\rho }_{\bullet }=10^{6.0}{\mathrm{M}}_{\odot }{\mathrm{Mpc}}^{-3},$$

Full text

††pagerange: LABEL:firstpage–LABEL:lastpage††pubyear: 2019

1 Introduction

The Cosmic Infrared Background (CIB) and Cosmic X-ray Background (CXB) contain the integrated emission of all photons in the universe from starlight and accretion powering active galactic nuclei (AGN) (Salvaterra & Ferrara, 2003; Cooray et al., 2004; Kashlinsky et al., 2004). Both the CIB and the CXB are composed of a resolved component detectable in deeper surveys, and an unresolved component arising from sources that are not individually detected. The unresolved component of the CIB is believed to originate primarily from undetected low luminosity galaxies (Kashlinsky et al., 2012; Cooray et al., 2012; Yue et al., 2013a), while the unresolved CXB likely originates from faint populations of as yet undetected AGN and star-forming galaxies (Gilli et al., 2007). This is inferred by evaluating the putative contribution from extrapolations of the following known source populations: star-forming galaxies, AGN, and the hot diffuse gas in groups and clusters.

In the near-infrared, the absolute intensity of the CIB is difficult to infer due to the dominant contribution of zodiacal light and the uncertain contribution of Galactic cirrus (Leinert et al., 1998). While zodiacal light dominates the total intensity, it is fortunately thought to be very structurally smooth (Kelsall et al., 1998). Therefore, studying CIB fluctuations provides an alternative approach to determine its origin (Kashlinsky et al., 1996). Fluctuation analyses of the measured CIB—its power spectrum or two-point correlation function—have revealed excess power on arcminute to degree scales on the sky that remain unexplained even with the most optimistic extrapolations of faint galaxy populations (Helgason et al., 2012; Yue et al., 2013a), pointing to an entirely new missing contribution to the inventory. This excess has now been seen in both the 3.6 and 4.5 micron bands by two different infrared telescopes (Kashlinsky, 2005; Kashlinsky et al., 2012; Matsumoto et al., 2011), is detected in multiple independent fields, and has persisted through years of improved measurements (see Kashlinsky et al., 2018, for a review).

Analyses of the missing fluctuations’ electromagnetic spectrum offer some additional clues to their origin. Recently, these excess CIB fluctuations have been shown to cross-correlate with the fluctuations in the CXB at the $5\sigma$ level (Cappelluti et al., 2013; Mitchell-Wynne et al., 2016; Cappelluti et al., 2017b; Li et al., 2018). This has prompted some authors to attribute this excess to accreting massive black holes (MBHs), by which we refer to the black holes at the centres of galaxies that power AGN. A later cross-correlation analysis performed by Li et al. (2018) found that the spectral energy distribution (SED) of the fluctuations is consistent with unobscured black hole accretion. On the other hand, broad-band SEDs of these fluctuations in the near infrared (NIR) are consistent with a Rayleigh-Jeans spectrum (Matsumoto et al., 2011). Hence, other authors have attributed the excess to stars in these extra-galactic haloes. Currently, there is consensus that this excess power is real and astrophysical in nature, but its origin remains unsettled and controversial.

Helgason et al. (2014) perform the most comprehensive study to date of the known source populations contributing to both the CIB and CXB. They determine that while extrapolations of known galaxy and AGN populations can explain the coherence of the backgrounds on scales $<1^{\prime}$, only 3 per cent of the measured cross-power between the CIB and CXB on larger angular scales can be obtained from known populations at redshifts $z\leq 5$. Based on their luminosity functions, known populations of AGN are expected to contribute negligibly to the CIB auto-power fluctuations. Theoretical estimates of galaxy populations at yet higher redshifts also fall short of explaining the excess by an order-of-magnitude (Helgason et al., 2012; Yue et al., 2013a).

The two predominant hypotheses for the explaining the proposed excess with new populations are intra-halo light (IHL) (Cooray et al., 2012) and $z>10$ Compton-thick accreting MBHs (Yue et al., 2013b). In the IHL interpretation, there exists uncharacterised starlight beyond the radius at which galaxies are masked. This picture is tantalizing due to the Rayleigh-Jeans-like spectrum of this component in the NIR. On the other hand, the background fluctuations do not appear to be sensitive to the size of the mask (Arendt et al., 2010). More importantly, this model does not provide a mechanism to generate a cross-correlation with the CXB. While X-ray emitting halo gas could hypothetically cross-correlate with the IHL, it would not create enough hard X-rays to be consistent with the amplitude of the fluctuations (Li et al., 2018). If it exists, the IHL component remains to be characterised by either cosmological simulations theoretically or low-surface brightness observations.

Alternatively, an early population of $z>6$ Compton-thick AGN could naturally explain the cross-correlation with the X-ray (Yue et al., 2013b; Cappelluti et al., 2013, 2017b). The Compton-thick nature of these objects is required to prevent them from exceeding the X-ray background and causing reionisation prematurely. However, as we show, in order to produce sufficient integrated intensity to match the observed excess, such a population of accreting MBHs must accrete an amount of matter even greater than the $z=0$ mass density of MBHs (e.g., Yue et al., 2013b) at these early epochs. This is not expected or predicted in any models that track early MBH growth, which in fact generally struggle to grow MBHs rapidly enough to explain $z\sim 6$ quasars without special conditions such as massive seeds, super-Eddington accretion rates, and galactic inflows that are not curtailed by feedback (e.g., Haiman & Loeb, 2001; Volonteri et al., 2015; Pacucci et al., 2017c; Woods et al., 2018). Another hypothesized solution also utilizes early black holes, but indirectly. If dark matter is composed of primordial black holes (Bird et al., 2016), extra Poissonian fluctuations in the early universe would lead to isocurvature density fluctuations that create additional low-mass haloes that could contribute to the background (Kashlinsky, 2016).

Previous work thus far has not employed self-consistent models of MBH assembly tracked over cosmic time that are constrained by the wealth of currently available multi-wavelength data. To decisively settle whether a population of early MBHs could indeed explain the observed excess, we adopt a new comprehensive, self-consistent, and empirically motivated approach that enables the calculation of correlation functions. The standard approach at present, referred to as a semi-analytic model (SAM), has been very successfully used to track the co-evolution of dark matter haloes, galaxies and their MBHs, standing as an established and robust method to interpret galaxy and AGN surveys (e.g., Kauffmann et al., 1993; Volonteri et al., 2003; Somerville et al., 2008; Benson, 2012). Using our comprehensive modelling, we find that a realistic $z>6$ MBH population falls many orders of magnitude short in flux to explain the observed excess. We then explore alternative models: first, a population of buried $z\sim 2-3$ AGN, and then a simple power-law which may be due to a missing local component. We rule out an AGN explanation to the missing fluctuations, but also argue that no other competing physical model can explain all of the observations.

In §2, we describe the data that we use, which we adopt without reanalysis. In §3, we establish our theoretical framework, which consists of a semi-analytic model for MBH growth §3.1, simulated spectra from hydrodynamical simulations §3.2, and a new formalism to describe the clustering of undetected sources §3.3. The results of our calculations are shown in §4, where we show that $z>6$ MBH accretion cannot be responsible for the excess. In §5, we explore alternate explanations: the more moderate-redshift, buried AGN as well as a new Galactic foreground. Our work is summarized in §6.

2 Data

We utilize the most recent and carefully estimated power spectrum measurements of Cappelluti et al. (2017b) without performing any reanalysis. These include measurements of the cross-power spectra of the unresolved Spitzer-IRAC cosmic infrared background and 0.5-2 Chandra-ACIS CXB fluctuations. These include data from the Chandra Deep Field South, Hubble Deep Field North, Extended Groth Strip/AEGIS field, and UDS/SXDF surveys, where a $>5-6\sigma$ cross-power signal is detected on angular scales $>20$” between [0.5-2] keV and the 3.6 and 4.5 $\mu$m bands, as well as 3.6+4.5 $\mu$m combined. No significant correlation was found with harder X-ray bands. Moreover, they find an excess of about an order of magnitude at $5\sigma$ confidence when compared to models for the contribution of the known unmasked source population at $z<7$ from Helgason et al. (2014). Mitchell-Wynne et al. (2016) detected a similar signal in the Chandra Deep Field South. Recently, Li et al. (2018) extended and confirmed these results and detected the signal on even larger scales - up to scales of $\sim$1000”-1500” using the Spitzer-SPLASH survey and the Chandra-COSMOS legacy survey data. The signal even on these scales is still in excess with respect to that produced by foreground sources.

IR maps were self-calibrated and Zodiacal light removed following Arendt et al. (2010). Resolved IR sources were subtracted down to AB magnitude $m\sim 25$ in both the 3.6 and 4.5 $\mu$m IRAC bands (Kashlinsky et al., 2012). 0.5-2 keV, X-ray maps of each individual field were particle subtracted and masked from sources brighter than a few $10^{-17}\ \mathrm{erg}\;\mathrm{cm}^{-2}\;\mathrm{s}^{-1}$. These maps were employed to measure the cross-power spectrum with Fourier analysis. Errors are estimated with a simple cosmic variance estimator (see Cappelluti et al., 2013).

In addition to the cross-power signal, Li et al. (2018) were able to determine that the X-ray spectrum of the CXB fluctuations is coherent with those in the CIB and found it can be fitted by either a power-law component with spectral index $\gamma=2$ or with a $z\sim 7-15$ population of absorbed Compton-Thick AGN or with a direct collapse black hole (DCBH) template from Pacucci et al. (2015).

3 Methodology

We develop a comprehensive theoretical framework to predict the MBH contribution to the power spectra of the CIB and CXB. We begin with a semi-analytic model for MBH growth from the seeding epoch to the present day (Ricarte & Natarajan, 2018a, b). For each individual MBH in the model, we compute its contribution to the CIB and CXB using hydrodynamical simulations to extract the multi-wavelength spectra of early MBHs. We begin with the standard halo occupation distribution (HOD) methodology, and derive a new formalism to determine the power spectrum of sources that are not individually detected. Our framework and approach are illustrated graphically in Fig. 1. Below, we explain each aspect of this modelling in more detail. Note that in this work we do not consider a hypothetical primordial black hole population as in Kashlinsky (2016).

3.1 Semi-analytic Modelling

We begin with a semi-analytic model for the assembly of MBHs that starts at the seeding epoch. We explore both the “light” and “heavy” seeding models of Ricarte & Natarajan (2018a, b). These models integrate the evolution of MBHs from $z=20$ to $z=0$ using simple prescriptions that relate the evolution of MBHs to their host galaxies. They have been tested against the local $M_{\bullet}-\sigma$ relation (relating MBH mass to host velocity dispersion Kormendy & Ho, 2013; Saglia et al., 2016; van den Bosch, 2016) and bolometric luminosity functions for $0<z<6$ (Hopkins et al., 2007). A hybrid approach is taken whereby empirical relations are used for the galaxy-halo connection, while analytic prescriptions are used to model the MBH-galaxy connection.

First, dark matter merger trees are constructed for a statistically representative sample of the universe, using an extended Press-Schechter algorithm (Press & Schechter, 1974) calibrated to the Millennium simulations (Parkinson et al., 2008). We first estimate the contribution of the undetected high redshift population of growing black holes in our the context of our model. Since we are interested only in the high-redshift universe, the parent nodes of these merger trees originate at $z=4$. We sample 35 halo masses between $10^{10-13.4}\ \mathrm{M}_{\odot}$, 100 times each in order to capture cosmic variance. This wide range of halo masses includes even highly biased $5\sigma$ peaks which are implicated as the sources of the most distant quasars (Fan et al., 2003).

With these trees in hand, we then integrate the assembly of MBHs from $z=20$ to $z=4$. Black hole seeding is permitted in the redshift range $15<z<20$, when the pristine nature of gas promotes the formation of objects with large Jeans masses (Bromm et al., 2002; Agarwal & Khochfar, 2015). In light seed models, any halo which represents a 3.5$\sigma$ peak or higher is randomly assigned a 30-100 $\mathrm{M}_{\odot}$ seed (Stacy et al., 2016). In heavy seed models, direct collapse black hole seeds in the mass range $10^{4-6}\ \mathrm{M}_{\odot}$ form in haloes with low spins and virial temperatures appropriate for atomic-cooling haloes (Lodato & Natarajan, 2006). A detailed description of our seeding model can be found in (Ricarte & Natarajan, 2018b). The origin of seed black holes is currently an open question and observations suggest that more than one channel to produce early seeds is likely to operate in the universe (Volonteri & Bellovary, 2012). Therefore, investigating potential differences between light and heavy seed populations that might enable discrimination between them is warranted (Natarajan, 2014; Pacucci et al., 2018; Natarajan et al., 2019). We note here that there are no numerical simulations currently available that can self-consistently track seed formation and growth over cosmologically relevant volumes and then compute measurable quantities like clustering statistics - therefore, we rely currently on semi-analytic modelling.

As in previous work with this fiducial model, major mergers trigger MBH growth until MBHs reach the the $M_{\bullet}-\sigma$ relation, where $\sigma$ represents the stellar velocity dispersion within the effective radius of the host galaxy. In order to directly compare with the observational data, it is important that $\sigma$ not be the velocity dispersion of the dark matter halo, which evolves differently with redshift. We instead estimate $\sigma$ based on a combination of the stellar-to-halo mass relation (Moster et al., 2013) and the stellar mass-to-size relation (Mosleh et al., 2013; Huertas-Company et al., 2013). In previous work, we had set the black hole growth rate to the Eddington rate during this period. In this work, we instead draw Eddington ratios randomly from the observed distribution of observed broad-line quasars (Kelly & Shen, 2013; Tucci & Volonteri, 2017). These distributions allow for some super-Eddington accretion at high-redshift, but again, only until MBHs reach the $M_{\bullet}-\sigma$ relation. In Appendix A, we show that these accretion prescriptions allow us to match luminosity functions at $0<z<6$ (for a separate set of low-redshift merger trees) as well as obtain reasonable values for the clustering of AGN at low-redshift that are in agreement with observational data.

Once these calculations are complete, we have a large population of haloes with known AGN bolometric luminosities. These bolometric luminosities are self-consistently related to MBH accretion rates via $L_{\mathrm{bol}}=\epsilon_{r}\dot{M}_{\bullet}c^{2}$, where $\epsilon_{r}$ is the radiative efficiency (set as usual to 0.1) and $c$ is the speed of light. To proceed, we now need to determine what fractions of these bolometric luminosities are emitted in the relevant bands, and additionally we develop a formalism for calculating the clustering of undetected sources.

3.2 From MBH Growth Rates to Fluxes

3.2.1 The Spectra of High-z Black Hole Seeds

The contribution of MBHs to the CXB/CIB depends on the underlying spectral template assumed for the output emission of these accreting sources. Rather than simply treating these bolometric corrections as free parameters, we employ the radiation-hydrodynamic code GEMS (Growth of Earth Massive Seeds Pacucci et al., 2015, 2017a) for this task. This code was specifically designed to simulate the accretion and irradiation processes on high-redshift black hole seeds.

GEMS simulates the accretion flow on to high-redshift black hole seeds, focusing on spatial scales of the order of the Bondi radius $R_{B}=2GM_{\bullet}/c_{s}^{2}(\infty)$ (Bondi, 1952) of the central source, where $M_{\bullet}$ is the black hole mass and $c_{s}(\infty)$ is the sound speed of the unperturbed gas. While the code does not simulate the formation of the seed itself, once the initial mass $M_{\bullet}(t=t_{0})$ is set GEMS is able to compute the output spectrum until complete gas depletion (Pacucci et al., 2017b). In particular, GEMS computes the frequency-integrated radiative transfer through the gas content of the host. The frequency-dependent radiative transfer is performed in a post-processing step with the publicly available spectral synthesis code CLOUDY (Ferland et al., 1998; Ferland et al., 2013). The final result is the calculation of the time-evolving spectrum emerging from the galaxy hosting the black hole seed. The time-dependent data used by CLOUDY includes the following: (i) spatial profiles for hydrogen number density and gas temperature, (ii) the shape of the emission spectrum of the central source of radiation, and (iii) the bolometric luminosity $L_{\mathrm{bol}}$ of the source. The latter is computed directly from the accretion rate. The SED of the source, ranging from far-infrared to hard X-ray, includes the three, classic, components: (i) a multicolour blackbody, (ii) a power-law, and (iii) a reflected component (see e.g. Yue et al., 2013b, for details).

In order to feed initial conditions into the semi-analytic models, the spectral templates for black hole seeds were computed for seeds of 5 different masses ranging from $10^{2}\ \mathrm{M}_{\odot}$ $10^{6}\ \mathrm{M}_{\odot}$, spanning the Pop III remnant to DCBH formation channels (Pacucci et al., 2017b; Woods et al., 2018). The templates are computed at different redshifts, to input into the semi-analytic models with the correct spectral properties depending on the cosmic time. With a redshift step $\Delta z=2$, the range of variation spans from $z=6$ (end of the epoch of reionization Kashikawa et al., 2006) to $z=20$ (formation of the first stars Miralda-Escudé, 2003). We also explore the effect of the hydrogen column density $N_{H}$ on the spectral energy distribution (Pacucci et al., 2017a), and evaluate the effect of values of $N_{H}$ between $10^{23}\ \mathrm{cm}^{-2}$ (i.e., Compton thin) and $10^{26}\ \mathrm{cm}^{-2}$ (i.e., extremely Compton thick) for these early accreting sources. We then select the simulation snapshot which maximizes the total output flux. An example of a typical spectral template for a $10^{5}\ \mathrm{M}_{\odot}$ black hole seed is shown in Fig. 2. The template extends from the far infrared to the hard X-ray bands.

From these spectral templates, we derive the bolometric corrections as a function of the redshift and mass of the accreting black hole. These bolometric corrections are needed in order to compute the variation $dF/dz$ of the flux $F$ in the infrared and X-ray bands as a function of the redshift. The bolometric correction is defined as the ratio between an object’s luminosity emitted in some band to its bolometric (total) luminosity. Fig. 3 shows the bolometric corrections for a fixed ($10^{5}\ \mathrm{M}_{\odot}$) black hole mass but at different energies with the 1-$\sigma$ uncertainty bands. Overall, the bolometric corrections for infrared bands increase with redshift, with a mean value of $\sim$2 per cent at $z=6$ and reaching $\sim 7$ per cent at $z=20$. The relevance of infrared bands for these sources increases with increasing redshift, confirming the importance of the upcoming James Webb Space Telescope (JWST) mission to search for high-redshift black hole seeds (Natarajan et al., 2017). The bolometric corrections for X-ray bands are higher, with an average value of $\sim 12$ per cent for both bands. The soft (0.5-2.0 keV) and hard (2.0-10.0 keV) X-ray bands have opposing trends with redshift, due to the fact that the bell-shaped X-ray spectrum shifts with redshift. At $z\sim 15$, the bolometric corrections for soft and hard X-ray bands equal each other at $\sim 12$ per cent of the total luminosity of the source.

3.2.2 Masking Sources

In the observations, resolved sources are masked in order to sift out the diffuse component of the cosmic backgrounds. We include the effect that this selection has on the power spectra in our model, removing any accreting MBH with an IR flux brighter than an AB magnitude $\sim 25$ or an X-ray flux brighter than $\sim 10^{-17}\ \mathrm{erg}\;\mathrm{s}^{-1}\;\mathrm{cm}^{-2}$. The infrared component of the masking is a simple cut on the apparent magnitude at 4.5 $\mu m$, while in fact the X-ray component has a slightly more complicated form due to Chandra’s variable sensitivity as a function of off-axis angle and epoch. Encapsulating these effects, we define the unmasked probability, $\eta(F_{x},F_{IR})$, to be

[TABLE]

where $\zeta(F_{x})$ is the probability that a source would be undetected in Chandra’s AEGIS field. $\zeta(F_{x})$ is reproduced in Fig. 4. We have experimented with an IR magnitude cut of 24 instead of 25, and found that this increases the total flux of the derived background by a factor of $\sim 2$. As we shall show, our main results are insensitive to uncertainties of this level.

3.3 The Power Spectra of Undetected Sources

We begin with the methodology outlined in (Helgason et al., 2014) in order to compute the clustering signal of high-redshift AGN, but make new and critical modifications in order to derive the power spectrum of undetected sources. Traditional halo occupation distribution (HOD) modelling is based on the number of central and satellite sources as a function of halo mass (Cooray & Sheth, 2002). This implicitly assumes that all accreting MBHs in the model have the same flux (at least at fixed redshift), which is generally untrue. A generic prediction of almost all MBH assembly models is that more massive haloes host more massive MBHs, and therefore more luminous MBHs cluster more strongly.

Instead of computing the number of central and satellite sources as a function of halo mass, we compute the flux attributed to central and satellite galaxies directly from the SAM. Reformulating in this way avoids issues defining an AGN based on its Eddington ratio or luminosity and allows us to correctly treat AGN at different redshifts. For example, if there is a fixed flux below which sources are part of the background, the power spectrum in that event probes more luminous sources as redshift increases.

The observed power spectrum can be expressed as a sum of the clustering and shot noise of the emitting population,

[TABLE]

where $q$ is the angular wavenumber in rad*-1*. Note that the shot noise component is a scale-independent floor. The first term can be computed by integrating contributions as a function of redshift,

[TABLE]

where $d_{c}(z)$ is the comoving distance to redshift $z$, $H(z)=H_{0}\sqrt{\Omega_{M}(1+z)^{3}+\Omega_{\Lambda}}$, and $dF_{i,j}/dz$ is the flux gradient of the unmasked sources in observed bands $i,j$. To proceed, we therefore compute the contribution to the angular power spectrum at fixed redshift intervals, then integrate these pieces.

3.3.1 Total Flux Per Redshift Interval

In this section, we compute the total flux in each relevant band from MBHs as a function of redshift, $dF_{i,j}/dz$ in Equation 3. The flux gradient in band $i$ is given by

[TABLE]

We can express $d^{3}N/dF_{i}d\Omega dz$ as

[TABLE]

Here, the first factor is the volume-weighting, which comes from standard cosmology,

[TABLE]

The second factor is the flux distribution, calculated directly from the semi-analytic model after applying bolometric corrections. The results of these calculations are shown in Fig. 5 for the soft X-rays, 3.6 $\mu m$, and 4.6 $\mu m$. Note that these calculations are performed at the observed, not intrinsic, wavelength. The light seed model produces a factor of 2.0 more background, due to the increased abundance of these seeds compared to the rarer heavy ones.

3.3.2 Flux-weighted HOD Modelling

In this section, we compute how the flux at a given redshift is spatially distributed, $P_{i,j}(qd_{c}^{-1},z)$ in Equation 3. In this section, all quantities are computed only for a single redshift slice and for a single band. We employ the halo occupation distribution (HOD) modelling (Cooray & Sheth, 2002) technique, but modify from previous approaches by weighting our accreting MBHs by flux. In this picture, cosmology provides the clustering of dark matter haloes, and the distribution of sources within these haloes modifies the dark matter clustering signal.

HOD modelling typically begins with the number of centrals and satellites as a function of halo mass as its basis. Again, this is not appropriate for the clustering of an unresolved background, because this implicitly assumes that AGN luminosities are independent of halo mass. Instead, we directly compute from our model the total amount of flux associated with centrals and satellites in a given band.

The power spectrum can be written as a sum of the one- and two-halo terms:

[TABLE]

These represent components arising from the intra- and inter-halo clustering respectively. In order to solve for these components, we first determine the average AGN flux $F$ in the band that is contributed by the central and satellite populations, weighted by the unmasked probability:

[TABLE]

and

[TABLE]

Since we use Monte Carlo merger trees in order to generate our population, the distinction between central and satellite is simple: satellites are those sub-haloes which have merged in the merger tree, but have not experienced enough time orbiting their central halo to completely disrupt (the dynamical friction time-scale, Boylan-Kolchin et al., 2008). We assume that the light of the satellite population is on average distributed in a Navarro-Frenk-White (NFW) profile (Navarro et al., 1997). In practice, the one-halo term is negligible for our study due to the high redshifts we consider, but we include it for completeness in our formalism. The mean flux density at a given redshift can then be calculated by integrating both components over the halo mass function (Murray et al., 2013).

[TABLE]

The flux-weighted central and satellite biases tell us how the distribution of sources is related to the distribution of dark matter at this epoch. This can then be computed via,

[TABLE]

and

[TABLE]

where $b(M_{h}z)$ is the linear halo bias using the ellipsoidal collapse formalism (Sheth et al., 2001), and $u(k|M_{h},z)$ is the Fourier transform of an NFW profile, normalized such that $u(0)=1$ (Scoccimarro et al., 2001). This depends on the mass-concentration relation of dark matter haloes, for which we adopt the $z=5$ relation from (Dutton & Macciò, 2014). Finally, having constructed these quantities for some bands $i$ and $j$, we can construct the one-halo component,

[TABLE]

which contains terms correlating centrals and satellites as well as satellites with themselves, and the two-halo component,

[TABLE]

where $P^{\mathrm{lin}}(k,z)$ is the linear $\Lambda$CDM power spectrum (Murray et al., 2013).

While this quantity is not used, it is instructive to compute the “effective” bias via

[TABLE]

We plot the effective bias for both heavy and light seeds as a function of redshift in Fig. 6. Note that the effective bias of the background at $z=6$ ($\approx 4$) is less extreme than that of quasars at $z\sim 3$ ($\approx 6$ e.g., Timlin et al., 2018). This is because the typical undetected source resides in a much lower mass dark matter halo than a quasar. Also note that the bias of heavy seeds is higher than that of light seeds, since heavy seeds preferentially occupy higher mass haloes, but this effect is slight.

3.3.3 Shot Noise

Computing the shot noise floor is the final missing piece of the power spectrum, $P_{\mathrm{SN}}$ in Equation 2. The shot noise power between any given bands $i$ and $j$ is given by

[TABLE]

where $d^{2}N/dF_{i}dF_{j}$ is the distribution of sources with observed fluxes $F_{i}$ and $F_{j}$ per unit solid angle. In practice, we compute $P_{SN}$ at each redshift slice and then integrate over redshift, weighting each component by $d^{2}V/dzd\Omega$.

3.3.4 One Final Modification: Beam Smearing

We have now developed the tools to calculate the intrinsic angular power spectrum of hidden/undetected sources on the sky. Yet in reality, the finite size of the instrument beam suppresses power on small scales. The final step to computing the observed power spectra is to convolve the signal with the instrument’s point spread function (PSF). Since we are working in Fourier space, we can simply multiply by its Fourier transform. We approximate the PSFs of Spitzer and Chandra as Gaussians with full-width half-maximum values of 1.2 arcseconds and 4 arcseconds respectively (the latter being the average over the field of view). Deviations from Gaussianity are not important for the large scales that we consider in this work.

We do not convolve our model with the mask. As shown in Kashlinsky et al. (2012), mask convolution is only important on small scales, and we exclude $<30$ arcsecond measurements from our analysis.

4 Power Spectra of Unmasked $z>6$ Black Holes

We now present the power spectra of unmasked $z>6$ MBHs as predicted by this self-consistent model of MBH assembly, bolometric corrections adopted from hydrodynamical simulations, and a new formalism for computing the power spectrum of undetected sources. The X-ray auto-power spectrum is shown in Fig. 7, while the IR auto-power spectrum and X-ray to IR cross power spectrum are shown in Fig. 8. The power spectrum from the heavy seeding model is shown in red, while that of the light seeding model is shown in blue. Known foregrounds are shown with a black line (Helgason et al., 2014), while the data from Cappelluti et al. (2017a) are shown as black points.

It is immediately clear from the IR auto-power spectrum (Fig. 8) that growing $z>6$ MBHs, assembling according to our current understanding of MBH evolution, fall many orders of magnitude short of explaining the excess. There is very little difference between the power of light and heavy seeds. Although we found in Fig. 6 that heavy seeds cluster slightly more strongly, the light seeds actually contribute more power. This is due to two factors: (i) light seeds populate more haloes than heavy seeds do, by construction, so there are simply more MBHs that are able to contribute to the background, and (ii) since light seeds are initialized with lower masses than heavy seeds, they must accrete luminously to make up the difference. The end result is that light seeds produce a factor of 1.7 greater total intensity, and very slightly overtake the heavy seeds.

As seen clearly, this comprehensive model fails dramatically in explaining the excess. We consider a number of possible explanations, and determine that constraints at lower redshifts on the total mass density of MBHs prevent our model from overgrowing high-redshift MBHs.

4.1 Higher Bolometric Corrections?

Is it possible for us to match the data if a larger fraction of the total flux is attributed to the relevant bands? The answer is no, because each IR band already contains $\approx 2$-$10$ per cent of the total flux. In the CIB power spectrum, the deficiency of 4 orders of magnitude that we predict corresponds to a deficiency of 2 orders of magnitude in total flux. However, by adopting a more optimistic bolometric correction, the total flux can only possibly be increased from 2-10 per cent to $\approx 50$ per cent in each IRAC band.

Previous work has invoked an unusually high obscuring column of $\sim 10^{25}\ \mathrm{cm}^{-2}$ in order to produce the observed ratio between the CIB and CXB fluctuations and prevent premature reionization of the universe (Yue et al., 2013b). Our bolometric corrections do not provide as extreme IR-to-X-ray flux ratios because they are taken from the simulation snapshot when the flux is maximized, which is more representative of the time-averaged SED. Even if simulations are initialized with extremely Compton-thick column densities, the maximum flux is not emitted until the obscuring column has diminished (Pacucci et al., 2015). One might then invoke a continuous renewal of the extreme column density, by means of some gas inflow from the outside environment. In this case, the reduced radiative efficiency would require proportionally more MBH accretion, increasing the tension with the MBH mass density which we explore in subsequent sections.

4.2 Larger Halo Biases?

Would it be possible to match the data if a larger bias is assumed? This is impossible due to the magnitude of the discrepancy. Since the power spectrum is proportional to the square of the bias, the bias would need to be increased by 2 orders of magnitude. While biases of flux-limited AGN samples tend to be a factor of $\approx 2$ higher than what we predict in Fig. 6 (Cappelluti et al., 2012), haloes with biases 2 orders of magnitude higher do not exist in the observable universe.

Note that we do not treat the halo bias as a free parameter, but calculate it directly based on the haloes which host MBHs in our model, which is constrained by the local $M_{\bullet}-\sigma$ relation and redshift-evolving luminosity functions. The effective bias reflects the haloes whose MBHs contribute most to the background, excluding ones with MBHs that are luminous enough to be masked.

4.3 Greater MBH Occupation Fractions?

Is it possible for us to match the data by increasing the total number of MBHs? We argue that it is not, because our model is constrained by the MBH occupation fraction of low-mass galaxies at $z=0$ (Miller et al., 2015). The ability for $z=0$ measurements to constrain $z>6$ predictions underscores the power of our self-consistent modelling. As shown in Ricarte & Natarajan (2018b), our heavy seeding prescription is consistent with current constraints based on X-ray emission, while the light seeding prescription optimistically assigns almost every low-mass galaxy a MBH seed.

The only way to increase the number of MBHs that contribute to the background is to assign additional MBHs to galaxies which already have them. While we do expect galaxies to host “wandering” MBHs (Tremmel et al., 2018), they could not exist in sufficient quantities to explain this deficit. The deficit of 2 orders of magnitude in total flux requires that every galaxy by $z=6$ contain 100 times as much mass in wandering MBHs as at the centre of galaxy. In our model, massive galaxies at $z=6$ already acquire an optimistic black hole-to-total stellar mass ratio of $\approx 1:100$. This means that every galaxy would need to accrete as much mass onto MBHs as exists in stars, and most of this would be off-centre. As we further discuss in §4.4, even if this enormous excess mass in wandering MBHs exists, it must somehow be prevented from merging with the central MBH and becoming part of the $z<6$ MBH mass budget. It is difficult to simultaneously argue that 100 times more mass is accreted onto MBHs, but these masses are also excluded from the centres of galaxies.

4.4 More Black Hole Accretion?

Is it possible for us to match the data by simply allowing for greater amounts of MBH accretion? Here we argue that constraints at lower redshifts do not permit this, and this is the main source of our disagreement with previous studies. Recall that MBH growth is capped by the $M_{\bullet}-\sigma$ relation in our model. It is this cap which limits the rate at which MBHs are allowed to grow on cosmic timescales, not the Eddington ratio distribution, which permits super-Eddington accretion rates in our model. The $M_{\bullet}-\sigma$ mass cap regulates the shape of the luminosity function at each epoch (see Appendix A), ensures that appropriate masses are obtained as a function of host mass at $z=0$, and determines the black hole mass density in place at $z=6$.

The total amount of mass locked up in MBHs in our model is shown in Fig. 9. We include a large selection of observational estimates based on MBH-host scaling relations at $z=0$, and accretion at higher redshift (Shankar et al., 2009; Hopkins et al., 2007; Treister et al., 2009; Salvaterra et al., 2012; Treister et al., 2013; Cappelluti et al., 2017a). The large range of upper limits at $z\sim 6$ is due to different assumptions regarding the fraction of accreted black hole mass that is converted into observable X-ray emission. The constraints from Shankar et al. (2009); Hopkins et al. (2007) should be considered more robust, however, as these are based on local empirical relationships between MBHs and their hosts and bolometric luminosity functions respectively. We plot the total amount of mass locked up in MBHs with dashed lines (this includes seed masses), and the accreted portion of the mass locked up on MBHs as a solid region. In green, we also plot the total integrated mass density required by the Yue et al. (2013b) model which explains the missing fluctuations with early MBH growth. We also provide these mass densities at $z=6$ in Table 1. Only $10^{3-4}\ \mathrm{M}_{\odot}\;\mathrm{Mpc}^{-3}$ is in place at this epoch in our model, while the Yue et al. (2013b) interpretation requires $>10^{6}\ \mathrm{M}_{\odot}\;\mathrm{Mpc}^{-3}$ to be in place by $z=10$.

In the Yue et al. (2013b) model, so much mass is accreted onto MBHs that the integrated mass density at $z=10$ already exceeds the mass density at $z=0$. Such large mass densities are required simply to produce enough flux, as explored in detail in Helgason et al. (2016). As discussed in §4.3, this leaves no room for MBHs to grow at lower redshifts. These accreting MBHs would somehow need to simultaneously accrete far more efficiently than expected, and yet remain displaced from the centres of galaxies throughout all of cosmic time. Note that burying these sources in large columns does permit such a model to agree with the intensity of the X-ray background (Cappelluti et al., 2017a)—however, it is bolometric luminosity density measurements at lower redshift (Hopkins et al., 2007) and local scaling relations (Shankar et al., 2009) which rule out this interpretation.

5 Discussion

The comprehensive model developed here fails by many orders of magnitude to reproduce the missing CIB/CXB fluctuations with $z>6$ accretion. A successful model would need to produce about 2 orders of magnitude more total flux in the IR. The main reason that our results disagree with previous attempts to explain this excess with accreting black holes (Yue et al., 2013b) is that this previous work requires more mass locked up in MBHs at $z=10$ than there exists at $z=0$, which is not permitted in our models of self-consistent growth over cosmic time. In this section, we explore potential alternative models. In §5.1, we consider missing AGN at lower redshift. In §5.2, we consider a missing Galactic foreground component.

5.1 Can we be missing buried moderate-redshift AGN?

We now explore the existence of missing AGN not at $z>6$, but rather at $z\sim 2-3$, unaccounted for in our census of SMBH growth. Missing SMBH accretion can be more easily accommodated at these epochs than at $z>6$ due to the amount of SMBH mass already in place. Interestingly, population synthesis models that reproduce the observed spectral shape of the CXB require the existence of a population of Compton-thick actively accreting black holes that may account for 30-50 per cent of SMBH growth (Treister et al., 2004; Ballantyne et al., 2006; Gilli et al., 2007; Aird et al., 2015; Ananna et al., 2019). Some evolutionary models of AGN triggering suggest that MBH growth may begin in a highly obscured phase before blowing away the obscuring column and becoming observable as a traditional AGN (Hopkins et al., 2008; Alexander & Hickox, 2012; Blecha et al., 2018). Such an evolutionary phase before each AGN episode could be accommodated in a self-consistent model by modifying the average radiative efficiency of an accretion disk (Soltan, 1982).

We hypothesize that the flux from such missing buried sources could be distributed as a function of redshift according to the luminosity densities of either obscured AGN (Aird et al., 2015) or red type I quasars (Lacy et al., 2015), whose redshift distributions are shown in Fig. 10. Since these redshift distributions peak around $z\sim 2$ and $z\sim 3$, we label these hypothetical missing populations AGN_z2 and AGN_z3 respectively. We further assume that they occupy a characteristic halo mass of $10^{13}\ \mathrm{M}_{\odot}$, the typical halo mass of X-ray selected AGNs (Cappelluti et al., 2012). This characteristic halo mass provides an optimistically high bias factor, since such high mass haloes are likely to host galaxies which would be masked. We do not include a satellite contribution to the flux of this population, and set the one-halo term to zero.

Nevertheless, we use these optimistic assumptions to simultaneously fit all power spectra for 4 free parameters: the total amount of isotropic flux in the IRAC1 band, flux ratios to the IRAC2 and soft X-ray bands, and the coherence between the CIB and CXB. Best fit values are provided in Table 2, with errors are estimated using the Markov Chain Monte-Carlo package emcee (Foreman-Mackey et al., 2013). Values for $\chi^{2}$, reduced $\chi^{2}$, and the Bayesian Information Criterion (BIC) are provided in Table 3. The BIC is a useful metric for comparing models with different numbers of parameters (Schwarz, 1978). Additional details of this optimization process can be found in Appendix B.

Power spectra for the best fit hypothetical populations are shown in Figs. 11 and 12. We find that these generic extra-galactic populations also provide a poor fit to the data, with $\chi^{2}/DOF\sim 20$. The incompatibility of this model with the data is evident in Fig. 12, in the computed auto-power spectrum of the CIB. The two-halo term of dark matter clustering provides too much power on degree scales relative to arcminute scales, and the best fit attempts to compensate by undershooting the data on small scales and overshooting on large scales.

In addition, as with Yue et al. (2013b), the best-fit isotropic flux value translates to an amount of mass in tension with the local black hole mass density. Following the “Soltan argument” (Soltan, 1982), allowing 10 per cent of the bolometric flux to be radiated in each of the IR bands, and assuming a typical radiative efficiency of 0.1, we find that this extra flux would correspond to an extra black hole mass density of $10^{6.0}\ \mathrm{M}_{\odot}\;\mathrm{Mpc}^{-3}$ for either model, once again in excess of current estimates of the $z=0$ total mass density (Shankar et al., 2009). (For more details, see Appendix C.) While more informed priors in the fitting process could reduce this value, this would most likely come at the expense of even worse fits to the data.

In principle, the underproduction of power at sub-arcminute scales could be alleviated with the addition of a one-halo term to the power spectrum. We therefore experimented with a model in which all of the flux is attributed to satellites distributed according to the appropriate NFW profile for a given mass and redshift. We found that this modification could boost the power on 10 arcsecond scales by up to 50 per cent for standard halo concentrations. This increases to up to a factor of 4 if a much higher concentration parameter of 50 is assumed, which some authors suggest is a value of the NFW concentration parameter required to match small-scale AGN clustering (Hennawi et al., 2006; Lin & Mohr, 2007; Kayo & Oguri, 2012; Eftekharzadeh et al., 2019). However, even in the latter case, the one-halo term is buried beneath the foreground of known sources at these scales.

In addition, the spectral properties of the fluctuations also do not favor an obscured AGN scenario at moderate redshift. Even obscured AGN do not typically emit strongly at 1 micron, which in fact is a good band to probe the stellar continuum (e.g., Collinson et al., 2017; Hickox & Alexander, 2018). Reprocessed emission usually emerges in the far infrared, where the optical depth falls. In addition, a recent analysis of the spectrum of the fluctuations favors unobscured AGN over an obscured population (Li et al., 2018). This surprising result is not predicted by any current explanations of the CIB excess.

5.2 Are we missing a Galactic foreground?

The power spectrum of the fluctuations resembles a simple power law much more than the shape of the extra-galactic two-halo term. In addition, the spectra resemble a Rayleigh-Jeans tail in the NIR and unobscured black hole accretion in the X-ray. Guided by these facts, we explore the hypothesis of a new Galactic foreground of warm gas and dust spatially distributed separately from colder Galactic cirrus. We show that although the NIR data would prefer a simple power law that may be due to a Galactic component, reflected Galactic X-ray light is insufficient to explain the coherence with the CXB.

Observationally, it has been difficult to estimate the contribution from Galactic cirrus to the NIR fluctuations, as discussed in detail in Kashlinsky et al. (2018). The main method with which to assess its contribution has been to attempt to find cross-correlations with bands where cirrus is known to be dominant. Using Spitzer data, Kashlinsky et al. (2012) find marginal cross-correlations with the NIR and 8 micron emission, while Matsumoto et al. (2011) find no cross-correlation with longer wavelength 100 micron maps with AKARI. In contrast, Thacker et al. (2015) do report a cross-correlation at longer 250, 350, and 500 micron wavelengths, but Kashlinsky et al. (2018) argues that the claimed signal may be derived from known Galactic populations. At the shorter wavelengths probed by HST, Mitchell-Wynne et al. (2015) observe a clear power-law on scales $>1^{\prime}$, which they attribute to diffuse Galactic light. Zemcov et al. (2014) report a cross-correlation between 1.1 and 1.6 micron fluctuations with the 3.6 micron fluctuations, which they attribute to IHL. However, Yue et al. (2016) argue that this cross-correlation is due to diffuse Galactic light based on its power-law behaviour reported by Mitchell-Wynne et al. (2015).

To explain the coherence with the CXB, we hypothesize that this foreground reflects Galactic X-ray emission, as speculated by Cappelluti et al. (2017b). We show, however, that such reflected emission is too faint for to plausibly explain the cross-correlation with the CXB. We first perform a simple power-law fit to the missing source of the fluctuations, $P(q)\propto q^{a}$. The free parameters of this model are P1000, which is the value of the power spectrum at $q=1000\ \mathrm{rad}^{-1}$ (or $0.36^{\circ}$), the slope $\alpha$, the flux ratios between different bands, and the coherence between the CIB and CXB. Best fit values are shown in Table 4.

The best fit model for a simple power law fit to the CIB and CXB has a power of $5.6\times 10^{-25}\ \mathrm{erg}^{2}\;\mathrm{s}^{-2}\;\mathrm{cm}^{-4}\;\mathrm{ster}^{-1}$ in the soft X-rays at $q=1000\ \mathrm{rad}^{-1}$. The root-mean-square flux of fluctuations at this scale is then equal to $(q^{2}P(q)/2\pi)^{1/2}$ (Kashlinsky & Odenwald, 2000), which yields a value of $3.0\times 10^{-10}\ \mathrm{erg}\;\mathrm{s}\;\mathrm{cm}^{-2}\;\mathrm{ster}^{-1}$. We equate this to the minimum amount of X-ray flux required by our hypothetical Galactic foreground. It is possible that only 20 per cent of this value is required, which is the measured coherence between the CIB and CXB.

We approximate the geometry of the Milky Way’s X-ray emitting population to be an infinite disk, above which lies a scattering screen of column depth $N_{H}\sim 10^{20}\ \mathrm{cm}^{-2}$, as is typical for the observational data derived from the four fields considered in this work. Grimm et al. (2002) determine that the total 2-10 keV luminosity of X-ray binaries in the Milky Way is approximately $2-3\times 10^{39}\ \mathrm{erg}\;\mathrm{s}^{-1}$. For an order-of-magnitude estimate, we assume that the 0.5-2.0 keV flux is the same and distribute these sources uniformly in a disk of radius 15 kpc, allowing us to estimate the soft X-ray flux of the Milky Way as follows:

[TABLE]

The fraction of this flux which encounters the scattering screen and is reflected back can be calculated on the basis of its optical depth, $\tau=N_{H}\sigma_{T}$, where $\sigma_{T}$ is the Thomson cross-section. This yields

[TABLE]

which for our estimated value of $F_{\mathrm{MW}}$ and $N_{H}=10^{20}\ \mathrm{cm}^{-2}$ results in

[TABLE]

This falls two orders of magnitude short of producing the minimum amount of flux required, which would need a column of at least $\sim 10^{22}\ \mathrm{cm^{-2}}$. Consequently, reflected Galactic X-ray emission is not a plausible mechanism for generating these coherent CXB fluctuations.

6 Conclusion

The CIB exhibits unexplained fluctuations on arcminute scales and higher that cross-correlate with CXB fluctuations. This has motivated some authors to explain the origin of these fluctuations as arising from early, obscured MBH growth (Yue et al., 2013b; Yue et al., 2016). For the first time, we examine this hypothesis using a self-consistent semi-analytic model of black hole populations accreting over cosmic time that are capable of incorporating observed constraints at lower redshifts, namely local MBH-to-host relations and bolometric luminosity functions. We combine these calculations with bolometric corrections from hydrodynamical simulations. By deriving a new formalism for the power spectrum of undetected sources, we then compute the power spectrum expected due to $z>6$ MBH growth.

We find that $z>6$ MBH accretion cannot explain the missing fluctuations due to the enormous amount of growth that would be acquired at such early epochs. Previous attempts to explain these fluctuations using early MBH accretion require more mass locked up in MBHs at $z=10$ than there exists at $z=0$.

We then examine alternate solutions to explain the observed excess. First, we hypothesize that buried AGN at more moderate redshifts, missing from our current census, are responsible. We find that even an optimized model for an AGN origin at low-redshift produces a poor fit on arcminute scales while simultaneously overproducing the local MBH mass density. Then, we explore the hypothesis that a new foreground of Galactic dust might produce the missing IR fluctuations by emission and the coherent X-ray fluctuations by reflection. We find that the Galaxy produces too little X-ray flux for this to plausibly explain the coherence between the two backgrounds.

We decisively rule out accreting MBHs as the sources of the missing background fluctuations. However, additional studies will be required to determine its true nature.

7 acknowledgements

We thank Alexander Kashlinsky for pointing out an error in the previous version of this manuscript. We are grateful for our fruitful conversations with Michael Tremmel, Nir Mandelker, and Vivienne Baldasarre. AR acknowledges funding from the NASA Earth and Space Science Fellowship (NESSF) grant number 80NSSC17K0459. FP acknowledges funding from the NASA Chandra award No. AR8-19021A and from the Nederlandse Onderzoekschool Voor Astronomie (NOVA). NC and FP acknowledge NASA support through ADAP grant NNX16AF29G.

Appendix A Model Validation

Our model reproduces bolometric luminosity functions for $0<z<6$, as well as clustering of AGN to $z\sim 3$. Recall that the set of merger trees used in the main text is calculated only for $z\geq 4$. The following is based on a separate set of merger trees, which allow us to follow MBH growth to $z=0$ using parent halos ranging from $10^{10.6-15}\ \mathrm{M}_{\odot}$. See Ricarte & Natarajan (2018a, b) for more details on this methodology.

Bolometric luminosity functions are displayed in Fig. 13, comparing to estimates from Hopkins et al. (2007) and Ueda et al. (2014) for the both the light and heavy seeds. Our model agrees well for both seeds. While our model overshoots somewhat at $z\geq 3$, here the observations are highly uncertain, especially regarding the contribution of obscured BH growth to the faint end. Since we are using a separate set of merger trees which go to $z=0$, rare massive peaks are missing at high-redshift. An orange dashed line demarcates the luminosity beyond which we estimate half of haloes are missing from the dark matter merger trees. At low redshift, artifacts appear and errors increase due to the decreased number of halos produced by the merger tree algorithm.

AGN clustering is shown in Fig. 14, comparing to a wealth of observations for both X-ray selected AGN and broad-line quasars (Timlin et al., 2018; Allevato et al., 2016; Eftekharzadeh et al., 2015; Allevato et al., 2014; Mountrichas et al., 2013; Koutoulidis et al., 2013; Starikova et al., 2012; Krumpe et al., 2012; Allevato et al., 2011; Cappelluti et al., 2010; Coil et al., 2009; Ross et al., 2009; Hickox et al., 2009; Gilli et al., 2009; da Ângela et al., 2008; Francke et al., 2008; Porciani & Norberg, 2006; Yang et al., 2006; Gilli et al., 2005; Croom et al., 2004). In red and orange, we show the clustering predicted by our model for AGN shining above $10^{46}$ and $10^{42}$ $\mathrm{erg}\;\mathrm{s}^{-1}$ bolometric luminosity thresholds. These curves are computed following a traditional HOD modeling technique whereby the number of sources above these luminosity cuts are counted. While an apples-to-apples comparison accounting for all of the selection effects that may affect AGN clustering measurements is beyond the scope of this paper, it is encouraging that our model fits the clustering of broad-line quasars well, as long as only luminous AGN above $10^{46}\ \mathrm{erg}\;\mathrm{s}^{-1}$ are included. As the threshold is lowered, the clustering strength decreases, which is a natural prediction from any model in which AGN luminosity is correlated with halo mass. While only heavy seeds are shown, the clustering of light seeds is similar, since these high-luminosity, high-mass AGN are insensitive to the seeding mechanism.

Appendix B Best Fits and Error Estimation

Best fits correspond to the values obtained via $\chi^{2}$ minimization to a joint fit of infrared, X-ray, and cross-power spectra. Since we lack a covariance matrix, all data points are treated independently. This should not affect the best fit values, but may have an effect on our estimated errors. We note that these fits are most strongly affected by the small error bars in the CIB power spectra on small scales. At the same time, this region of the power spectrum is expected to have artifacts due to masking effects. Therefore, we fit only on scales above 30 arcseconds and comment that including data on smaller scales leads to worse values of $\chi^{2}$ and a steeper slope $(\alpha=1.8)$ for the simple power law.

We estimate errors using the emcee MCMC package (Foreman-Mackey et al., 2013). Uninformative priors are used for each of the parameters, ensuring only that the coherence does not exceed unity, for example. Convergence is confirmed by visually inspecting both the evolution of both the median and standard deviation values of each chain’s walkers. Error regions are obtained by encapsulating 68 per cent of the posterior distribution and comparing to the best fit from the $\chi^{2}$ minimization.

Appendix C The Black Hole Mass Density

Since black holes only radiate while they are accreting gas, any proposed flux attributed to AGN implies active black hole growth. A hypothetical population of missing AGN implies that the overall black hole mass density must increase by this new growth channel. Here, we derive the increase in the black hole mass density that is produced by a hypothetical population of AGN at a fixed redshift.

Suppose that the population emits a background in some band with a flux in units of $\mathrm{erg}\;\mathrm{s}^{-1}\;\mathrm{cm}^{-2}\;\mathrm{ster}^{-1}$

[TABLE]

and that at fixed redshift, the comoving (observed) flux per unit volume is

[TABLE]

where $f_{\mathrm{bol}}$ is the fraction of the bolometric luminosity radiated in the observed band. The observed flux can be expressed

[TABLE]

Meanwhile, the Soltan argument (Soltan, 1982) links the population of AGN to the mass density of black holes via

[TABLE]

Combining these results, the mass density can be linked to the observed flux and the population via

[TABLE]

where $z_{\mathrm{eff}}$ is the effective redshift of the population, defined as,

[TABLE]

From this equation, we can see that for any missing AGN population, the amount of extra black hole mass necessary to generate a fixed flux is substantial for all redshifts, scaling only as $(1+z)$. AGN_z2 and AGN_z3 have $z_{\mathrm{eff}}$ of 2.37 and 3.10 respectively. We assign our hypothetical populations of missing AGN $dF/d\Omega$ equal to the best fit in Table 2, $f_{\mathrm{bol}}=0.1$, and $\epsilon=0.1$. For both AGN_z2 and AGN_z3, the amount of extra black hole mass density gained during these AGN episodes is equal to

[TABLE]

This is in excess of current estimates of the black hole mass density (Shankar et al., 2009), and similar to that predicted by (Yue et al., 2013b).