Exploring the High-redshift PBH-ΛCDM Universe: Early Black Hole Seeding, the First Stars and Cosmic Radiation Backgrounds

Dark matter (DM) represents the most abundant form of matter in the universe and dominates the dynamics of collapsed objects. It also offers the scaffolding within which all visible matter is structured into galaxies. Thus far, in the context of the cold dark matter (CDM) paradigm, it has been widely assumed that DM exists in the form of still unknown particles that interact primarily through gravity and perhaps through weak interactions (e.g., Feng 2010). However, despite several decades of targeted experimental searches aimed at uncovering weakly interacting massive particles (WIMPs) as potential DM candidates, these efforts have all come up empty. At present there seem to be no sign of particle DM candidates in the mass interaction cross-section parameter space and energy ranges where they have been predicted (Baudis 2012; Boveia & Doglioni 2018).

Meanwhile, the discovery of gravitational waves from merging black hole (BH) binaries by LIGO and VIRGO reveal surprisingly large masses for the individual merging BHs with, on average, low pre-merger spins. Typical inferred masses of the merging sources are higher than expected from astrophysical formation channels (Abbott et al. 2016). Consequently, the community revived the hypothesis originally proposed by Hawking (1971) that DM could be constituted by primordial black holes (PBHs) that formed in the infant universe (Jedamzik 1997; Carr 2003; Bird et al. 2016; Kashlinsky 2016; Clesse & García-Bellido 2017; Jedamzik 2021).

Early models assumed that all DM is comprised of PBHs that formed with a monochromatic mass function, but observational constraints firmly rule out this hypothesis (e.g., Belotsky et al. 2019). Later work and further refinements, notably by Carr et al. (2019), García-Bellido (2019), Carr et al. (2020), and Carr et al. (2021) showed that DM PBHs can, in principle, have a broad birth mass spectrum ranging from 10⁻¹⁰–10⁷ M_⊙ while accounting for all the DM without violating current observational constraints. In their model, PBHs are created in the early universe during QCD phase transitions (around 100 MeV, which corresponds to ∼10¹⁰ K) involving different particle families freezing out of the primordial quark–gluon plasma within the first 2 s after the inflationary phase. When W^+/−, Z bosons, baryons, and pions are created, and e⁺e⁻ pairs annihilate, they leave an imprint in form of a significant reduction of the sound speed at the corresponding phase transitions, thereby causing regions of high curvature to collapse and form PBHs. The typical mass scale of these PBHs is defined by the size of the horizon at the time of the corresponding phase transition. In this model, four distinct populations of PBHs in a wide mass range are expected to form: planetary mass BHs at the W^+/−–Z transition; PBHs of around the Chandrasekhar mass when the baryons (protons and neutrons) form from three quarks; PBHs with masses of order 30 M_⊙ (these correspond to the suggested LIGO BHs), when pions form from two quarks; and finally PBHs with masses corresponding to those of supermassive black holes (SMBHs) with M ≥ 10⁶ M_⊙ that form at the e⁺e⁻ annihilation. If PBHs form with a broad mass distribution, the DM they constitute is expected to strongly cluster, which would help alleviate some of the more stringent observational constraints on the allowed contribution of PBHs to the DM (Clesse & García-Bellido 2017; Belotsky et al. 2019) budget.

Intriguingly, an excess of small-scale DM substructure compared to CDM predictions has been recently reported from gravitational cluster lensing studies with the deepest Hubble Space Telescope observations (Meneghetti et al. 2020). Clustering of DM in excess of what is predicted by the standard WIMP CDM paradigm as expected with PBH DM, could possibly account for this discrepancy. This excess concentration of mass on small scales is revealed in the discrepancy between the observed and predicted event rates for galaxy–galaxy strong lensing (GGSL) events. The internal structure of sub-halos with masses ∼10¹¹ M_⊙ are implicated for these GGSL events, and ΛCDM simulations simply do not produce enough sub-halos in this range with the requisite central concentrations.

In addition to possibly providing a simple resolution of the persistent DM problem, PBHs it appears could also serve to account for early massive BH seed formation and address the intriguing origin of the SMBHs with mass of the order 10¹⁰ M_⊙ powering detected luminous quasars already in place by z > 7 when the universe was <0.8 Gyr old (see, e.g., Lodato & Natarajan 2006; Li et al. 2007).

Quite a large number of other recent observational results also strengthen the conjecture that PBHs could contribute to the overall DM budget. The latest GWTC-2 catalog of LIGO-Virgo-Kamioka Gravitational Wave Detector (KAGRA) gravitational merger events (Abbott et al. 2021) has widened the observed BH mass distribution considerably. In particular, it includes the most massive merger detected as yet GW190521 (Carr et al. 2019; Clesse & Garcia-Bellido 2021; Abbott et al. 2020a), in which at least one of the two components is more massive than the upper mass gap expected for pair instability supernovae (SNe), and thus could signal a PBH origin (Clesse & Garcia-Bellido 2021; De Luca et al. 2021a). It also includes the event GW190814 (Abbott et al. 2020b), in which one of the components likely falls into the lower SN mass gap between neutron stars and BH, and which has a surprisingly large mass ratio of ∼1:9, not entirely easily compatible with known astrophysical production channels (De Luca et al. 2021a). Wong et al. (2021) analyzed the whole GWTC-2 catalog and concluded that the observed event rate is fully consistent with the assumption that all LIGO-VIRGO detected merging BHs are of primordial origin, contributing a fraction of f_PBH ∼ 0.3% overall to the DM budget in the mass range 1 < M_BH < 100 M_⊙.

Constraints on the fraction of DM contributed by PBHs from X-rays emitted by accretion onto them producing measurable effects on the spectrum and anisotropies of the cosmic microwave background (CMB) were calculated previously by Ricotti et al. (2008). New results from the 5 yr Optical Gravitational Lensing Experiment (OGLE) microlensing campaign (Niikura et al. 2019a) present the discovery of a sizeable population of long-duration microlensing events, which, together with Gaia parallaxes, point to the presence of putative PBHs in the mass range 1–10 M_⊙ (Wyrzykowski & Mandel 2020). Another mass-gap BH candidate was recently discovered in the nearby nearly edge-on ellipsoidal variable binary star V723 Mon (Jayasinghe et al. 2021). OGLE has also detected six ultrashort-timescale microlensing events, which may, in fact, indicate the existence of planetary mass PBHs (Niikura et al. 2019a). The North American Nanohertz Observatory for Gravitational Waves (NANOGrav) pulsar timing observatory has recently reported interesting upper limits on the stochastic gravitational wave background in the nanohertz band in their data, which does not show statistically significant quadrupolar spatial correlations expected for a cosmic gravitational wave background (Arzoumanian et al. 2020) yet. Meanwhile, the Parkes Pulsar Timing Array (PTA) Collaboration has also reported a tentative detection consistent with the NANOGrav finding (Goncharov et al. 2021). However, the derived upper limit they provide is compatible with several PBH formation processes operating over a wide mass range (De Luca et al. 2021b; Vaskonen & Veermäe 2021; Kohri & Terada 2021; Domènech & Pi 2020; Sugiyama et al. 2021). Recent indications from GRB microlensing constraining f_PBH ∼ 3 × 10⁻³ for 10⁶ M_⊙ PBHs (Kalantari et al. 2021), are consistent with the mass spectrum considered here. Finally, microlensing events from multiply lensed quasar images have been interpreted to indicate that the DM in galaxy halos and clusters could be potentially composed of PBH with masses of around 1 M_⊙ (Hawkins 2020a, 2020b). All of these currently reported observational constraints and bounds are plotted in Figure 1. However, we note that these are mostly derived for a monochromatic initial mass spectrum for PBHs and not for the broader mass spectrum we investigate in this work.

Zoom In Zoom Out Reset image size

In addition to this suggestive circumstantial evidence, there are several other open problems in cosmology for which the PBH-DM hypothesis might proffer explanatory power (Carr et al. 2019). The exciting implications for structure formation and evolution at early cosmic epochs in a PBH-DM universe motivates our current detailed exploration of this model. A notable unsolved mystery in observational astrophysics is the origin of the large-scale excess fluctuations in the unresolved cosmic infrared background (CIB) discovered by Kashlinsky et al. (2005), confirmed by Kashlinsky et al. (2007), Cooray et al. (2012b), and Kashlinsky et al. (2012) and their coherence with the cosmic X-ray background (CXB) (Matsumoto et al. 2011; Cappelluti et al. 2013, 2017; Mitchell-Wynne et al. 2016; Li et al.2018). Careful modeling of the shape of these fluctuations, suggests two possible origins: (i) that they are consistent with being produced in the young universe from early stellar emission or (ii) from more local intra-halo light. However, the coherence with the CXB is even more intriguing because it suggests an origin from accretion powered emission from the first BHs, believed to be the progenitors of observed SMBHs. There is no clear consensus on how these SMBH seeds formed and several possible scenarios have been invoked, including the collapse of Population III stars (see, e.g Kashlinsky & Rees 1983) and direct collapse black holes (DCBHs; see reviews by Volonteri 2012; Natarajan 2014). Ricarte et al. (2019) showed that neither DCBHs nor Population III stars could in fact produce the required amount of radiation to explain the observed cross-power spectrum and excess. In addition, the origin of the empirical scaling between the masses of central SMBHs and the stellar velocity dispersion of their host galaxies and DM halos is poorly understood at present, and whether they arise as a result of the seed formation process (nature) or emerge over time as a consequence of growth and assembly in tandem (nurture) is debated (Lodato & Natarajan 2006).

Another potential problem arises for models of emission from Population III stars in the early universe. Extrapolating the star formation density observed at z < 8 to higher redshifts, Cooray et al. (2012a) point out that the expected emissivity of Population III stars falls short in reproducing the required amplitude of the CIB power spectrum. Helgason et al. (2016) studied the physical conditions required of early star-forming galaxies to produce enough flux to explain the observed CIB flux and power spectrum and note that in order to reproduce the observed signal, stars need to form either with an unreasonable baryon conversion efficiency f_⋆ > 0.1 and/or with extremely top-heavy and tilted initial mass functions (IMFs) with all stars being born with masses larger than 500 M_⊙. Even considering other proposals, problems still persist in reconciling the observed source number counts and Planck reionization constraints (Yue et al. 2013).

According to Kashlinsky (2016), if DM is made of PBH then we should expect that at redshifts z > 15 they could accrete a substantial amount of matter to emit enough in the IR (rest-frame UV) and X-ray wavelengths to significantly contribute to the CIB and CXB. Afshordi et al. (2003) pointed out that in a PBH-DM ΛCDM cosmology, the fraction of collapsed halos at high redshifts, specifically at z > 7, is significantly higher than in the standard ΛCDM universe. Due to this excess, even a modest efficiency in converting baryons into radiation could give rise to significantly higher emissivity either through accretion or star formation.

Therefore, computing the detailed history of emission from star formation and BH accretion in a PBH-DM universe over cosmic time is warranted. In a recent paper, Hasinger (2020) predicted the amount of extragalactic background flux produced by accretion onto PBHs by Bondi capture via advection dominated flows (ADAFs). In this scenario, PBHs account for all the DM with a broad mass spectrum peaking near the Chandrasekhar mass of 1.4 M_⊙, following García-Bellido (2019) (see Figure 1). In this event, PBH accretion can produce about 1% of the total [0.5–10] keV X-ray background and about 0.5% of the CIB necessary to explain the observed large-scale, near-infrared (NIR) surface brightness fluctuation power spectrum excess. This emission arises from the recombination process that is spread over redshift from z ∼ 1100, peaking at z ∼ 20 and rapidly fading by z ∼ 7–8. This amount of radiation produced by PBHs and its redshift distribution appears to be consistent with the overall shape of the CIB power spectrum and the intensity required to explain the observed large-scale cross-correlation excess of the CIB and CXB from known populations. These very same accreting PBHs, it turns out, can also produce sufficient radio background emission at high redshift to explain the surprisingly strong sky-averaged redshifted 21 cm line signal observed by the Experiment to Detect the Global EoR Signature (EDGES; Hasinger 2020; Bowman et al. 2018). The signal is reported to be centered at a frequency around 78 MHz and covers a broad range in redshift z = 15–20 (Bowman et al. 2018).

Here, we investigate in detail the contribution of PBHs and consequences thereof at cosmic dawn for the mass assembly history of DM; star formation and reionization history; and BH growth and the production of cosmic radiation backgrounds. While the CXB versus CIB cross-correlation excess might not be produced at high redshift (see, e.g., Cooray et al. 2012b) here we model and investigate in addition to the global properties enumerated above, the observed auto- and cross-power spectra of the CIB and CXB using the predictions of Hasinger (2020). We adopt the standard Press–Schechter formalism to derive the DM halo mass function and combine this with a halo occupation distribution (HOD) approach to compute the clustering strength of PBHs in DM halos. In this paper, we assume a PBH-DM ΛCDM (hereafter referred to as the PBH-ΛCDM model) cosmology with parameters Ω_Λ = 0.7, Ω_m = 0.3, and H₀ = 67 h⁻¹ km s⁻¹ Mpc⁻¹ when needed.

2.1. Mass Power Spectrum and Mass Function of Collapsed Halos

The linear power spectrum of density fluctuations in the PBH-ΛCDM universe has been described by Afshordi et al. (2003) as follows:

$$\begin{eqnarray}&&{P}_{\mathrm{lin}}(k)={P}_{{\rm{\Lambda }}\mathrm{CDM}}(k)+{P}_{\mathrm{Poiss}}(k),\end{eqnarray} \tag{ 1 }$$

where k is the spatial frequency, P_ΛCDM is the power spectrum of the DM, and P_Poiss is an additional component introduced by the discrete nature of PBH DM (Meszaros 1975). This component has little or no effect on large scales but becomes dominant on very small scales in the form of shot noise. Note that with this assumed form, the PBH component is not affected by bias. The Poisson piece can be approximated as

$$\begin{eqnarray}&&{P}_{\mathrm{Poiss}}=\displaystyle \frac{9}{4}{\left(1+{z}_{\mathrm{eq}}\right)}^{2}{n}_{\mathrm{BH}}{\left[g(z)\right]}^{-2},\end{eqnarray} \tag{ 2 }$$

where g(z) is the linear growth factor of fluctuations from z to today, with g(0) = 1; and where n_BH = $\tfrac{{f}_{\mathrm{BH}}}{{M}_{\mathrm{BH}}}$ is the number density of PBHs. As PBHs have a wide range of birth masses, the Poissonian component needs to be evaluated as the density weighted number density of PBHs, which in our assumed case, can be safely represented by the mass scale M_BH=1.4 M_⊙ which corresponds roughly to the peak observed in the mass spectrum plotted in Figure 1. Per standard methodology, with P_M (k, z) the matter density fluctuation power spectrum written out as function of the spatial frequency k and redshift, we can evaluate the rms density contrast over spherical region of co-moving radius r_M and of mass M(r_M ) as

$$\begin{eqnarray*}&&{\sigma }_{M}(z)=\left[\int {P}_{M}(k,z){W}_{{\rm{TH}}}({{kr}}_{M}){k}^{2}{dk}\right],\end{eqnarray*}$$

where W_TH is the standard top-hat window function. Then, using the Press–Schechter formalism (Press & Schechter 1974), we evaluate the fraction of mass in a given volume that is contained in halos of a given mass at redshift z, which is given by

$$\begin{eqnarray}&&{f}_{\mathrm{halo}}(\gt {M}_{h},z)=\displaystyle \frac{1}{2}\,{erfc}\left(\displaystyle \frac{{\delta }_{\mathrm{col}}}{\sqrt{2}{\sigma }_{M}(z)}\right),\end{eqnarray} \tag{ 3 }$$

where δ_col is the over-density required for spherical top-hat collapse.

In Figure 2, we compare the fraction of collapsed halos as function of z in PBH-ΛCDM and in the classic ΛCDMcosmology, using the temperature as the proxy for mass. Conversion from one to the other is straightforward; for example, at z = 10, the halo masses corresponding to the temperature range studied span from M_h = 2.6 × 10⁶ M_⊙–M_h = 2.4 × 10⁸ M_⊙. Notably, we see that in the early universe, at z > 15, the PBH-ΛCDM model predicts a significantly larger number of mini-halos (i.e., M_h ∼ 10⁶ M_⊙–10⁷ M_⊙), than ΛCDM that correspond to T_vir > 2000 K. In the PBH-ΛCDM scenario, mini-halos wherein the first stars form through molecular cooling assemble much earlier than in the ΛCDM universe, shifting the beginning of star formation to a much earlier epoch. Simultaneously, by forming earlier, the first light emitted by these highly biased halos imprints a stronger clustering signal in the power spectrum of the diffuse background fluctuations. This early activity would also inject ionizing photons into the IGM earlier than in the standard star formation models assumed for ΛCDM. We compute the detailed consequences of astrophysical processes that occur in this excess population of early mini-halos in the following sections.

Zoom In Zoom Out Reset image size

3.1. Early Star Formation

First, we evaluate how the significantly higher mini-halo space density at high-z predicted in the PBH-ΛCDM scenario noted above impacts the star formation rate (SFR). To do so, we follow the prescription of Greif & Bromm (2006), as done previously by Helgason et al. (2016). The virial temperature of a halo is related to its mass via

$$\begin{eqnarray}&&M={10}^{8}{M}_{\odot }{\left(\displaystyle \frac{\mu }{0.6}\right)}^{-3/2}{\left(\displaystyle \frac{{T}_{\mathrm{vir}}}{{10}^{4}{\rm{K}}}\right)}^{3/2}{\left(\displaystyle \frac{1+z}{10}\right)}^{-3/2},\end{eqnarray} \tag{ 4 }$$

where μ is the mean molecular weight. Here, we adopt μ = 1.2 consistent with the fact that the universe was mostly neutral at high-z. According to this prescription, the star formation in early galaxies is taken to be proportional to the halo collapse rate at a given redshift as

$$\begin{eqnarray}&&\dot{{\rho }_{\star }}({M}_{h},z)={f}_{\star }\displaystyle \frac{{{\rm{\Omega }}}_{b}}{{{\rm{\Omega }}}_{M}}\displaystyle \frac{d}{{dt}}M\displaystyle \frac{{dn}}{{{dM}}_{h}}(\gt {M}_{\min },t).\end{eqnarray} \tag{ 5 }$$

where f_⋆ is baryon conversion efficiency (i.e., the fraction of baryons that are converted into stars) and M_min is the minimum halo mass that is permitted to shock heat to T_vir and consequently allow gas cooling. Consistent with the treatment of Helgason et al. (2016), we consider two ranges of halo temperatures, each with its own dominant coolant. One population, with 10³ K < T_vir < 4 × 10⁴ K where most of the cooling is due to molecular hydrogen and the other atomic cooling halos with T_vir > 4 × 10⁴ where cooling is mostly due to mono-atomic hydrogen. While the efficiency of star formation, f_⋆ at high-z is unknown (it is scarcely observationally constrained even at more recent epochs), guided by the work of Sun & Furlanetto (2016) and Mirocha et al. (2017), we model the relation between f_⋆ and M_h as

$$\begin{eqnarray}&&{f}_{\star }({M}_{h})=\displaystyle \frac{0.05}{{\left(\tfrac{{M}_{h}}{2.8\times {10}^{11}}\right)}^{0.49}+{\left(\tfrac{{M}_{h}}{2.8\times {10}^{11}}\right)}^{-0.61}}+{f}_{{\star }_{\min }},\end{eqnarray} \tag{ 6 }$$

where ${f}_{{\star }_{\min }}$ is floor for star formation efficiency at low mass. Unlike Helgason et al. (2016), who assumed an universal value for f_⋆, here we adopt a M_h dependent value with ${f}_{{\star }_{\min }}$, that is left as a free parameter.

While in a PBH-ΛCDM universe mini-halos form early, star formation does not necessarily occur in them promptly but is in fact initially suppressed. According to criteria for star formation studied by Stacy et al. (2011), Tseliakhovich & Hirata (2010), and Fialkov (2014) star formation occurs only when baryons can collapse and fragment within the halo. At very high-z, this might not be possible if the circular velocity of the halo v_cool (computed using the analytic prediction of Fialkov 2014) is smaller than the relative streaming velocity v_s of the IGM and dark matter. If v_cool(M_h , z) is the circular velocity (or cooling velocity) of a halo of mass M_h at redshift z, we define a function Θ(M_h , z), which is 0 when v_cool(M_h , z) < v_s and 1 of v_cool(M_h , z) > v_s . Basically, this step function is a simplified representation of the trigger probability of star formation in a halo. Note that v_cool(M_h , z) depends on the IGM temperature and therefore also on the ionization history. As this is also unknown at present, for the purpose of our calculation in this paper, we adopt the fiducial best-fit values proposed by Fialkov (2014). It is worth noting that Kashlinsky (2021) studied the timing of halo collapse in a PBH-ΛCDM universe in the case of a monochromatic PBH mass function. While this proposal would be more pertinent for the physics of PBH DM, with our model we will adapt the simplified approach described in this section because of the large number of concurring variables that might provide a different star formation history (see below).

With this assumption, the star formation density becomes

$$\begin{eqnarray}&&\dot{{\rho }_{\star }}({M}_{h},z)={f}_{\star }({M}_{h})\displaystyle \frac{{{\rm{\Omega }}}_{b}}{{{\rm{\Omega }}}_{M}}\displaystyle \frac{d}{{dt}}M\displaystyle \frac{{dn}}{{{dM}}_{h}}(\gt {M}_{\min },t){\rm{\Theta }}({M}_{h},z).\end{eqnarray} \tag{ 7 }$$

As the very first star-forming halos in the universe are pristine and unpolluted by metals, it is expected that the first stars to form in these Population III mini-halos might be substantially more massive, hotter, and bluer than the subsequent generation of stars. Population III stars are therefore expected to be the main constituents of these lower mass mini-halos until, after a few million years, they explode producing supernovae that pollute the environment with metals thereby modifying the initial mass function for the next generation of star formation. To implement this scenario, we include a function motivated by cosmological simulations of the early star formation reported in Greif & Bromm (2006), that takes into account the metallicity evolution of halos with a parameter, p_pris(z), that represents the fraction of pristine, metal-free halos as function of redshift. We use the model of Greif & Bromm (2006) that demarcates the regimes dominated by Population II and Population III halos, and where the star formation density is taken to be

$$\begin{eqnarray}&&\begin{array}{l}{\dot{{\rho }_{\star }}}_{\mathrm{PopIII}}({M}_{h},z)=\dot{{\rho }_{\star }}{p}_{\mathrm{pris}}(z)\\ {\dot{{\rho }_{\star }}}_{\mathrm{PopII}}({M}_{h},z)=\dot{{\rho }_{\star }}\left[1-{p}_{\mathrm{pris}}(z)\right].\end{array}\end{eqnarray} \tag{ 8 }$$

As the true metallicity evolution of these halos is also unknown, in order to explore the parameter space exhaustively, here we assume three different scenarios by imposing a redshift cut at which a specified fraction of each family of halos is 50% enriched: z_1/2 ∼ = 15, 13, and 11 motivated by the excess number of PBH-DM halos that are available compared to ΛCDM at these epochs, as shown in Figure 2. Going forward, we refer to these three options as early, mid, and late enrichment scenarios, respectively.

3.2. Flux production Rate from Accretion onto PBH DM

After setting up the framework for including star formation, we now proceed to compute the emission from accretion onto PBH DM. In order to evaluate the flux production rate of PBHs, we use the Hasinger (2020) PBH-DM model to evaluate the luminosity of each halo of mass M_h and apply the appropriate bolometric and k-corrections to estimate it in all photometric bands of interest. We note that these factors are typically of the order of ∼10% and are derived, at each redshift from an accretion model—in this case, the ADAF model that aptly describes the expected obscured accretion mode at these early epochs. If f_DM(M_BH) is the fraction of DM in PBHs of a given mass, the occupation fraction of a PBH of mass M_BH in a halo of mass M_h is

$$\begin{eqnarray}&&{n}_{\mathrm{BH}}({M}_{h},{M}_{\mathrm{BH}})={f}_{\mathrm{DM}}({M}_{\mathrm{BH}})\displaystyle \frac{{M}_{h}}{{M}_{\mathrm{BH}}}.\end{eqnarray} \tag{ 9 }$$

The luminosity of each PBH of mass M_BH is then ${L}_{\mathrm{BH}}=\dot{m}\eta {L}_{\mathrm{Edd}},$ where η the radiative efficiency is ∼0.1. Following Hasinger (2020) $\mathrm{Log}(\dot{m})$ is a function of the BH mass and redshift, peaking at around $\mathrm{Log}(\dot{m})$∼−1.6 at z > 20 where the relative streaming velocity between DM and the IGM reaches its minimum value, as demonstrated in high-resolution cosmological simulations by Stacy et al. (2011). We refer the reader to Hasinger (2020) and Stacy et al. (2011) for elaboration on details of the emission properties of PBHs and coupling between baryons and DM at the earliest epochs respectively.

Once the accretion luminosity has been evaluated, the halo luminosity from PBHs can be computed and is given by

$$\begin{eqnarray}&&{L}_{{\nu }_{\mathrm{halo}}}({M}_{h},z)=\int {f}_{\mathrm{DM}}({M}_{\mathrm{BH}}){L}_{{\nu }_{\mathrm{BH}}}({M}_{\mathrm{BH}},z){M}_{h}\,{{dM}}_{\mathrm{BH}},\end{eqnarray} \tag{ 10 }$$

where ${L}_{{\nu }_{\mathrm{BH}}}$ is the normalized bolometric PBH luminosity. The halo, PBH conditional luminosity function (CLF) is then given by

$$\begin{eqnarray}&&{\phi }_{{L}_{{\nu }_{\mathrm{BH}}}}({M}_{h},z)=\displaystyle \frac{{dN}}{{{dM}}_{h}}(z)\end{eqnarray} \tag{ 11 }$$

with $\tfrac{{dN}}{{{dM}}_{h}}(z)$ being the halo mass function (Sheth et al. 2001). The flux production rate is then simply

$$\begin{eqnarray}&&{f}_{{\nu }_{i}}({M}_{h},z)=\displaystyle \frac{{L}_{{\nu }_{\mathrm{BH}}}}{4\pi {D}_{L}{\left(z\right)}^{2}},\end{eqnarray} \tag{ 12 }$$

where D_L (z) is the luminosity distance.

The peak of the expected emission from these mini-halo populations is centered around z ∼ 20, and it extends out to z ∼ 100 in the X-ray. The bulk of this emission arises from low mass (10⁵⁻⁶ M_⊙) halos while no flux is observed in the IR at z > 40 due the Lyman absorption. This early generation of the CXB in the PBH-ΛCDM model has important implications for direct collapse of gas required to form massive BH seeds in conventional ΛCDM.

3.3. The Formation and Growth of Mini-quasars and AGN

With the accretion model adopted above, we note that PBH-ΛCDM provides a recipe for the swift formation of SMBH seeds as they are naturally in place right at the time of halo assembly. In this framework, the most massive BH in the halo sinks to the center and as soon as baryons are available it starts accreting like an AGN (see schematic in Figure 12). As the mass of the largest BH is simply a function of the halo mass by construction in the PBH-ΛCDM model, their masses are inherently correlated (see Equation (9)). This offers a natural qualitative explanation for the origin of the locally observed M_bh–σ relation and accounts for co-evolution between halo assembly, galaxy assembly, and BH growth. This empirical scaling relation informs all current models of the cosmic growth history of BH populations over time (e.g., in Ricarte & Natarajan 2018). The very existence and origin of such a scaling relation, whether it reflects the initial conditions of BH seed formation or is a consequence of co-evolution over time is currently debated. In the PBH-ΛCDM model, the scaling relation naturally arises out of the initial conditions and is sustained by the very nature of DM as PBHs. Since the halo is composed of PBHs, from which the central SMBH originated, and depending on how we model the density profile of the halo, it will scale as σ³⁻⁴. The amplitude of the relation is not predicted explicitly by the model, which is why as noted below, the Bandara et al. (2009) value is adopted.

As the mass of the most massive BH in the halo is proportional to the mass of the halo in PBH-ΛCDM, we calibrate this correlation with the same relation as done by Bandara et al. (2009) using local observational data, where ${\rm{log}}({M}_{\mathrm{BH}})=8.18+1.55* {\rm{log}}({M}_{h})-13$. With the mass of the central BH settled in this manner, we impose growth at the Eddington limit and compute the emission using the template spectral energy distribution (SED) for an obscured SMBH as found in simulations by Pacucci et al. (2015). This choice is driven by the fact that at high-z the fraction of obscured AGN is essentially of order unity (see, e.g., Vito et al. 2014). However, in the case of PBH-ΛCDM, we do not expect all these AGN to accrete simultaneously; therefore, at each redshift, we impose the condition that the fraction of halos hosting an active AGN is given by

$$\begin{eqnarray*}&&{N}_{\mathrm{AGN}}({M}_{h})={\left(\displaystyle \frac{{M}_{h}}{2.19\times {10}^{12}}\right)}^{0.9}+0.023\end{eqnarray*}$$

as modeled by Hassan et al. (2018). This number can be interpreted as either the fraction of SMBHs accreting simultaneously or alternatively as the average Eddington rate if all the SMBHs are accreting, or most likely, a combination of the two. During accretion at these early epochs, as all the AGNs are reasonably assumed to be obscured, we do not consider their contribution to reionization as the bulk of the UV light is absorbed in the torus/envelope. Therefore, stellar sources still contribute the bulk of the photons needed for reionization and even in the PBH-ΛCDM universe replete with copious early accretors, accretion onto BHs does not significantly add to the budget.

With these model assumptions, we find that these early AGN contribute at the level of about 2%–3% of the CXB. In Figure 3 we plot the predicted X-ray luminosity function (XLF) at z = 6 compared with data from the extended ROentgen Survey with an Imaging Telescope Array (eROSITA) (Wolf et al. 2021; Medvedev et al. 2020) and data derived from an extrapolation of the XLF obtained by Vito et al. (2014) at z < 5 using deep Chandra observations. Our model shows remarkable agreement with both these observational determinations, though it predicts a slight excess of sources below L_⋆, where L_⋆ ∼ 10⁴⁵ erg s⁻¹. This likely arises due to lack of data for such faint sources in observational samples that are therefore missing from the extrapolation or it could be due to our general assumptions about the accretion rate and the halo occupation statistics in the model. While further discussion of these faint sources is beyond the scope of this paper, we point out that our model provides a simple solution to the problem of SMBH seed formation without invoking hitherto unobserved physical mechanisms like extreme super-Eddington accretion or the rapid direct collapse of pristine gas.

Zoom In Zoom Out Reset image size

3.4. Stellar and AGN Emission in the PBH-ΛCDM Model at z > 6

3.5. Calculating the Ionization History

In order for the PBH-DM model to be viable, we need to ensure that reionization constraints are strictly met. The enhanced star formation at large redshifts might overproduce the ionizing photon injection rate. Using our SFR model developed in Section 3.1, and combining it with the predictions from the SEDs as outlined before, the injection rate can be written as

$$\begin{eqnarray}&&\dot{n}(z)={\int }_{13.6{h}^{-1}{\rm{eV}}}^{\infty }\int \displaystyle \frac{{J}_{\nu }({M}_{h},z)}{h\nu }\,{dM}\,d\nu ,\end{eqnarray} \tag{ 13 }$$

where J_ν (M_h , z) is the emissivity per unit halo mass and redshift. We estimate the IGM ionization fraction using

$$\begin{eqnarray}&&{\dot{x}}_{\mathrm{ion}}=\displaystyle \frac{{f}_{\mathrm{esc}}\dot{}n(z)}{\langle {N}_{{\rm{H}}}\rangle }-\displaystyle \frac{{x}_{\mathrm{ion}}}{{t}_{{\rm{rec}}}},\end{eqnarray} \tag{ 14 }$$

where 〈N_H〉 is the mean hydrogen column density and t_rec is the recombination timescale. From these, we can derive the Thompson optical depth via

$$\begin{eqnarray}&&{\tau }_{e}=c{\sigma }_{T}\langle {N}_{{\rm{H}}}\rangle =\int {x}_{\mathrm{ion}}(1+\displaystyle \frac{\eta Y}{4X}{\left(1+z\right)}^{2})\displaystyle \frac{{dt}}{{dz}}\,{dz},\end{eqnarray} \tag{ 15 }$$

where X and Y are the hydrogen and helium abundance (X = 0.76, Y = 0.24), respectively, and η = 1 at z > 3 (helium ionized once) and η = 2 at z < 3 (helium fully ionized). In this work, we focus exclusively on sources at z > 6; therefore, in order to account for later ionization history, we derive ${\dot{n}}_{\mathrm{ion}}$ for z < 6 from the measured star formation histories collated in Robertson et al. (2015).

3.6. Computing IR Source Counts

Having developed a model for growth of structure and emission, we can ask ourselves how many of these faint early sources in PBH ΛCDM we can expect to detect in IR surveys. Accretion onto PBHs produces very faint X-ray emission, so we estimate the number counts of proto-halos/galaxies by integrating the CLF in the IR. First we derive the luminosity function

$$\begin{eqnarray}&&\phi ({L}_{\nu },z)={\int }_{{M}_{\mathrm{Min}}}^{{M}_{\mathrm{Max}}}\phi ({L}_{\nu }({M}_{h},z)){dM},\end{eqnarray} \tag{ 16 }$$

from which we then obtain the number counts

$$\begin{eqnarray}&&N(m)={\int }_{0}^{\infty }\phi ({m}_{\nu },z){dV}(z){dz},\end{eqnarray} \tag{ 17 }$$

where dV(z) is the co-moving volume element at redshift z. Observed source counts from the soon to be launched James Webb Space Telescope (JWST), in the the [2–5] μm band with Near Infrared Camera (NIRCAM) offer the prospect for placing tight constraints on our model as reported in the Results section of this paper.

3.7. The Production Rate of Cosmic Backgrounds

In order to evaluate the integrated cosmic backgrounds produced by these additional high-z mini-halo populations in the PBH-ΛCDM universe, we first need to compute the emissivity of individual star-forming halos for which we deploy spectral templates from the Yggdrasil model (Zackrisson et al. 2011). Yggdrasil creates synthetic spectra that include a single age stellar population and a mix of Population II and Population III stars. The parameters of the model are simply the initial mass function (IMF) and metallicity. With these parameters the spectra are modeled to include nebular emission and extinction/reprocessing from dust. As done in Helgason et al. (2016), we assign different IMFs to the halos according to their masses: for halos with T_vir > 4 × 10⁴ K, and if they are already enriched, we adopt a Population II Kroupa IMF in the mass range 0.1–10 M_⊙ with a metallicity Z = 0.0004 (hereafter, our Model IMF-2). The other free parameters of the model are the escape fraction f_esc and the star-forming efficiency floor ${f}_{{\star }_{\min }}$. We refer to Helgason et al. (2016) and Zackrisson et al. (2011) for a detailed description of the spectral templates. Cooling in the lower mass halo population is driven by H₂ cooling, so for these halos, we assume that the IMF is the same as assumed for the formation of Population III stars. To model the IMFs in these lower mass mini-halos, we explore three zero metallicity IMFs, a top-heavy IMF with a cutoff at 500 M_⊙ (IMF-3A), a lognormal with a characteristic mass of 10 M_⊙ and σ_M = 1 (IMF-3B), and a Kroupa IMF with a cutoff at 100 M_⊙ (IMF-3C), with these models referred to as IMF-3A, IMF-3B, and IMF-3C respectively. For each SED, the emissivity per unit DM halo mass M_h per unit redshift is given by

$$\begin{eqnarray}&&{j}_{\nu }({M}_{h},z)={\int }_{z}^{\infty }\dot{{\rho }_{\star }}({M}_{h},{z}^{{\prime} }){L}_{\nu }(t-{t}_{z}^{{\prime} })\displaystyle \frac{{dt}}{{{dz}}^{{\prime} }}\,{{dz}}^{{\prime} },\end{eqnarray} \tag{ 18 }$$

where ${L}_{\nu }(t-{t}_{z}^{{\prime} })$ is derived from the evolving spectral templates from Yggdrasil. For a given wavelength, we can now derive the flux production rate,

$$\begin{eqnarray}&&f({M}_{h},z)=\displaystyle \frac{c}{4\pi }{j}_{\nu }({M}_{h},z)\displaystyle \frac{\nu }{1+z}\displaystyle \frac{{dt}}{{dz}}.\end{eqnarray} \tag{ 19 }$$

We now have the f(M_h , z) necessary for Equation (31) to derive P₂(q). Similarly, we can derive the CLF by deriving the halo luminosity

$$\begin{eqnarray}&&L({M}_{h},z)=\displaystyle \frac{\nu {J}_{\nu }({M}_{h},z)}{({M}_{h}{dn}/{{dM}}_{h}(z))},\end{eqnarray} \tag{ 20 }$$

from which we obtain that $\phi (L({M}_{h},z))=\tfrac{{dN}}{{{dM}}_{h}(z)}$. Since Yggdrasil produces only results for rest-frame wavelengths longer than the UV, we derive the X-ray emissivity and CXB production rate using an extrapolation of the L_X–SFR relation of Aird et al. (2017):

$$\begin{eqnarray}&&{\rm{log}}({L}_{{\rm{X}}})=39.48+0.83\,{\rm{log}}(\mathrm{SFR})+1.34\,{\rm{log}}(1+z)\end{eqnarray} \tag{ 21 }$$

and convert it into flux. This allows us to compute the CXB and CIB produced by this source population in the PBH-ΛCDM universe.

3.8. Calculation of Fluctuations in the Cosmic Backgrounds

With the PBH-DM halo populations and their emissivity from star formation and AGN activity computed, we then use an appropriately modified version of the HOD formalism to compute the integrated surface brightness and fluctuations therein. At a given frequency ν, we define the surface brightness fluctuations of the diffuse component of background radiation as

$$\begin{eqnarray}&&\delta {F}_{\nu }(x)={F}_{\nu }(x)-\langle {F}_{\nu }\rangle ,\end{eqnarray} \tag{ 22 }$$

where F_ν is the surface brightness at the position x in the sky and 〈F_ν 〉 is the mean integrated flux of the background at the frequency ν. We evaluate the power spectrum of fluctuations defined as P(q) = 〈∣δ F _ν (q)∣〉 with q as the angular wavenumber and where δ F _ν (q) is the Fourier transform of the 2D fluctuation field. Here, we adopt the formalism derived by Ricarte et al. (2019). The power spectrum can be written as a sum of the one- and two-halo terms:

$$\begin{eqnarray}&&P(k)={P}^{1h}(k)+{P}^{2h}(k).\end{eqnarray} \tag{ 23 }$$

These two components represent contributions to flux from within and from clustering outside the halo, respectively. For each source population (see below) we define

$$\begin{eqnarray}&&{F}_{\nu }^{c}=\langle {f}_{\nu }^{c}({M}_{h},z)\rangle ,\end{eqnarray} \tag{ 24 }$$

and

$$\begin{eqnarray}&&{F}_{\nu }^{s}=\langle {f}_{\nu }^{s}({M}_{h},z)\rangle .\end{eqnarray} \tag{ 25 }$$

where ${f}_{\nu }^{c}$ and ${f}_{\nu }^{s}$ are background production rates for each source population of PBHs within halos of given mass in central (denoted by the superscript c) and satellite sub-halos (denoted by the superscript s) at a frequency ν, respectively. For the populations described here, we assume that all PBH sources, even the ones that constitute the parent halo are essentially satellites except for the single most massive central BH that accretes like an AGN. We assume that PBHs are distributed within individual halos in accordance with the expected Navarro–Frenk–White profile (NFW; Navarro et al. 1997) predicted by CDM. The background production rate in a given M_h bin and redshift is described in Section 3.2. As each population traces dark matter differently, their two-halo terms are then given by

$$\begin{eqnarray}&&{P}_{i,j}^{2h}{(k,z)=b{({M}_{h},z)}_{i}\,\ast \,b({M}_{h},z)}_{j}\,\ast \,{P}^{\mathrm{lin}}(k,z),\end{eqnarray} \tag{ 26 }$$

where P^lin(k, z) is the PBH-ΛCDM power spectrum from Equation (1) and $b{({M}_{h},z)}_{i}$ is the linear halo bias (Sheth et al. 2001), for the ith population. As we are dealing with diffuse emission rather than individual point sources, bias is now computed by flux weighting central and satellite emissivities such that

$$\begin{eqnarray}&&{B}^{c}=\int \displaystyle \frac{{dn}}{{{dM}}_{h}}\displaystyle \frac{{F}_{\nu }^{c}}{\bar{{F}_{\nu }}}b({M}_{h},z){{dM}}_{h}\end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray}&&{B}^{s}=\int \displaystyle \frac{{dn}}{{{dM}}_{h}}\displaystyle \frac{{F}_{\nu }^{s}}{\bar{{F}_{\nu }}}b({M}_{h},z)u(k| {M}_{h},z){{dM}}_{h},\end{eqnarray} \tag{ 28 }$$

where u(k∣M_h , z) is the Fourier transform of the NFW profile. We can then write out the one-halo term as

$$\begin{eqnarray}&&\begin{array}{l}{P}_{i,j}^{1h}(k,z)=\displaystyle \int \displaystyle \frac{{dn}}{{{dM}}_{h}}\\ \quad \times \displaystyle \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}}{{dM}}_{h},\end{array}\end{eqnarray} \tag{ 29 }$$

and the two-halo term as

$$\begin{eqnarray}&&{P}_{i,j}^{2h}(k,z)={P}^{\mathrm{lin}}(k,z)({B}_{i}^{c}+{B}_{i}^{s})({B}_{j}^{c}+{B}_{j}^{s}).\end{eqnarray} \tag{ 30 }$$

The indices i, j label the photometric bands so that if i = j we obtain the auto-power spectrum in a single band while if i ≠ j, one obtains the cross-power spectrum. Next, we reconstruct the power spectrum using Limber's equation, which projects the 3D power spectrum folded in with the emissivity of each halo

$$\begin{eqnarray}&&\begin{array}{l}{P}_{2}(q)=\displaystyle \int \displaystyle \frac{H(z)}{{{cd}}^{2}(z)}\displaystyle \int {f}_{{\nu }_{i}}({M}_{h},z){f}_{{\nu }_{j}}({M}_{h},z)\\ \quad \times \,{P}_{{3}_{i,j}}({{qd}}_{c}{\left(z\right)}^{-1}){dM}\,{dz},\end{array}\end{eqnarray} \tag{ 31 }$$

where, q is the angular wavenumber in rad⁻¹; d_c (z) is the co-moving distance; H(z) is the expansion rate of the universe at redshift z; ${f}_{{\nu }_{i,j}}({M}_{h},z)$ is the flux produced per unit halo mass and redshift at a frequency ν_i,j ; and ${P}_{{3}_{i,j}}({{qd}}_{c}{(z)}^{-1})$ is the DM power spectrum described in Equation (1) expressed a function of q. We now have all the requisite machinery in place to estimate the unresolved background fluctuations in the X-ray [0.5–2] keV and [2–5] μm bands, their auto- and cross-power spectra.

3.9. Auto and Cross-correlations of Cosmic Backgrounds

A powerful set of observables predicted by our model are the CIB and CXB fluctuation power spectra and their cross power. We have computed these by including the flux production rate in Equations (5) and (6) by assuming that stars and the general population of growing PBHs are distributed as satellites and that the central sources—the most massive PBHs in their parent halos—grow like an AGN. Our final angular power spectra are then obtained as follows:

$$\begin{eqnarray}&&\begin{array}{l}{P}_{i,j}(q)={P}_{i,j}^{s\mathrm{tars}}(q)+{P}_{i,j}^{\mathrm{QSO}}(q)\\ \quad +{P}_{i,j}^{\mathrm{PBH}}(q)+{{XP}}_{i,j}(q)+{P}^{\mathrm{SN}}.\end{array}\end{eqnarray} \tag{ 32 }$$

We enclose in XP(q) the cross-power terms for each pair of populations and bands.

We have added as foregrounds shot noise and clustering from unresolved local galaxies, AGN, and clustering both in the IR and X-ray bands using the predictions of Helgason et al. (2012, 2014) who extrapolated known luminosity functions of galaxies, AGN, and clusters to fluxes below current survey flux limits. Here, we have assumed a fiducial flux limit in the X-ray of ${S}_{\mathrm{lim},0.5-2}={10}^{-17}$ erg cm⁻² s⁻¹ in the [0.5–2] keV band and m_AB = 24.5 in the [2–5] μm band.

As with all treatments that involve extrapolating current observations to higher redshifts, to z > 7, be it in the context of the standard ΛCDM model or variants like the PBH-ΛCDM that we consider here, assumptions need to be made regarding astrophysical processes relevant to star formation, BH accretion, and the corresponding emitted flux at these extremely early epochs. As a result, unsurprisingly, we find that there exist classes of PBH-ΛCDM models that are viable and consistent with current observational constraints. Our first key finding, modulo our assumptions detailed in Section 3, is that there are two parameters that drive the best fits for the class of feasible models and they are f_⋆, the efficiency floor for star formation, and f_esc the escape fraction of ionizing photons, i.e, quantities that are empirically ill-constrained even at much later cosmic epochs.

In order to determine the optimal combination of parameters for our models that are in agreement with observed constraints, we sampled the parameter space using the Python package emcee (Foreman-Mackey et al. 2013) that uses the affine-invariant ensemble sampler for the Markov chain Monte Carlo (MCMC). We run these chains assuming flat prior conditions on all our parameters (i.e., all the values sampled have the same prior probability). The priors for our models are listed in Table 1 in the form of intervals outside which the likelihood becomes−∞.

Table 1. Best-fit Results for Our MCMC Realizations for the Three Assumed Enrichment Histories and the Three Assumed Population III IMFs

IMF-3A	Early	Mid	Late	Priors
Log(ηIII)	−${3.41}_{-1.07}^{+1.44}$	−${3.63}_{-0.91}^{+1.40}$	−${3.34}_{-1.19}^{+1.56}$	(−5, −1)
Log(fescIII)	−${2.61}_{-0.90}^{+1.62}$	−${2.58}_{-0.98}^{+1.49}$	−${2.41}_{-1.09}^{+0.97}$	(−4, 0)
Log(fescII)	−${2.23}_{-0.99}^{+0.34}$	−${2.27}_{-0.64}^{+0.33}$	−${2.79}_{-0.83}^{+0.75}$	(−3.5, −1)
Log(ηII)	−${2.22}_{-0.28}^{+0.22}$	−${2.21}_{-0.23}^{+0.17}$	−${2.17}_{-0.47}^{+0.20}$	(−4, 0)
IMF-3B
Log(ηIII)	−${3.46}_{-1.09}^{+1.44}$	−${2.94}_{-1.33}^{+1.22}$	−${3.21}_{-1.22}^{+1.44}$	(−5, −1)
Log(fescIII)	−${2.43}_{-0.77}^{+1.04}$	−${2.43}_{-0.77}^{+1.04}$	−${2.55}_{-0.88}^{+1.14}$	(−4, 0)
Log(fescII)	−${2.11}_{-0.37}^{+0.32}$	−${2.22}_{-0.97}^{+0.37}$	−${2.26}_{-0.64}^{+0.45}$	(−3.5, −1)
Log(ηII)	−${2.24}_{-0.38}^{+0.21}$	−${2.24}_{-0.30}^{+0.19}$	−${2.25}_{-0.65}^{+0.27}$	(−4, 0)
IMF-3C
Log(ηIII)	−${3.14}_{-1.19}^{+1.43}$	−${3.01}_{-1.25}^{+1.11}$	−${3.17}_{-1.24}^{+1.50}$	(−5, −1)
Log(fescIII)	−${2.40}_{-1.11}^{+1.60}$	−${2.43}_{-1.03}^{+1.3}$	−${2.12}_{-1.19}^{+1.44}$	(−4, 0)
Log(fescII)	−${2.10}_{-0.34}^{+0.26}$	−${2.06}_{-0.33}^{+0.26}$	−${2.19}_{-0.73}^{+0.34}$	(−3.5, −1)
Log(ηII)	−${2.24}_{-0.33}^{+0.18}$	−${2.28}_{-0.24}^{+0.21}$	−${2.20}_{-0.44}^{+0.25}$	(−4, 0)

Download table as:  ASCIITypeset image

To create our MCMC chains we employed 30 walkers and 1000 steps, from which we then define a likelihood function to determine the posterior probability in the following usual form:

$$\begin{eqnarray}&&\mathrm{ln}P=-\displaystyle \frac{1}{2}+\displaystyle \frac{{\left(\tau -{\tau }_{m}\right)}^{2}}{{\sigma }_{\tau }}+\displaystyle \frac{{\left({\rm{SFR}}-{{\rm{SFR}}}_{m}\right)}^{2}}{{\sigma }_{{\rm{SFR}}}},\end{eqnarray} \tag{ 33 }$$

where the index m represents the mth realization of the model; τ is the Thompson optical depth measured by the Planck Collaboration et al. (2020); and SFR is the star formation density at z = 7 estimated by Bouwens et al. (2020). Additionally, we impose the condition that the total X-ray background CXB produced by our sources in the PBH ΛCDM does not exceed the unresolved component in the [0.5–2] keV band constrained by Cappelluti et al. (2017). In the event that a model exceeds the observationally constrained CXB level, we set $\mathrm{ln}P=-\infty$. We then repeat the sampling for each combination of choice of Population III IMFs and enrichment history with the following models: IMF-3A (top-heavy Population III IMF with a cutoff at 500 M_⊙), IMF-3B (lognormal Population III IMF with a characteristic mass of M ∼ 10 M_⊙), IMF-2 (Population II Kroupa IMF in the mass range 0.1–10 M_⊙ with a metallicity Z = 0.0004), and IMF-3C (Population III Kroupa IMF with a cutoff of 100 M_⊙). The free parameters of the overall model are f_⋆ and f_esc for each population: η_III, f_esc,III, η_II and f_esc,II for the Population III and and Population II halos, respectively. We ran independent fits for each combination of choice of Population III IMF template and enrichment history, for a total of nine explorations. The best-fit results for each combination are summarized in Table 1. The posterior distributions for our MCMC fits are shown in Figure A1 in the Appendix.

Overall, since the fit relies on and is driven only by two parameters, several scenarios are feasible and can reproduce the observed universe and are comparable with current observational constraints. Our results provide clarity on the possible permitted scenarios. Now, we focus on key predictions of the model scenarios that are viable and consistent with the current scant, high-redshift observational constraints. From Figures 4–11, it is clear that the fits provide significant constraints for each parameter with the notable exception of the late enrichment scenario, in which η_III and f_esc,III are degenerate. The region with the parameter space at high η_III and high f_esc,III can be excluded conclusively by our model, likely due to the fact the such a combination would produce too many ionizing photons. For η_II and f_esc,II our fits point to a feasible scenario where the star formation efficiency is of the order of a fraction of a percent and the escape fraction is significantly lower than 10%. By construction the PBH-ΛCDM universe explored here patches onto the ΛCDM universe at z ≤ 6. In order to keep the model economical we do not include a halo mass dependent escape fraction, and therefore these reported values must be considered as an average for the population. While constraining the properties of early star formation is an open question in astrophysics, addressing it any further other than with these simple parameterizations is beyond the scope of this work. Our goal has been to to explore the possible permitted realizations of an early universe within a simple PBH-ΛCDM cosmology that has explanatory power.

4.1. Detailed Model Predictions

4.1.1. SFR Density and Accreted Mass Density

In Figure 4 we show the predicted SFR density as function of redshift for all realizations with a 1σ dispersion as drawn from our fits. In the main text we show the case for IMF-3A while the other three cases IMF-2, IMF-3B, and IMF-3C are reported in Appendix B. We also show the three cases with different enrichment histories—the early, mid, and late considered here that all smoothly connect with the current SFRD plotted in green, determined from observations of extragalactic background light (EBL) (Fermi-LAT Collaboration et al. 2018) and high-z galaxy surveys (best fit and data from Madau & Dickinson 2014). We note that by construction, all model variants are tuned to reproduce the z ∼ 6–7 star formation rate density (SFRD); therefore, our predictions are primarily for z > 7. The density of model lines plotted represents the model sampling density. As can be seen, while all the three IMFs produce a similar dispersion overall, the specific track families predicted for each enrichment model differ. The generic prediction is the existence of two clear peaks in the SFRD between 10 < z < 30. Each enrichment model has a distinct peak. The realizations suggest that in the case of early enrichment the bulk of the star formation happens in mini-halos at very high-z and Population III halos become dominant at z < 10. In the late enrichment case, on the other hand, most of the Population III activity is delayed and therefore it dominates the reionization epoch.

Zoom In Zoom Out Reset image size

4.1.2. Thompson Optical Depth

For the assumed star formation history in the PBH-ΛCDM model, we can compute the predicted redshift evolution of the Thomson optical depth for the three enrichment models and for the class of explored IMFs. For IMF-3A, we show these in Figure 5, where current constraints from Planck (Planck Collaboration et al. 2020) are shown as the green band. As noted above, many classes of models with various combinations of values of f_esc and f_* considered here are, in principle, consistent with current observations.

Zoom In Zoom Out Reset image size

4.1.3. CIB and CXB Flux Density

With the emission from star formation and BH accretion modeled as described in Sections 3.3 and 3.4, we calculate the cumulative flux densities for the resultant CIB and CXB from the entire population at z > 7, as outlined in Section 3.7. The contribution from accreting AGN is shown with the black dashed line and that from Bondi accretion onto PBHs is shown with the black dotted line in Figure 6 for the model IMF-3A. The other two Population III IMF cases, IMF-3B, and IMF-3C are reported in the Appendix. The clear picture that emerges from these PBH-DM models is that the bulk of the flux contribution to the CIB is provided by star formation activity and accretion processes are a subdominant component. All the model predictions produce <1 nW m⁻² sr⁻¹ of the CIB. This value is what is generally assumed to be necessary to safely account for the observed fluctuations. However, our model does not rule out a scenario in which first light commences with an early highly efficient star formation episode, followed by a less efficient Population II episode. In this scenario, the escape fraction tends to be extremely low in Population III stars, of the order 10⁻³, suggesting that such an episode must occur in an optically thick environment from the point of view of UV light production.

Zoom In Zoom Out Reset image size

Similarly, we compute the contribution to the CXB from these z > 7 sources, and these are shown in Figure 7, where the black dashed line again shows the contribution from AGN and the black dotted line the contribution from streaming baryons. The current limits on the unresolved component of the CXB is shown as the cyan band. At z > 15 the bulk of the CXB is produced by PBH accretion, while at z < 15 AGN hosted at the centers of PBH-DM halos become the dominant contributor.

Zoom In Zoom Out Reset image size

4.1.4. Source Counts

Applying appropriate bolometric corrections, we compute the predicted NIR source counts. At faint magnitudes m_AB ∼ 28–30 the contributions from PBH-ΛCDM sources completely dominates. In Figure 8, we show our model predictions along with high-redshift extrapolations of the faint end of current observed IR luminosity functions derived using population synthesis models produced by Helgason et al. (2012). It turns out that these computed faint IR source counts offer the most stringent test of the PBH-ΛCDM model. Observations from NIRCAM on board JWST are expected to detect sources in precisely this magnitude range m_AB ∼ 28–30 where our models clearly predict a significant excess in counts as well as an unambiguous, strong steepening of the number count slope.

Zoom In Zoom Out Reset image size

In Figure 8, we compare the deep Spitzer-Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (S-CANDELS) results (Ashby et al. 2015) with a prediction of the 4.5 μm number counts from our model for the Population III IMFs considered here and the three enrichment histories and assuming f_⋆ = 0.005 and an escape fraction f_esc = 0.1. The predictions are compared with the extrapolations of the population synthesis model of Helgason et al. (2012) based on the faint end of the luminosity function. Star formation in the PBH-ΛCDM model does not produce enough flux to significantly change the shape of the source counts at m_IRAC < 28–30, but at the very faint end (i.e., m ≳ 28–30) we predict, that due to the rapidly steepening $\mathrm{log}N-m$ relation, early star-forming objects become the dominant population. Interestingly, this steepening lies right at the expected depth of the JWST 10 ks survey like the JWST-NEP (Jansen & Windhorst 2018). Therefore, this important model prediction stands to be tested very soon.

4.1.5. CIB and CXB Angular Auto- and Cross-power Spectra

The predicted power spectrum of the unresolved [2–5] μm CIB fluctuations from our model are shown in Figure 9 and that of X-ray sources in the [0.5–2] keV range modeled here as accreting PBHs + X-ray binaries are shown in Figure 10. These are compared with the expected power from the hitherto unknown population producing the CIB and CXB joint fluctuations as derived from their coherence by Kashlinsky et al. (2018). Once again, we show the realizations of the early, mid, and late enrichment models for IMF-3A, while the case for models IMF-3B and IMF-3C are reported in Appendix B. We have focused on the IMF-3A model and chosen to elaborate on it as the observed excess can be accounted for only in realizations of IMF-3A, regardless of the enrichment history. We point out that for IMF-3C, the power spectrum predictions in the late enrichment case produce a sizeable fraction of realizations where the shot noise component exceeds the measure so we can safely conclude that IMF-3C provides the least likely scenario.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The cross-power spectrum of the unresolved [2–5] μm CIB versus [0.5–2] keV CXB fluctuations from our model added to the foreground galaxy estimate from Helgason et al. (2012) is shown in Figure 11. The data points overplotted are combined results from Cappelluti et al. (2013, 2017). We plot the three enrichment models and the three IMFs as before. Only a few realizations from IMF-3A and a few from IMF-3B with the late/mid enrichment model can explain the detected excess power on large angular scales in CIB versus CXB cross power spectrum. As a side note, our modeling shows that without the AGN component introduced into our model, the X-ray auto power and the cross power would have been significantly weaker than what is shown in Figures 10 and 11. Therefore, our model strongly suggests the existence of such a population of sources rising at z < 15.

Zoom In Zoom Out Reset image size

In this paper, we have explored the high-redshift properties of a universe in which DM matter is composed of a broad mass spectrum of PBHs. The main picture that we postulate in this paper is summarized in Figure 12.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The key results of our investigation of the PBH-ΛCDM model can be summarized as follows:

• 1.  
  At z > 6 low mass halos are more abundant than in the classical ΛCDM model even though a large fraction of these halos do not satisfy the physical conditions required for commencing star formation. The relative velocity of the IGM and DM suppresses star formation at M_h < 10^5–6 M_⊙. However, depending on the enrichment history, our model predicts a secondary peak of star formation at z ∼ 15–20 beyond the well-established observed peak at low redshifts, at z ∼ 3. This second peak at very early epochs is found to be driven by mini-halos that most likely host the first episode of Population III star formation. In terms of the overall star formation history and its evolution during the period that we have observational constraints, our model successfully reproduces the SFRD observed in deep surveys by, e.g., Bouwens et al. (2020) and modeling by Madau & Dickinson (2014). This means that in terms of SFRD at z < 6–7 our model predicts a scenario that is basically indistinguishable from predictions and measurement of ΛCDM.
• 2.  
  In order to explain the observed fluctuations in the CIB and CXB as a corollary of the PBH-ΛCDM scenario, we assumed that SMBHs are simply the most massive objects in each halo and by construction their mass scales as the mass of the host halo across cosmic time. These SMBHs then grow as obscured AGN with a product of their Eddington rate times the halo occupation fraction of a few percent. These assumptions lead to significant boosting of the X-ray fluctuations and hence of the unresolved CIB versus CXB cross-power spectrum with a modest effect on the CIB fluctuations alone. In terms of accounting for the power in fluctuations, accreting AGN turn out to be lead players, as they inhabit higher mass halos, which are extremely rare and biased, and are at the same time, extremely luminous. The combination of these two properties makes AGN the objects of interest in our quest for the missing CXB power.
• 3.  
  In the PBH-ΛCDM cosmology we observe early production of the X-ray background due to the early commencement of Bondi accretion from the streaming IGM onto PBHs. This emission accounts for about 1% of the unresolved CXB. A negligible fraction of the CXB is produced by XRBs during the star formation bursts associated with early star formation peaks. At z < 15 the unresolved CXB becomes dominated by the onset of AGN activity that is a natural consequence of the PBH nature of DM.
• 4.  
  A hard constraint that any model of the early universe has to comply with is the Thompson optical depth that must not exceed measurements from current cosmic microwave background (CMB) experiments. Our model does not follow any detailed prescription for the reionization history, and we simply use Planck Collaboration et al. (2020) observations directly to drive the fit. While our model in general satisfies the observational limits, it produces several possible reionization histories depending on the combination of SFRD and f_esc with high early star formation efficiency and low f_esc less favored. In fact, a high star formation efficiency paired with extremely low f_esc in Population III halos would trigger extremely early reionization. Further observations, in particular at 21 cm are a required addition to the model to predict detailed ionization histories to constrain the model.
• 5.  
  A polarization signal is predicted to be produced in the CMB during reionization. The large-scale E-mode CMB polarization can therefore serve as a very sensitive probe of the reionization history (Reichardt 2016). If reionization starts significantly earlier than is classically assumed via instantaneous recombination around z = 6–7, there will be a detectable peak in the EE component of the polarization at higher multipoles than predicted in the case of ΛCDM (see, e.g., Figure 8 of Kashlinsky et al. 2018). Similarly, the additional radio background produced by accreting PBHs at higher-z can imprint a signal in the spin temperature of neutral hydrogen that could be discerned in the integrated, redshifted 21 cm signal from high-precision low frequency radio wave measurements extending the EDGES experiment (Bowman et al. 2018).
• 6.  
  A common feature of our predictions across all models is a steepening of the NIR source counts at m_AB > 28−30. Our model clearly shows that if the IMF of Population III stars is tilted to be top-heavy, the steepening occurs at brighter magnitudes than in the case of classic Population III Kroupa IMF or the lognormal Population III IMF. This steepening appears at the 10 ks exposure limit for JWST. However, deeper observations and exploiting the lensing magnification of sources behind clusters will be extremely powerful in detecting these fainter counts.
• 7.  
  The PBH DM with the wide mass distribution assumed here will cluster around supermassive PBHs in the centers of DM halos. Similar to the dynamical hardening processes in stellar globular clusters, many of these PBHs will merge through triple interactions while others will be ejected. In particular, mergers between intermediate-mass BHs and SMBHs, as well as extreme mass ratio inspirals will be more frequent and will occur at earlier redshifts than expected for classical ΛCDM. The gravitational wave observatories LISA and PTA will be sensitive to the signal from these expected high redshifts and can also provide strong constraints on the PBH-DM scenario presented here.
• 8.  
  As demonstrated, the excess CIB fluctuations can be explained by a stellar origin in some realizations of IMF-3A, and some additional amplitude can be recovered by adding an accreting AGN component. However, while we showed that by adding an AGN component we satisfy the requirements posed by the CIB versus CXB cross power, full accounting for the required integrated 1 nW m⁻² sr⁻¹ measured signal (see, e.g., Helgason et al. 2016; Kashlinsky et al. 2018) is not achievable with a stellar component alone. Some of this can be mitigated by including modeling potential additional unaccounted components like a more detailed measurement of the cirrus and/or an IHL component. In the case of the PBH-ΛCDM scenario explored here, where the typical mass of the PBHs is assumed to be of the order 1 M_⊙, the excess number of mini-halos is insufficient to produce enough stars to boost the CIB fluctuations on large angular scales. However, other scenarios with larger average PBH masses that have not been explored here, might produce more star-forming halos at even higher masses at lower redshift (and hence with larger bias) and therefore impact the fluctuations significantly. Such high masses for PBHs, for instance, can be obtained by assuming that PBH-PBH mergers shift the mass distribution shown in Figure 1 to higher masses over time and therefore lead to the collapse of more massive halos—a modifying effect that we plan to pursue in future work. Finally, it bears noting that because of the high shot noise predictions, the majority of IMF-3C realizations are ruled out.
• 9.  
  It turns out that the X-ray spectral energy distribution of the CXB/CIB cross-correlation signal also contains information about their production mechanism, for example, on the amount of absorbing material around the accreting black holes. Li et al. (2018) show a model comparison for the averaged signal from all Chandra deep fields, emphasizing that in particular the signal in the softest X-ray band (0.5–1 keV) can provide useful discrimination of the sources of production. However, the statistical quality of these data are insufficient to do so at the present time. Future observations, in particular the cross-correlation of the data from eROSITA deep fields—and those from the planned Advanced Telescope for High ENergy Astrophysics (ATHENA) satellite, along with Euclid NIR observations, will in combination provide powerful constraints and tests of our model (Kashlinsky et al. 2019).
• 10.  
  We note that the models explored here look like DCBH models with early BH formation as AGN turn on very early in this scenario too. The PBH-DM models predict an initial occupation fraction of unity—one BH per galaxy naturally. This of course, stands to be modified at later times due to dynamical interactions arising from binary black hole mergers. The DCBH formation channel stands to be tested with data from the NIRCAM and Mid-Infrared Instrument (MIRI) detectors on board JWST, that will soon be available Natarajan et al. (2017).
• 11.  
  In terms of further AGN predictions of the PBH-ΛCDM model, we note that extrinsic AGN variability on multiple timescales due to ubiquitous microlensing is expected to be more prevalent than in ΛCDM. If and when high-z microlensing signatures are detected, further support for the PBH-ΛCDM model is likely. Though a calculation is hard to set up, given the granularity of PBHs in every DM halo in this model, we can expect many more gravitational wave coalescence events arising from three-body interactions (see, e.g., Gow et al. 2020; García-Bellido 2017).

To conclude, in this paper we show that a broad mass spectrum for PBH DM can be easily accommodated in a model of the early universe, producing features comparable to ΛCDM. One of the main limitations of this is work is the widespread lack of observational data to tightly constrain our parameter space at early epochs, times during which this model maximally differs from ΛCDM. In particular, the reionization history of the universe is largely unknown and future instruments like JWST and Square Kilometer Array (SKA) will likely provide a clearer window into this epoch. At the same time, CIB and CXB fluctuation measurements are currently limited to angular scales of a few tens of arcminutes, making them very susceptible to low counts statistics and sample variance. Future wide field surveys like those planned with the EUCLID, Nancy Grace Roman Wide Field InfraRed Survey Telescope (WFIRST-Roman), eROSITA, and ATHENA missions, will finally permit the measurement of the fluctuations on several degree scale with high accuracy of a few percent. The lack of knowledge of the very high-z SFRD is likely filled in soon by the forthcoming launch of JWST. Deep JWST data will provide a brand new window and allow us to explore star formation and early growth of AGN up to z ∼ 15. The detection of these high-z sources stands to inform our understanding of not just f_PBH but also the efficiency of star formation and the escape fraction at these early epochs, as these are the directly observable signatures predicted by modeling work here.

Finally, as mentioned above, PBH mergers are expected to be very frequent in this scenario and forthcoming LIGO runs and planned LISA, PTA, and third-generation ground-based gravitational wave detectors will help in further testing this alternative DM scenario on this front. While in some respects this DM model might appear more baroque than standard ΛCDM as it involves multiple-mass components—in contrast to the single mass particle CDM—the explanatory power of the scenario makes it compelling. For instance, this is the only model that successfully accounts for the natural formation of early BH seeds; the existence of a correlation between properties of the central SMBH and the host galaxy and underlying DM halo and explains the CXB-CIB cross-correlation and its excess while being consistent with all the multiwavelength observations at z < 6. Other alternate DM models, like those that include dissipative DM with self-interactions that lead to core-collapse (Self Interacting Dark Matter (SIDM) variants) are among non-CDM proposals for BH seed formation (see, for instance, the recent proposal by Xiao et al. 2021). In the absence of any detection of the putative DM particle after several decades of direct and indirect experimental searches, PBH DM offers an economical scenario that is well motivated by physics and is one that couples early universe physics with phenomena on cosmological scales in the late universe.

We thank the anonymous referee for the timely reading and the useful comments that really improved this paper. N.C. acknowledges support from the Chandra-SAO Grant TM9-20008X and kindly acknowledges the University of Miami College of Arts and Science for support. N.C. kindly acknowledges the Cantera, Moscetti, and Larrea families for support. N.C. thanks Kari Helgason for valuable help in developing the model and Fabio Vito for providing extrapolations of the luminosity function. N.C. and G.H. thank INAF-OAS Bologna for kind hospitality in Summer 2021 and for providing convenient office space during the preparation of this paper. The authors thank Sasha Kashlinsky for insightful discussions and for providing the CXB fluctuation estimates. N.C. and G.H. kindly acknowledge the LIBRAE team for fruitful discussions. We also thank Juan García-Bellido and Fernando Atrio-Barandela for helpful discussions. P.N. acknowledges the Black Hole Initiative (BHI) at Harvard University, which is supported by grants from the Gordon and Betty Moore Foundation and the John Templeton Foundation for acting as hosts. P.N. is grateful for the Zoom platform for enabling continued collaborative work with N.C. and G.H. during this past difficult pandemic year. G.H. thanks the City of Remscheid, Germany, for the award of the Röntgen Medal and their hospitality during the days when this manuscript was submitted.

In this appendix, we show the posterior distributions derived by our MCMC chains. In Figure A1 we plot the confidence contours for the 16%, 50% and 84% quartiles. Each panel represents each combination of the possible nine IMF and enrichment models.

In this section we provide in Figures B1–B7 the plots with results of the IMF-3B and IMF-3C model realizations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size