Distinguishing the impact and signature of black holes from different origins in early cosmic history

We use semi-analytical models to study the effects of primordial black hole (PBH) accretion on the cosmic radiation background during the epoch of reionization ($z\gtrsim 6$). We consider PBHs floating in the intergalactic medium (IGM), and located inside haloes, where star formation can occur. For stars with a mass $\gtrsim 25 \rm\ M_{\odot}$, formed in suitable host haloes, we assume they quickly burn out and form stellar remnant black holes (SRBHs). Since SRBHs also accrete material from their surroundings, we consider them to have similar radiation feedback as PBHs in the halo environment. To estimate the background radiation level more accurately, we take into account the impact of PBHs on structure formation, allowing an improved modeling of the halo mass function. We consider the radiation feedback from a broad suite of black holes: PBHs, SRBHs, high-mass X-ray binaries (HMXBs), and supermassive black holes (SMBHs). We find that at $z\gtrsim 30$, the radiation background energy density is generated by PBHs accreting in the IGM, whereas at lower redshifts, the accretion feedback power from haloes dominates. We also analyze the total power density by modeling the accretion spectral energy distribution (SED), and break it down into select wavebands. In the UV band, we find that for $f_{\rm PBH} \lesssim 10^{-3}$, the H-ionizing and Lyman-$\alpha$ fluxes from PBH accretion feedback do not violate existing constraints on the timing of reionization, and on the effective Wouthuysen-Field coupling of the 21-cm spin temperature of neutral hydrogen to the kinetic temperature of the IGM. However, in the X-ray band, with the same abundance, PBHs contribute significantly and could account for the unresolved part of the cosmic X-ray background.


INTRODUCTION
Dark matter (DM) has been a long-standing mystery in astrophysical cosmology (e.g.Feng 2010).Although it comprises a significant portion of the cosmic mass-energy density, we still have no knowledge of its micro-physical nature and origin.Whereas most explanations invoke extensions of the standard model of particle physics (e.g.Bertone et al. 2012), an alternative DM scenario is based on primordial black holes (PBHs), formed via the collapse of overdense regions in the early Universe, prior to matter-radiation equality (Zel'dovich & Novikov 1967;Hawking 1971;Sureda et al. 2021).The recent detection of gravitational wave (GW) events (Abbott et al. 2016;Abbott et al. 2020) from merging binary black holes (BBHs) has triggered a new wave of interest in PBHs.In particular, the detection of BBH merger events with total mass ∼150 M ⊙ (Abbott et al. 2020) challenges current stellar evolution models (e.g.Liu & Bromm 2020b), raising the possibility of a primordial origin (e.g.Bird et al. 2016;Kashlinsky 2016;Clesse & García-Bellido 2017;De Luca et al. 2020c,d).
However, it is unlikely that PBHs within the mass range  PBH ∼ 1 − 100 M ⊙ , typical for stars and their remnants, provide all of the ★ E-mail: szhangphys@utexas.educold dark matter (CDM) measured in the Universe, based on multiple empirical constraints (e.g.Carr & Kühnel 2020).Their contribution can be expressed as  PBH = Ω PBH /Ω DM , the ratio of mass fraction in PBH DM to the total CDM.A prominent observational signature is the GW signal generated from mergers of BBHs.The statistical study of BH merger rates captured by the LIGO-Virgo-KAGRA network (LVK; Ali-Haïmoud et al. 2017;Abbott et al. 2021) constrains the abundance of PBHs to within  PBH ≲ 0.1 − 10 −3 .Furthermore, the non-detection of the stochastic gravitational wave background (SGWB) by LVK gives a similar but weaker constraint of  PBH ≲ 0.1 (Vaskonen & Veermäe 2020;Bagui & Clesse 2022).The recent evidence for a SGWB from pulsar timing array data (e.g.Agazie et al. 2023) implicates a cosmological population of supermassive black hole (SMBH) binaries, and provides promising constraints on supermassive PBHs (  PBH ≲ 0.01) as well (e.g.Acuña & Tseng 2023;Gouttenoire et al. 2023).
Even though stellar-mass and supermassive PBHs only make up a small fraction of dark matter, they can play important roles in early cosmic history.One of the main physical processes that could affect cosmic evolution is PBH accretion, as the focus of this paper.When PBHs travel through gaseous environments, by the pull of their gravitational well, they will attract and accrete the surrounding gas.In this study, for most cases, we will assume steady accretion of gas onto BHs, as found in Bondi (1952).The accreted gas will lose potential energy and emit the heat generated during this process.The resulting feedback radiation will affect the surrounding medium that PBHs are located in (e.g.Takhistov et al. 2022), and will change the overall cosmic thermal history (e.g.Mena et al. 2019).More specifically, multiple studies show that the energy feedback from this accretion will contribute to cosmic radiation backgrounds at different wavebands (e.g.Ali-Haïmoud & Kamionkowski 2017;Inoue & Kusenko 2017;De Luca et al. 2020a, 2022b;Cappelluti et al. 2022;Ziparo et al. 2022).The change in cosmic microwave background (CMB) power spectra derived from PBH accretion limits their abundance, where stellar-mass PBHs are found to comprise no more than a small fraction of DM,  PBH ≲ 0.1 − 10 −3 (e.g.Poulin et al. 2017;Serpico et al. 2020), similar to the constraints derived from GW observations.Further constraints from accretion onto PBHs are summarized in Carr et al. (2021).
Moreover, with the recently launched JWST, observing in the nearinfrared (NIR) electromagnetic (EM) wave band [0.7 − 5 m], astronomers start to unveil an increasing number of galaxies at  ≳ 12 (e.g.Finkelstein et al. 2022), discover SMBHs during their early formation stages (e.g.Natarajan et al. 2017;Larson et al. 2023), and aim to study anisotropies in the cosmic near-infrared background (CIB; Kashlinsky et al. 2015;Helgason et al. 2016).Some of the galaxies, if confirmed by spectroscopic follow-up, are so massive that they could challenge the standard ΛCDM model of cosmological structure formation (Labbé et al. 2023;Boylan-Kolchin 2023).PBHs can potentially address these latest observations in a broader context with their effects on the matter power spectrum and the subsequent structure formation history, as well as on the halo mass distribution (e.g.Afshordi et al. 2003;Carr & Silk 2018;Kashlinsky 2021;Atrio-Barandela 2022;Liu & Bromm 2022a).This change to the halo mass function will also affect the total accretion feedback from PBHs located in haloes, and consequently the cosmic radiation background generated by the aggregate emission (Kashlinsky 2016;Cappelluti et al. 2022;Liu et al. 2022).
In this paper, we will assess the impact of PBHs in the mass range  PBH ∼ 1 − 100 M ⊙ on early cosmic history, considering all BHs generated from different mechanisms as possible power sources contributing to the cosmological background radiation.We will consider the redshift range  ∼ 100 − 6, starting from an almost homogeneous matter distribution to the end of reionization.Overall, we will divide our treatment into PBHs floating in the intergalactic medium (IGM), and those residing in DM haloes.Calculating the resulting feedback energy from the accretion of the surrounding gas, in the former case we consider the effects of Hubble flow and Compton scattering, which will effectively slow down the accretion process (see Ricotti 2007;Ricotti et al. 2008).For the latter case, for simplicity, we assume that PBHs sit near the center of the halo following an isothermal density profile, as previously simulated in Liu et al. (2022).In addition, we consider astrophysical BHs generated from the death of massive stars (≳ 30 M ⊙ ) in virialized haloes (e.g.Husain et al. 2021), as well as SMBHs in more massive haloes with  h ≳ 10 8 M ⊙ (e.g.Jeon et al. 2022).This broad suite of accreting BHs in the early Universe is likely contributing to the sources for the unresolved part of the fluctuations in the cosmic X-ray and infrared backgrounds (CXB/CIB), as summarized in the review by Kashlinsky et al. (2018).
A similar topic has been discussed in recent works (Cappelluti et al. 2022;Ziparo et al. 2022), where the radiation background from PBH accretion has been evaluated.However, different from the assumptions of those studies, we will focus on the radiation background contribution from all BH sources, including stellar-remnant black holes (SRBHs), such as high-mass X-ray binaries (HMXBs).Therefore, we will integrate over haloes with a broad mass range of ∼ 10 5 −10 13 M ⊙ , whereas Ziparo et al. (2022) only considered minihaloes without star formation (≲ 10 7 M ⊙ ).Also, we will consider the case of enhanced accretion by PBH-seeded DM haloes within the IGM, whereas Cappelluti et al. (2022) only considered PBH accretion in star-forming haloes.With our new implementations, we will reach a broader understanding of the PBH-feedback-generated radiation backgrounds.
The structure of our paper is as follows.Section 2 will address the impact of PBHs on structure formation and introduce the different sources of black holes, to be considered in the following treatment.In Section 3, we will discuss the accretion physics in the IGM and DM haloes separately.Subsequently, we derive the resulting generation of radiation and heat, assessing the corresponding observational signatures in Section 4. We summarize and offer conclusions in Section 5. Throughout this paper, we adopt a PBH-DM ΛCDM cosmology with parameters: Ω m = 0.3153, Ω b = 0.04930,  s = 0.96, and ℎ = 0.6736 (Planck Collaboration et al. 2020).The density fraction of CDM is then given by: Ω

BLACK HOLES IN THE EARLY UNIVERSE
In this section, we will address the overall cosmic evolution at high redshifts in the presence of PBHs, as summarized in Figure 1, where we display the entire suite of BHs from different origins.First, in Section 2.1, we briefly summarize previous results on how cosmological structure formation is modified by PBH DM, including changes to the halo mass function, as compared to standard ΛCDM.Then, in Section 2.2, we touch on select PBH formation mechanisms and the resulting mass distribution.We similarly discuss relevant aspects for BHs of astrophysical origins in Section 2.3.

Cosmological Structure Formation
Following the general formalism introduced in previous work (Afshordi et al. 2003;Inman & Ali-Haïmoud 2019), we consider density perturbations as a function of redshift , consisting of adiabatic and isocurvature terms,  ad () and  iso (), respectively.Adiabatic perturbations uniformly affect all mass-energy components, such that  ad () applies to all forms of DM.In this context, the isocurvature term derives from the discreteness of individual PBHs (Kashlinsky 2016;Desjacques & Riotto 2018;De Luca et al. 2020d), and can be written as Here,   is the Dirac delta function, x  marks the position of the -th PBH, and nPBH is the average PBH number density, which depends on the average PBH mass.This formalism can be extended to a broader mass spectrum, with a characteristic mass scale of  c , as stated in, e.g., Carr (1975); Choptuik (1993).In general, we thus write the average density as Given the limits on the abundance of PBHs, implying that they do not compose the totality of DM, we assume that the remainder of (cold) DM is provided by particle DM (PDM).We thus consider an additional term imprinted by the PBH isocurvature perturbations on In certain cases, some of the stars form in binaries, where the more massive component becomes a black hole.Through accretion from the companion star, the black hole emits strong X-rays, giving rise to low/high mass X-ray binaries (L/HMXBs), depending on the mass of the BH (see the lower right panel).As some of the haloes have grown increasingly massive, in excess of ≳ 10 8 M ⊙ , the black hole in the center of the galactic bulge starts to accrete at nearly Eddington rate to grow into supermassive black holes (SMBHs), which exert vigorous radiative feedback during this process.This more complex stage in cosmic evolution, with the full panoply of BH sources, is summarized in the lower left panel.
the PDM, such that the combined evolution of the density perturbations in both DM (PDM) and PBHs is given as follows: The evolution of adiabatic perturbations with respect to  is expressed via  ad () =  ad () 0 ad .Here,  0 ad and  0 iso correspond to the primordial adiabatic and PBH isocurvature perturbations, whereas  ad () and  iso () are the linear transfer functions for the two modes, both normalized to 1 as  → ∞.The perturbation in PDM induced by PBHs is  iso () = [ iso () − 1]  PBH  0 iso , and we can write the transfer function for PBH-induced isocurvature fluctuations as  PBH iso () = 1 + [ iso () − 1]  PBH .We note that the transfer function is in turn related to the structure growth factor,  () =  ().To a good approximation, we can use the results from solving the cosmic fluid equations with embedded point masses (Meszaros 1974), and adopt the following analytical fit for the growth factor (Inman & Ali-Haïmoud 2019): where  eq ∼ 3400 represents the redshift at matter-radiation equality.
Given the two-point correlation function for the density perturbations from PBH DM: where  PBH () represents the reduced correlation function as a function of separation distance  = |x| (e.g.Chisholm 2006;Desjacques & Riotto 2018;Dizgah et al. 2019;De Luca et al. 2022a), we can calculate the power spectrum by taking the spatial Fourier transform of the correlation function of all perturbation modes.Specifically, the isocurvature term  iso () from PBHs approximates to  2 PBH / nPBH .We can then express the total power spectrum1 at  = 0 as follows: Initially, the adiabatic and iso-curvature modes are not correlated.However, at a later stage ( ∼ 10), we need to consider the correlation between the isocurvature and adiabatic perturbations since PBHs are more abundant in overdense regions on larger scales, but disrupt DM structure around them on smaller ones.Here, similar to the prescription in Liu et al. (2022), we introduce a mode mixing term governed by the transfer function  mix .Considering the inter-PBH separation scale, , we use the heuristic formula calibrated to simulation results in Liu et al. (2022): where  0 =  ( = 0) is the growth factor for the isocurvature mode, and  ad,0 =  ad ( = 0)/ ad ( =  eq ) for the adiabatic mode 2 , both evaluated at redshift 0.
From the total power spectrum, we calculate the density contrast as a function of the window mass : where  TH ( M ) is the top-hat Fourier transform of the spherical window function, and  M corresponds to the comoving radius associated with the mass window .We note that for non-linear scales, where () ≳ 1, this mass corresponds to a virialized halo structure.Therefore, in the following, we will assume  ℎ =  for the halo mass.Employing the Press-Schechter formalism for the halo distribution (Press & Schechter 1974) 3 , the number density of virialized structures (halo mass function) between  h and  h +  h is: where ρ is the average comoving matter density, and the critical overdensity factor,  c (), as a function of redshift is given by Integrating the halo number density over the mass range of interest  Figure 2. Collapse fraction in haloes larger than 10 5 (pentagon), 10 6 (triangle), and 10 8 M ⊙ (circle), vs. redshift .Here, we consider the changes in the formation of structures in the presence of PBHs (red, cyan and yellow, compared to the standard ΛCDM case (blue).
and dividing by the average mass density, the collapse fraction of the cosmic matter is given by and is plotted in Figure 2.
As can be seen, in the presence of PBHs, more minihaloes ( h ∼ 10 5 M ⊙ ) form than in the standard ΛCDM case.Also, the abundance of minihaloes will depend on the characteristic PBH mass,  c , but more strongly on the PBH abundance.This behaviour can be explained by the shift of total power spectrum  tot in the presence of PBHs by the isocurvature term  2 PBH / nPBH in Equ.(6), which changes the density contrast subsequently.In addition, we have examined the effect of the mode mixing term in Equ. ( 7), finding that it does not significantly change the collapse fraction within the PBH mass scale considered here.However, as we increase the mass scale to ≳ 100 M ⊙ , minihaloes of mass ≲ 10 6 M ⊙ would be less abundant, if mode mixing were excluded.Going forward, the collapse fraction allows us to calculate the PBH abundance in the IGM, whereas the halo mass function enables us to evaluate the prevalence of PBHs in DM haloes.

Primordial Black Holes
Several previous works provide reviews of the PBH formation mechanisms and abundance constraints from phenomenological studies (e.g.Carr & Kühnel 2020;Carr et al. 2021;Escrivà et al. 2023).In this paper, we consider PBHs with different mass distributions and study how these changes affect the radiation background output.To begin, the PBH mass distribution ( PBH ) with respect to PBH mass  PBH is related to the differential number density as follows: Here, ( PBH ) is the mass probability distribution function, normalized to 1 when integrating over ln  PBH , and it reflects the specific PBH formation mechanism.In the most common scenario, PBHs follow a monochromatic mass distribution when all PBHs form from the collapse of overdensities during the radiation-dominated era (e.g.Carr 1975).The mass function then takes the simple form: where  c is the characteristic PBH mass scale.
Another scenario suggests that PBHs form in smooth peaks of the inflation power spectrum, generating an initial mass function with a lognormal shape (Dolgov & Silk 1993).The functional form is given by: where  denotes the width of the PBH distribution.Constraints on the parameter space (  PBH , ,  c ) for the lognormal distribution have been calculated in Kühnel & Freese (2017), specifically concluding that 0 <  ≤ 1.We note that as  → 0, the lognormal distribution reduces to the monochromatic distribution.We also note that many other shapes of PBH mass functions are discussed in the literature, like the power-law or critical collapse functions (Carr et al. 2021).Here, we mostly restrict our attention to the monochromatic mass function, but our calculation also applies to lognormal distributions with a narrow dispersion.A comprehensive investigation with respect to different shapes of the initial PBH mass function will require another study.

Supermassive Black Holes
We will consider SMBHs as powerful sources of radiation in difference from other astrophysical black holes.For a comprehensive discussion of the formation and impact of SMBHs, we refer the reader to recent reviews (e.g.Smith & Bromm 2019;Inayoshi et al. 2020).Our focus here is on the mass relationship between the SMBH and its host halo.For the sake of our analysis, we operate under the optimistic 4 assumption that SMBHs are located at the center of all 4 Recent JWST observations have discovered a population of active galactic nuclei (AGN) at  ≳ 5 with a much higher abundance than previously expected, which implies a surprising ubiquity of SMBHs in high- galaxies (e.g.Furtak et al. 2023;Greene et al. 2024;Harikane et al. 2023b;Kocevski et al. 2023;Maiolino et al. 2023;Matthee et al. 2023).Motivated by this discovery, we assume that all haloes above a certain mass threshold host SMBHs.
massive haloes with  h ≳ 10 8 M ⊙ .This assumption allows us to directly derive the SMBH mass function from the halo mass function.Specifically, the mass of the central black hole is related to that of its host halo, as follows: where the BH formation efficiency, ( h , ), is a function of halo mass and redshift.In Jeon et al. (2022), it is calculated using a global energy balance argument as: where  is the fraction of the BH rest mass energy deposited into the host halo gas as heat.However, this idealized formalism scales as  ∝ (1 + ), and it does not converge to 0 at high redshifts ( ≳ 15), where SMBHs begin to form.Therefore, we also adopt the fitted value of ( h , ) provided by the UNIVERSE MACHINE, which is based on observational data (Zhang et al. 2023).This is discussed in more detail in Appendix A.

Isolated Black Holes from Stellar Remnants
We also consider stellar-remnant black holes (SRBHs) as the eventual product of the steller evolution.The resulting abundance of SRBHs hinges on modeling of both the underlying stellar evolution and starformation history.To proceed, we need to apply constraints on BHs derived from their progenitor stars.The eventual fate of a single star with an initial mass  i can be summarized as follows (Heger et al. 2003): Here,  ★,c is the critical mass, before mass loss, above which a single star will eventually produce a BH, with a value of approximately 25 M ⊙ .According to Heger et al. (2003), the products of stellar evolution also depend on the metallicity of the stars.For extremely low metallicity stars, such as Population III (Pop III) stars, the summary above accurately captures the fate of dying stars.However, for higher metallicity Pop I/II stars with  ≳ 10 −4 Z ⊙ , where mass loss is expected to be severe, the PISN stage of evolution may never occur, and is replaced by a BH remnant fate.Therefore, we assume that all Pop I/II stars with an initial mass larger than the critical mass will evolve into BHs.
The literature offers various forms for the initial mass function (IMF) for different stellar populations, often characterized by a single or broken power-law (e.g.Chabrier 2003).Since our primary focus here is on the more massive stars that have the potential to form BHs, we adopt a simple power law for the stellar IMF: where  is a normalization constant.For simplicity, we assume that any BH with a progenitor mass larger than 260 M ⊙ contributes to the SMBH mass budget, and is thus implicitly accounted for in the SMBH tally.For simplicity, we assume a universal lower cutoff for the stellar mass, independent of metallicity, such that  i,min = 0.1 M ⊙ and  i,max = 260 M ⊙ represent the minimum and maximum stellar mass for our idealized IMF, respectively.Integrating, we obtain the total stellar mass: The derivative of the initial-to-remnant (SRBH) mass relation is evaluated by interpolating figure 2 from Husain et al. ( 2021).In the case of PISNe, for stars with an initial mass larger than 140 M ⊙ or smaller than the critical mass  ★,c , this function will return 0.
Overall, the total stellar BH mass in a halo of mass  h at a given redshift can be estimated by integrating over the halo star formation history (SFH), as follows: In the above expression,  loss (| i ≥  ★,c ) denotes the mass fraction of the stellar population contained in stars with initial masses exceeding the critical mass  ★,c that have lifetimes shorter than .This fraction is a function of the time difference between the current redshift  and the redshift  ′ at which the stars initially form.Since there is still uncertainty regarding the onset of first star formation (e.g.see fig. 2 in Klessen & Glover 2023), for definiteness, we assume that Pop III star formation begins at  ini = 50.We note that the lifetime of the massive stars (≳ 10 M ⊙ ) capable of forming BHs is less than ∼ O (10 Myr).Therefore, the fraction  loss (| i >  ★,c ) will reach the maximum at  ≳ 10 Myr.Since in most cases, the cosmic time difference between the moment of star formation ( ′ ) and a given reference redshift () will be ≫ 10 Myr, we remove the time dependence and approximate  loss with its final maximum value for simplicity5 : In Equation ( 20),  rem signifies the fraction of the progenitor star's mass converted into (astrophysical) BH remnants, given by the ratio of the final BH mass,  SRBH , to the initial stellar mass,  i .Upon integrating over the relevant mass interval, we arrive at the remnant mass fraction, as follows: To model the star-formation history in a halo at redshift , we separately consider the formation history of Pop III and Pop I/II stars as a function of redshift, represented as SFH( h , ) = SFH III ( h , ) + SFH I/II ( h , ).In a first order approximation, we can express this as the product of the comoving volume corresponding to a given halo, and the comoving cosmic star formation rate density (SFRD): Here,   ( h ()) = 8  ℎ ()/3 2 0 represents the comoving volume corresponding to the halo mass  h (), and the subscript  refers to the specific population of stars.
We note that the SFRD is an average over cosmological volumes, and is thus not accurately representing the conditions in a given halo.Star formation is dominated by Pop I/II stars at redshift  ≲ 10, a phenomenon well-constrained by available data (for a recent review, see Klessen & Glover 2023).Accordingly, we employ the fit from the Universe Machine (Behroozi et al. 2013) for the formation history of Pop I/II stars, denoted as SFH I/II ( h , ).The parametric fit, along with the SMBH masses, are summarized in Appendix A. However, the formation rate at  ≳ 10 remains largely uncertain due to the scarcity of available data, and the fit from the Universe Machine predicts almost no star-formation within minihaloes.As such, we use the fit for Pop III SFRD based on a cosmological simulation presented in Liu & Bromm (2020a), rephrased for the typical redshift  ∼ 20 of Pop III star formation: Evaluating Equ.(20-24), we arrive at Figure 3, which summarizes the total mass of SRBHs as a function of redshift and halo mass.We observe that, generally, as the mass of the dark matter (DM) halo increases, more astrophysical black holes (SRBHs) are produced.Moreover, the total BH mass approximately scales as ∼ 10 −4  h at  ≃ 6.
On the other hand, we also consider two separate fits, namely Madau & Dickinson (2014) and Harikane et al. (2022), for the total star formation rate density (SFRD) as a first-order approximation to the star formation history (SFH) of Pop I/II stars in a galaxy using Equ.( 23).Recent studies have found these two fits to be reasonable bounds for the observed SFRD of galaxies at  = 9 − 16 (Harikane et al. 2023a), with the work by Madau & Dickinson (2014) serving as the upper bound.Therefore, at  ≳ 10, we can express the upper bound and extrapolate it to even higher redshifts as follows: (25)

X-ray Binaries
One important source of radiation feedback arises from the remnants of massive stars that form in binaries.Specifically, we focus on highmass X-ray binaries (HMXBs, reviewed by e.g.Sana et al. 2013;Clark et al. 2023).Therefore, we apply this factor of 0.5 to Pop I/II stars.We acknowledge that not all massive binaries will evolve into HMXBs due to many factors, such as common envelope evolution and metallicity (e.g.Power et al. 2009;Mirabel et al. 2011).Regarding those massive binary stars, for simplicity, we assume that all of them experience stages of X-ray activity and use a simple phenomenological model to calculate their radiation, thus providing an upper limit for the HMXB production and radiation from the corresponding accretion.Specifically, we consider a duty cycle during which HMXBs accrete mass and strongly radiate.To bracket the uncertainties, we assign a value of 0.1 − 1 to represent this duty cycle, which has a weak correlation with HMXB luminosity, according to observations (e.g.Sidoli & Paizis 2018).Furthermore, we assign a lifetime of  acc ∼ 1 − 10 Myr to all HMXBs to allow them to accrete a substantial amount of mass (∼ O (M ⊙ )) from their companions (Mirabel et al. 2011).We estimate the total mass of active HMXBs in a given halo of mass  ℎ as follows: Here,  HMXB refers to the aforementioned fraction of massive stars that form in binaries.The second equality is valid since  acc is much smaller than the cosmic time .We note that a comprehensive assessment of the energy feedback from HMXBs would require a dedicated study with binary population synthesis (see e.g.Fragos et al. 2013a,b;Sartorio et al. 2023;Liu et al. 2024), beyond the scope of the current paper.
With the prescriptions above and integrating over the halo mass function, we can derive the average BH density for different species and redshifts.An example is shown in Figure 4 for a monochromatic PBH mass distribution with  c = 10 M ⊙ , and a mass fraction in DM of  PBH = 10 −3 .We calculated the total mass density of PBHs to be ∼ 10 8 M ⊙ Mpc −3 .We conclude that PBHs will dominate the BH mass budget throughout the redshift range considered here, unless  PBH ≲ 10 −4 .However, as more stars form and die, around  ∼ 6, the mass density in SRBHs (including HMXBs) becomes non-negligible in dark matter haloes, at about ∼ 10 −2 times the mass density of the PBHs.Moreover, we notice that the SRBH fraction in the form of active X-ray binaries is comparable to the abundance of SMBHs, both reaching a mass density level of around ∼ 10 3 M ⊙ Mpc −3 , generally consistent with the result from Jeon et al. (2022).

ACCRETION PHYSICS
In this section, we briefly describe the BH accretion physics as we aim to study the corresponding energy feedback.We consider BHs of different origins, located either in the intergalactic medium (IGM) (Section 3.1) or in haloes (Section 3.2).
In general, regarding the accretion rate of a BH, we adopt a modified version of the Bondi-Hoyle steady accretion formula (Tremmel et al. 2017): Here,  gas =  p represents the average gas density of the environment, given by the product of the average molecular weight  = 1.22, the gas number density , and the proton mass  p .Furthermore,  s represents the sound speed determined by the gas temperature, and  rel is the relative velocity of the BH with respect to the gas.Thus, the effective velocity is given by  eff = √︃  2 s +  2 rel .We also note that for binary BHs the accretion rate could be modified compared to single BHs with the same mass.However, the change in the magnitude of the accretion rate is still under debate, as some simulations identify suppressed accretion rates because of the gravitational torque exerted by binaries (e.g.MacFadyen & Milosavljević 2008;Miranda et al. 2017), whereas others find that binary accretion rates are not significantly modified (e.g.Farris et al. 2014;Shi & Krolik 2015).Therefore, for simplicity, we assume a single BH accretion rate for all BHs.
Next, we introduce a dimensionless accretion factor: where  Edd represents the Eddington accretion rate given by: Here,  0 refers to the radiative efficiency.

Black Holes in the IGM
Since there is no process that allows for star formation in the IGM, we here only consider the case of PBH accretion.For the convenience of the reader, we provide a brief summary of the relevant analytical formalism, as developed in previous studies.For PBH accretion in the IGM, we adopt the self-similar model that incorporates spherical Bondi accretion with the damping effect of Compton drag and Hubble flow at high redshifts (Ricotti 2007;Ricotti et al. 2008;Ali-Haïmoud & Kamionkowski 2017).For isolated PBH accretion in the IGM, the accretion factor, normalized to typical values, can be expressed as follows 6 : Here,  represents a scaling constant that represents the damping effect, numerically fitted as follows: The effective viscosity β is in turn determined by the Compton drag and Hubble flow, given by: Here,  e represents the ionization fraction as a function of redshift.For simplicity, we approximate  e as a step function, in line with the CMB constraint from Planck Collaboration et al. (2020).For 6 ≲  ≲ 1000, which extends to epochs long after recombination, we take its value to be a constant7 , specifically  e ∼ 10 −3 .For  ≲ 6, after reionization, this factor is set to 1.In all cases considered, we have β ≪ 1.Then, the scaling constant can be approximated as  ≃ exp(3/2)/2 ≃ 2.24.
For the gas effective velocity  eff in the IGM, we adopt the analytical approximation from Ricotti et al. (2008) and Ali-Haïmoud & Kamionkowski (2017) by averaging over the Maxwell velocity distribution for PBHs: Here, M = ⟨ rel ⟩ / s represents the Mach number, and ⟨ rel ⟩ is the root-mean-square of the relative velocity between dark matter (DM) and baryons, approximated as8 : ⟨ rel ⟩ ≈ min 1, /10 3 × 30 km s −1 . (34) The gas sound speed,  s , is calculated from the CMB temperature, which is a function of redshift.We use the parametric fit given in De Luca et al. (2020b): Here,  = 1.72 is a fitting parameter, and  dec ≃ 130 the decoupling redshift for baryons from the radiation field.Combining Equ. ( 29) and ( 30), we can approximate the mass growth rate as: Using the above equations to calculate the accretion rate and integrating over cosmic time, we find that for PBH masses of interest,  PBH ≲ 100 M ⊙ , the change in mass is  PBH / PBH ≪ 1, which is negligible.Therefore, to make conservative estimates and simplify our calculations, for PBH accretion in the IGM, we assume a mass function () that is independent of redshift.However, we also note that, although yet to be observed, the accretion rate of PBHs, , could exceed the Eddington limit, if DM haloes are seeded by more massive PBHs (with ≫ 100 M ⊙ ) in isolation (e.g.Ricotti 2007;Mack et al. 2007;Ricotti et al. 2008;De Luca et al. 2020a).In this scenario, the PBH can accrete a significant amount of mass by  ≲ 6, thus significantly changing the shape of the PBH mass function.To fully address this effect, a dedicated simulation study would be required following up on the work by Liu & Bromm (2022b).

Black Holes in haloes
For the case of accretion within virialized haloes, we encounter more complexity compared to the previous case, since there will be BHs of various origins within a given DM halo.To simplify the calculations, we assume that both SRBHs and PBHs undergo Bondi accretion: Here, for a BH at a given radius  relative to the center of the halo, () is the local gas number density.For the case of monochromatic PBHs, we use the simulated results from Liu et al. (2022), suggesting a nearly isothermal gas density profile 9 : () ∝  −2.2 .Additionally, we assume the same scaling for the PBH number density.For SRBHs, we assume a similar density profile for simplicity, since the formation of SRBHs should approximately track the gas density profile, which in turn has a similar shape as the DM density profile in this regime, according to previous simulation results (see fig. 7 in Liu et al. 2022).
Further simulation studies will be needed to more realistically track the location of SRBHs vs. PBHs in a halo, including any dynamical segregation effects at later stages.Meanwhile, we use the virial velocity as the effective velocity  eff between BHs and the gas: Here,  vir is the virial radius of the halo, which is related to halo mass and redshift by: where Δ c is the characteristic overdensity of the halo with respect to the cosmic average, and we adopt a fiducial value of 200.We note that for certain cases of BH accretion, Equ.(37) could yield rates larger than the Eddington limit.Conservatively, we enforce  ≲ 1 throughout, and thus do not consider any super-Eddington accretion regime, which we cannot reliably represent with our idealized modeling here.
The dynamics and density distribution of PBHs with an extended mass function can be complex, and its effect on SRBH and SMBH formation would require future simulations.Additionally, SRBHs will also have a non-trivial mass function, adding further complexity to the situation.For now, we can derive order of magnitude estimations.We assume that the extended PBH mass distribution will not significantly change the density profiles of both gas and BHs in the host halo.To calculate the gas number density at different radii, we normalize to the value at the virial radius: Similarly, since no previous simulations have considered the coevolution of PBHs and SRBHs, for simplicity, we apply the same formula to calculate the accretion rate for both BH classes.To further simplify our calculations, we impose additional assumptions, as follows.For more massive SRBHs (with progenitor stellar mass ≳ 260 M ⊙ ), we assume that they gradually sink towards the halo center to become seeds for SMBHs, or get captured by the pre-existing central SMBH.Thus, we cut off the progenitor stellar mass functions for SRBHs at ≳ 260 M ⊙ and restrict their radii to be inside of 0.001 vir .Typically, for a 10 9 M ⊙ halo at  ≃ 10, the virial radius is  vir ∼ 3 kpc.The inner region,  < 0.001 vir , thus has a size ∼ O (pc), and can only contain a SMBH.For PBHs following a lognormal mass function, we numerically integrate from 0.1 − 10 c to cover the plausible mass range.Furthermore, we assume a time independent mass function for both PBHs and SRBHs (excluding HMXBs).We find that any evolution of the BH mass function in a halo is likely small, as we argue below.Since the halo will also capture PBHs during its growth, we can compare the total PBH accretion rate with the capture rate of PBHs by DM halo accretion.A fiducial value for the capture rate is the accretion rate calculated for the Millennium simulation (McBride et al. 2009), modified for a contribution from PBH DM, as expressed in the PBH mass fraction  PBH : This is of the same order of magnitude as the accretion rate obtained by integrating over all PBHs with a characteristic mass of 100 M ⊙ in a 10 12  h halo at  = 0, as given by Equ.( 36)-(40).Therefore, at higher redshifts, we would expect a much larger PBH capture rate than the total accretion rate.For SRBHs, we find that the increase in total SRBH mass is proportional to the star formation rate, i.e. ∼ 0.01 SFH( h , ), as calculated using Equ.( 20), which is also significantly larger than the total accretion rate.In both cases, the newly increased BH population will likely preserve the mass function at formation.In general, however, the overall BH mass function in a given halo will result from the competition between the BHs that already sit in the halo versus newly accreted or born BHs.Although we do not explore such cases here, a more complete treatment should broaden the discussion.A more detailed simulation study of the mass function evolution within a DM halo would be necessary in future work.For now, we can accommodate only minor changes to the overall mass function of SRBHs and PBHs.Finally, for HMXBs and SMBHs, we note that feedback could affect the accretion efficiency in complex ways. 10For simplicity and consistent with the prescription in Jeon et al. (2014), we impose Eddington-limited accretion throughout (  = 1).

Cosmic Radiation Background
In this section, we aim to study the contribution from all BHs, as summarized in Figure 1, to the cosmic radiation background.Here, we focus on cases with monochromatic PBH mass distribution, and defer fully explored results for an extended mass distribution to future work.We begin by splitting the energy density  from BH feedback into two parts based on the BH habitat: where  IGM represents the energy feedback from PBHs accreting in the IGM, and  halo that from all BHs in DM haloes.Based on the DM halo mass, we further decompose the feedback in collapsed structures into contributions from SRBHs (including HMXBs), PBHs, and SMBHs, as follows: Here, several DM halo mass limits are introduced.Specifically,  min () refers to the smallest halo mass that could host a PBH, given by  min = max(  Ω m /  PBH Ω DM ,  vir ), where  vir refers to the smallest mass that can virialize at a given redshift into a halo structure with a virial temperature ≳ 100 K.We also notice some other pertinent mass scales, such as ∼  c (1 +  eq )/(1 + ), describing a halo seeded by a single PBH (see e.g.Ricotti 2007;Ricotti et al. 2008).This effect will further increase the radiation output since more minihaloes are present at high redshift.In our modeling, we derive a conservative estimate with the current choice of mass threshold and address the PBH-seeded halo case in Appendix C. We defer a more complete understanding of halo accretion at high redshifts to future study with simulations. H 2 −cool is the smallest mass of a DM halo capable of forming Pop III stars via H 2 cooling, and it is fitted by (Trenti & Stiavelli 2009). max () is the maximum allowed mass of DM haloes at redshift , given by the progenitor of ∼ 10 16 M ⊙ haloes at  = 0, following their growth history 11 .Haloes with mass larger than  max () are excluded due to low number density.The growth trajectories of all haloes are modeled by Behroozi et al. (2013) and Behroozi & Silk (2018) (see Appendix A).For PBHs in haloes, since this halo growth model might not be applicable at high redshift, we set a uniform upper limit of 10 13 M ⊙ for numeric approximation.
To calculate each term, we need to revisit the accretion physics and study the radiative efficiency based on the accretion rate.As a good estimation, we adopt the sub-grid model from Negri & Volonteri (2017) to capture the transition from a radiatively efficient thin disk to an optically thick and advection-dominated accretion flow (ADAF) case.The relationship between BH accretion rate  and bolometric luminosity is given by: where  EM is the electromagnetic efficiency,  = 100, and  0 = 0.057 the radiative efficiency for optically thin accretion.Note that  EM reduces to  0 for higher accretion rates,  ≫ 0.01.
To evaluate the total luminosity from BH accretion in a halo of mass  h at redshift , we sum over the mass distributions of BHs, and normalize with respect to the BH mass function.In a general form, the total luminosity from BHs of a single origin is given by: where Δ where  ( h , ) is the comoving number density of haloes per unit halo mass, given by Equ.( 9).The factor of (1 + ) 3 is applied to derive the physical number density at a given .We note that, if the maximum halo mass  max is smaller than the lower limit  min at some redshift , there is no contribution from that particular type of BH.The lower integration limit  min depends on the BH type.
For PBHs located in a halo, For SRBHs, we choose  min =  H 2 −cool , whereas for SMBHs, we adopt  min = 10 8 M ⊙ as the minimum halo mass that can host an AGN.The proper bolometric power density  at redshift  from BHs can be obtained by integrating the power density at a given redshift  ′ from different sources over cosmic time: where  i and  are the initial and final redshift of consideration, respectively.Here, we consider the most distant PBH source emitting at  i ≃ 1000.For the case of SRBHs, we start the integration from  i ≃ 50.
The results are shown in Figure 5 for different PBH mass fractions:  PBH = 10 −3 (upper panel) and  PBH = 10 −5 (lower panel), both assuming a monochromatic mass distribution with  c = 10 M ⊙ .The figure confirms our previous statement that PBHs can only make up a small fraction of dark matter, as otherwise reionization would have occurred too early (see below for further discussion).If we assume a PBH mass fraction of 10 −3 , the feedback power from PBHs in haloes starts to dominate at  ≲ 30.We also note that, since SMBHs and HMXBs are accreting at a much higher rate than SRBHs, their emission power at lower redshift ( ≲ 10) starts to surpass that from SRBHs.In addition, we discuss the energy output from the enhanced accretion rate within PBH-seeded halos in Appendix C.
For context, we compare our predicted bolometric luminosities to the lower limits on the cosmic radiation density required for key cosmological signatures.To render the 21-cm hyperfine structure line of neutral hydrogen detectable in the form of absorption, its internal spin temperature needs to be coupled to the IGM kinetic temperature.An effective way to accomplish this coupling is through the Wouthuysen-Field effect (Wouthuysen 1952;Field 1959), where Lyman- photons scatter with H atoms, thus indirectly modifying the hyperfine level populations (Madau et al. 1997).We can estimate the required Lyman- energy density, in proper units and assuming a flat spectrum, as follows:  Ly ≃ 4    / ≳ 3.5 × 10 −14 [(1 + )/7] erg cm −3 (Ciardi & Madau 2003).Similarly, we can estimate the minimum energy density in H-ionizing UV radiation required to reionize the IGM: Madau et al. 1999).However, to gauge the possible impact of PBH DM on cosmic reionization, we need to go beyond bolometric quantities, towards a more careful spectral modelling for accreting BHs (see below).Furthermore, any reionization constraints have to take into account the escape fraction of UV radiation from the source host haloes.
To obtain the (angle-averaged) radiation intensity at frequency  and a given redshift , we need to calculate the specific emission coefficient,  •, ′ ( ′ ), from sources at higher redshifts  ′ >  with  ′ = (1 +  ′ )/(1 + ).Applying the same procedure as for Equ.(45-46) but replacing the bolometric luminosity with the specific luminosity, i.e.  →   (per BH),  • →  •, (per halo), and  • →  •, , we first sum over the contributions from all BHs in a halo, and .Background intensity at the Lyman limit,  UV = 13.6 eV/ℎ, vs. redshift.Here, we again consider various BH sources (PBHs, SRBHs, HMXBs, and SMBHs), differentiating between halo and IGM PBH contributions.For the sake of simplicity, we assume a monochromatic mass function for PBHs with a characteristic mass of 10 M ⊙ and a mass fraction of 10 −3 .We show the minimum intensity,  UV , required for reionization (Madau et al. 1999), as well as the minimum Lyman- background,   , needed for efficient Wouthuysen-Field (WF) coupling (Ciardi & Madau 2003).For the comparison to the latter, we for simplicity assume   ≃  UV , and for both limits, we indicate uncertainties with a factor-of-ten range (grey shaded areas).We also indicate the timing constraint for effective WF coupling at  ∼ 15 − 20 (purple shaded area), inferred from the EDGES 21cm-cosmology measurement (Mirocha & Furlanetto 2019).Additionally, the brown region denotes the reionization limit at  ∼ 6.
then integrate the halo mass function to obtain the specific emission coefficient.Here, the features of the emerging spectrum from a single BH depend mostly on the normalized accretion rate , determining whether accretion operates in the thin-disk or ADAF regime.We employ the prescription of Takhistov et al. (2022), as summarized in Appendix B, to calculate the specific luminosity   ( • ) for any BH.Finally, integrating all contributions from preceding redshifts, the resulting specific intensity   at  is (e.g.Schauer et al. 2019): We note that both Equ. ( 47) and ( 48) have no absorption term, thus assuming a transparent IGM and maximizing the radiation level from BH sources.Therefore,  i denotes the largest redshift that a given BH source starts to emit light.To calculate the intensity for a given frequency band, we simply integrate   over the corresponding band, and  i denotes the highest redshift below which radiation produced by BHs can reach  in the given band.In the following sub-sections, we will specifically discuss the intensity in the UV and X-ray bands , comparing them with observational constraints.We emphasize that when evaluating the resulting feedback on the IGM, the escape fraction  esc of radiation from the host haloes should be taken into account, depending on the waveband considered.For now, and unless otherwise noted, we set the escape fraction to 1, thus maximizing the radiative output from the accreting BH sources (further discussed below).

UV Intensity
In Figure 6, we present our calculations of the (average) UV intensity from BHs as a function of redshift, assuming a monochromatic mass function for PBHs with a characteristic mass  c = 10 M ⊙ and a mass fraction  PBH = 10 −3 .We evaluate the radiation intensity at  UV = 13.6 eV/ℎ, and assume for simplicity the same value at   , i.e.   ≃  UV , which is sufficiently precise for our purpose here.It is evident that while the contribution from PBHs in the IGM initially dominates, it becomes significantly smaller compared to the radiation from BHs in haloes for  ≲ 30.Our estimation of the UV power from SMBHs is about an order of magnitude lower than the findings in Jeon et al. (2022).This discrepancy arises because of differences in the spectral modeling, employing a more detailed approach here (see Appendix B).Another contributing factor might be our model's prediction of a BH mass that is smaller by roughly a factor of ∼ 10 when compared to observed high-redshift quasars (e.g.Decarli et al. 2018;Wang et al. 2021).
According to a confluence of observations, the IGM was fully reionized by  ∼ 6 (reviewed in Robertson et al. 2010), with early evidence based on the Gunn-Peterson absorption trough in high- quasar spectra (e.g.Becker et al. 2001).When juxtaposing the reionization limit with the radiation background from all BHs, it is clear that PBHs consistently dominate over other BH sources prior to reionization, as long as  PBH ≳ 10 −3 .For our fiducial value of cosmic PBH density (0.1% of all DM), we can see that at  ≲ 15, the UV intensity from PBHs (within haloes) becomes comparable to the level required for reionization, suggesting that they could be a significant source of ionizing photons for the neutral hydrogen in the IGM.Taken at face value, such large UV intensities would allow for reionization much earlier than the canonical redshift of  ∼ 6.However, when calculating the resulting UV intensity, we have assumed that all ionizing photons produced via BH accretion inside a halo can escape from the host (  esc = 1).More realistic treatments, despite the inherent uncertainty, have constrained this factor to be significantly smaller, around 0.1 for typical host systems in the pre-reionization Universe (e.g.Gnedin et al. 2008;Khaire et al. 2016).Given our ballpark calculations, as presented in Figure 6, PBH DM could contribute significantly to the sources of reionization.Depending on the (uncertain) UV escape fraction, this yields another constraint on the PBH DM fraction, roughly giving  PBH ≲ 10 −3 /  esc .
Similarly, we assess whether the PBH contribution to the cosmic Lyman- background can induce Wouthuysen-Field coupling, thus imprinting a global absorption signature in the redshifted 21-cm radiation (e.g.Pritchard & Loeb 2012).Such a signature was reported by the Experiment to Detect the Global Epoch of Reionization Signature (EDGES), suggesting that the spin temperature of neutral hydrogen was coupled via the Wouthuysen-Field effect to the kinetic gas temperature around  ≃ 15 − 20 (Bowman et al. 2018;Mirocha & Furlanetto 2019).This claimed detection remains a topic of vigorous debate, with more recent spectral radiometer measurements that have found no such signal (Singh et al. 2022).Prior research (Mena et al. 2019) modeled the alterations in the kinetic temperature of IGM gas due to PBH heating, constraining the PBH abundance to  PBH ≲ 10 −3 for  c = 10 M ⊙ .Using the same characteristic mass, we find a similar constraint on the PBH abundance, as  PBH > 10 −3 would imply effective Wouthuysen-Field coupling at redshifts larger than the EDGES range (see Fig. 6).We note that in arriving at this conclusion, we have assumed an escape fraction for Lyman- photons of close to unity, which is plausibe for the dust-poor galaxies at cosmic dawn (e.g.Jaacks et al. 2018;Smith et al. 2019)  all BH sources considered here, applying the same set of parameters as used in Figure 6.The black dash-dotted line represents the unresolved Xray background (Cappelluti et al. 2017(Cappelluti et al. , 2022)).The upper panel shows the soft X-ray band with a photon energy in the [0.5 − 2keV] range, while the lower panel depicts the integrated hard X-ray band with photon energies of potent sources for Lyman- background radiation at high redshifts, whereas the other BH sources considered here fall short in producing a significant number of Lyman- photons.For a discussion of the contribution of stellar sources to the cosmic Lyman- radiation field, we refer to the extensive literature on this topic (e.g.Fialkov & Barkana 2014;Schauer et al. 2019;Gessey-Jones et al. 2022).

X-ray Intensity
In Figure 7, we illustrate the contribution of BHs from various sources to the unresolved portions of the present-day hard ([2 − 10 keV]) and soft ([0.5−2 keV]) X-ray backgrounds, again considering PBHs with  c = 10 M ⊙ for a monochromatic mass function and  PBH = 10 −3 .The X-ray background intensity in a given energy band , where ℎ is the Planck constant.We emphasize that any unresolved sources are preferentially expected at  ≳ 6, given the low luminosity of typical high-redshift systems.It is thus natural to consider accreting BHs at cosmic dawn when accounting for the unresolved CXB.Note that when assessing the contribution of sources at  ′ to the radiation observed at  = 0, we use  ′ = (1 +  ′ ) for the integration involved in Equ.(48).
Our results for PBHs agree with the findings of Ziparo et al. (2022), when using identical parameter choices.However, by incorporating the effects of DM-baryon streaming and Compton drag, which diminish the accretion rate, we determine that PBHs in the IGM contribute only a minuscule fraction to the overall X-ray radiation.For  ≲ 30, the BH feedback power from haloes becomes dominant, similar to the behaviour in the (rest-frame) UV band.Comparing our predicted radiation levels to the unresolved CXB inferred from deep Chandra observations (Cappelluti et al. 2017(Cappelluti et al. , 2022)), we conclude that the accretion feedback from PBHs, located within haloes, may account for a substantial part of the missing X-ray flux.More specifically, in the soft X-ray band ([0.5 − 2 keV]), radiation from PBHs at  ≳ 6 almost reaches the intensity level of the unresolved X-ray background, whereas in the hard X-ray band ([2 − 10 keV]), radiation from these PBHs approaches the unresolved background intensity, but falls somewhat short.In both bands, the radiation from PBHs dominate over all other BH sources by at least an order of magnitude.
We note that neutral H and He in the IGM absorb soft X-ray photons ([0.5 − 2 keV]) more efficiently than hard X-ray radiation ([2 − 10 keV]), due to the approximately  −3 dependence of the photo-ionization cross sections (e.g.Wilms et al. 2000).As a result, the upper panel in Figure 7 represents the upper limit for the soft X-ray contribution from PBHs.In summary, our calculations suggest that radiation from accreting PBHs at high redshifts may play a significant role in accounting for the unresolved portion of the cosmic X-ray background, given a PBH abundance parameter of  PBH ≳ 10 −3 .

SUMMARY AND CONCLUSIONS
We have investigated the impact of PBHs as a dark matter component on structure formation and radiation backgrounds in the early Universe.Taking into account the presence of PBHs, we have analyzed the energy feedback resulting from the accretion of matter by all BHs of various origins across cosmic history.Our study has led us to the following key findings and conclusions: (i) From the constructed PBH-ΛCDM Universe, we find that even in the absence of initial clustering of PBHs, their presence will accelerate the formation of virialized structures (haloes), as the Universe becomes clumpy with the emergence of minihaloes of mass ≲ 10 6 M ⊙ .Increasing both PBH mass abundance and characteristic mass parameters will accelerate structure formation, with the former being more significant, as shown in Figure 2. (ii) In the IGM, only PBHs are present.However, within DM haloes, we need to consider BHs formed through astrophysical channels in addition to PBHs.As the first stars begin to form, the remnants of these massive stars give rise to SRBHs.However, even assuming that PBHs constitute only a small fraction of the DM (  PBH ≲ 10 −3 ), we find that HMXBs and SMBHs start to generate comparable accretion power as PBHs within halo environments only at lower redshifts  ≲ 7 (see Figs. 4 and 5).(iii) Considering the evolving bolometric radiation energy density produced by BH accretion, we have observed that the power density from PBH accretion in DM haloes begins to dominate over that from accretion in the IGM at  ≲ 30.Initially, the radiative output from PBHs dominates, but it becomes comparable to that from SRBHs around  ∼ 7.However, if the PBH abundance were much greater than  PBH ∼ 10 −3 /  esc , as shown in Figs. 5 and 6, the Universe would undergo reionization too early, in violation of existing constraints.This further supports the conclusion that PBHs within the solar mass range cannot be the primary component of dark matter.(iv) We have found that by changing the mass function of PBHs from a monochromatic distribution to a log-normal one, there is a minor change in the total energy density, by a factor of approximately 2.Moreover, with a larger mass dispersion parameter , the overall energy density further increases.This can be explained by the fact that the scaling of the accretion power  with respect to BH mass is non-linear with a power-law index > 2, as shown in Equ.(28-29), ( 37) and ( 44).However, since the combination of characteristic mass and PBH abundance parameters is highly constrained by Kühnel & Freese (2017), any effect from changing the mass function from monochromatic to log-normal is likely small.(v) By modeling the spectral energy distribution (SED) of BH accretion feedback (see Appendix B), we can calculate the energy feedback at specific frequencies (bands).From Figure 6, we conclude that accretion onto PBHs could be an important source for generating sufficient UV and Lyman- background radiation to affect cosmic reionization or to produce detectable absorption in the (redshifted) 21-cm line via efficient Wouthuysen-Field coupling.However, for  PBH ≲ 10 −3 , no existing constraints would be violated.On the other hand, when comparing the X-ray power from accreting BHs with the limits on the unresolved X-ray background, as shown in Figure 7, we find that the total power from BH sources in DM haloes is comparable to the present-day unresolved X-ray background limit.Again, the accretion power from PBHs in the IGM is negligible.For PBHs within haloes, the X-ray power generated is always dominant over that by other BHs, in both soft and hard X-ray bands, with the difference in the latter being even stronger.The PBH contribution may thus reach the intensity level of the unresolved soft X-ray background, but falls somewhat short of the level of the hard X-ray background, for plausible values of  PBH .We note that any conclusions for the soft X-ray background depend on the uncertain modeling of the optical depth to photo-ionization in the neutral IGM.
We would like to acknowledge several limitations in our study.Firstly, concerning the mass function of PBHs, we note that the predicted accretion rates in the IGM and virialized halos are negligibly small (see Equ. 37), rendering the mass function effectively time independent, within the characteristic mass range of  c ∼ 1 − 100 M ⊙ .Previous studies (e.g., Hasinger 2020;De Luca et al. 2020a;Yuan et al. 2023) have attempted to address this issue using a semi-analytical approach with accretion models from Mack et al. (2007); Ricotti et al. (2008), where the accretion rate could be significantly enhanced within haloes seeded by PBHs in the IGM environment.We note that in the case of extended or monochromatic PBH mass functions with a large characteristic mass  c ≫ 10 M ⊙ , the time-independent assumption will likely no longer be a good approximation.For PBHs in a virialized halo environment, changes to the mass function become more complex.Here, additional factors need to be considered, such as the distribution within haloes of PBHs with an extended mass function, the case of infalling PBHs through halo accretion from the cosmic web or through mergers with other haloes, and the merger of PBHs with other BHs.Selected recent studies have considered the possibility that BH binaries or massive BH seeds could form though PBH captures and mergers (see e.g.Hayasaki et al. 2016;Ali-Haïmoud et al. 2017;De Luca et al. 2020b, 2023).In this work, for the sake of simplicity, we ignored the effect of gravitational capture or merger, and assumed a time-independent PBH mass function, for both IGM and halo environments throughout cos-mic history.Therefore, a more accurate modeling of the change in the mass function can be achieved through cosmological simulations, which we leave as a topic for future studies.Any predictions for the radiation feedback from all BHs are subject to the broader uncertainties in modeling BH number densities and halo gas density profile (see footnote 9).As an example, not all galaxies will host a central SMBH, giving rise to AGN feedback.Furthermore, for HMXBs we may have overestimated the accretion rate and its duty cycle.Additionally, we have not considered the scattering/absorption attenuation caused by hydrogen and helium atoms due to the challenges in calculating source distributions.However, taking this additional factor into account would not significantly alter our main conclusions since it equally applies to BHs in haloes.Overall, we can still conclude that distinguishing PBHs from other BHs solely based on the cosmic radiation background is challenging.
In future work, we plan to test the robustness of our modeling, such as how halo density profiles may be modified in the presence of BHs, with cosmological simulations.Additionally, we will consider the PBH-generated emissions in other bands, in particular the nearinfrared and radio bands.Such studies will serve to provide a more complete interpretation for ongoing or upcoming sky surveys, such as with Euclid (Amendola et al. 2018) or the Square Kilometre Array (SKA) (Weltman et al. 2020), to detect the imprints left on the early IGM by the presence of PBH dark matter. .

(B1)
Here,   is the inner radius of the disk, which is equal to three times the Schwarzschild radius   .The temperature   at the inner radius is derived from the balance between gravitational energy release and radiative cooling: where  SB is the Stefan-Boltzmann constant, and  B the Boltzmann constant.
The maximum temperature  max =  ( max ) of the disk is attained at  max = 1.36  , as given by the temperature profile.The minimum temperature  o =  ( out ) is reached at the outer radius of the disk: (B3) This outer radius is determined by equating the angular momentum of the gas from the surrounding environment to the angular momentum of the accretion disk.Here, eff is the same effective velocity as in Equ. ( 27).
The overall shape of the spectrum is the result of the blackbody spectrum contributing from all parts of the accretion disk, taking the following form:  5, we show here the integrated (bolometric) energy density from BH accretion feedback of different origins vs. redshift: PBHs in the IGM (blue), in haloes (green), SRBHs (red), HMXBs (cyan), and SMBHs (yellow).For PBHs, we assume a monochromatic mass distribution with  c = 10 M ⊙ , and mass fractions of  PBH = 10 −3 (upper panel) and 10 −5 (lower panel).Here, we consider the same range for the abundance of BHs as shown in Figure 4, when deriving possible radiation outputs, but include the effect of halo enhancement on PBH accretion.

Figure 1 .
Figure1.Overall role of PBHs in early cosmic evolution.Initially, before the onset of structure formation, PBHs reside in the quasi-homogeneous intergalactic medium (IGM), as shown in the top left panel.As the Universe expands and cools, haloes of small masses (∼ 10 5 − 10 6 M ⊙ ) start to emerge from the gravitational instability, and PBHs cluster in those minihaloes, as illustrated in the top right panel.As haloes keep growing, the first stars form in sufficiently massive haloes from the collapse of primordial gas clouds.Stars with mass below the pair-instability limit (∼ 10 − 130 M ⊙ ) quickly die and collapse into stellar-remnant black holes (SRBHs).In certain cases, some of the stars form in binaries, where the more massive component becomes a black hole.Through accretion from the companion star, the black hole emits strong X-rays, giving rise to low/high mass X-ray binaries (L/HMXBs), depending on the mass of the BH (see the lower right panel).As some of the haloes have grown increasingly massive, in excess of ≳ 10 8 M ⊙ , the black hole in the center of the galactic bulge starts to accrete at nearly Eddington rate to grow into supermassive black holes (SMBHs), which exert vigorous radiative feedback during this process.This more complex stage in cosmic evolution, with the full panoply of BH sources, is summarized in the lower left panel.
) in turn determining the normalizing constant.Here, we will consider both top-heavy ( = −1.35;Stacy & Bromm 2013) and Salpeter ( = −2.35;Salpeter 1955) IMF slopes for Pop III stars, but only the Salpeter value for Pop I/II stars.Given the stellar mass function, a rough estimate for the BH mass function is:

Figure 3 .
Figure 3. Cumulative BH mass from star formation throughout cosmic history.The overall star formation rate was generated by the parametric fit from the UNIVERSE MACHINE (Zhang et al. 2023).We explicitly add Pop III star formation as the dominant contribution at high redshifts ( ≳ 15) (Liu & Bromm 2020a).The grey area represents the effective absence of star forming haloes, with negligible halo number densities at those redshifts.The white lines reproduce the halo mass growth histories of select haloes.

Figure 4 .
Figure 4.The evolution of BH density from various sources: all PBHs (blue), PBHs in haloes (green), SRBHs (red), HMXBs (blue), and SMBHs (yellow).Here, we show results for a monochromatic PBH mass distribution with  c = 10 M ⊙ , and a mass fraction in DM of  PBH = 10 −3 .For SRBHs, we assume instantaneous formation from the remnants of massive stars.To comprehensively calculate the BH mass budget for SRBH and HMXBs, we employ the star formation history of both Pop III and Pop I/II stars (Madau & Dickinson 2014; Liu & Bromm 2020a; Harikane et al. 2022, 2023a; Zhang et al. 2023), assuming a Salpeter IMF for Pop I/II stars, and a top-heavy IMF for Pop III stars.We show the maximum possible HMXB contribution (cyan dotted line), representing the limiting case of full duty cycle and accretion time of ∼ 10 Myr.For the other cases, we assign a duty cycle of 0.1 and an accretion time of ∼ 1 Myr.For the case of SMBHs, we use the models from Jeon et al. (2022) (dashed) and Zhang et al. (2023) (solid).Here, Mass densities are given in comoving units.

Figure 5 .
Figure 5. Integrated (bolometric) energy density from BH accretion feedback of different origins vs. redshift: PBHs in the IGM (blue), in haloes (green), SRBHs (red), HMXBs (cyan), and SMBHs (yellow).For PBHs, we assume a monochromatic mass distribution with  c = 10 M ⊙ , and mass fractions of  PBH = 10 −3 (upper panel) and 10 −5 (lower panel).Here, we consider the same range for the abundance of BHs as shown in Figure 4, when deriving possible radiation outputs.For comparison, we also plot the minimum Lyman- flux required for Wouthuysen-Field coupling (dot-dashed line; from Ciardi & Madau 2003), as well as the minimum UV flux required for reionization (dashed line and grey-shaded area; from Madau et al. 1999).
Figure6. Background intensity at the Lyman limit,  UV = 13.6 eV/ℎ, vs. redshift.Here, we again consider various BH sources (PBHs, SRBHs, HMXBs, and SMBHs), differentiating between halo and IGM PBH contributions.For the sake of simplicity, we assume a monochromatic mass function for PBHs with a characteristic mass of 10 M ⊙ and a mass fraction of 10 −3 .We show the minimum intensity,  UV , required for reionization(Madau et al. 1999), as well as the minimum Lyman- background,   , needed for efficient Wouthuysen-Field (WF) coupling(Ciardi & Madau 2003).For the comparison to the latter, we for simplicity assume   ≃  UV , and for both limits, we indicate uncertainties with a factor-of-ten range (grey shaded areas).We also indicate the timing constraint for effective WF coupling at  ∼ 15 − 20 (purple shaded area), inferred from the EDGES 21cm-cosmology measurement(Mirocha & Furlanetto 2019).Additionally, the brown region denotes the reionization limit at  ∼ 6.

Figure 7 .
Figure7.Contribution of accreting high- BH sources to the present-day ( = 0) X-ray background.X-ray intensities are plotted against redshift for all BH sources considered here, applying the same set of parameters as used in Figure6.The black dash-dotted line represents the unresolved Xray background(Cappelluti et al. 2017(Cappelluti et al. , 2022)).The upper panel shows the soft X-ray band with a photon energy in the [0.5 − 2keV] range, while the lower panel depicts the integrated hard X-ray band with photon energies of [2 − 10keV].

Figure A2 .
Figure A2.Average SMBH mass as a function of redshift and halo mass, from the parametric fit generated by the UNIVERSEMACHINE(Zhang et al. 2023).Here, we truncate the mass that could serve as SMBH host at 10 8 M ⊙ .Similar to Figure3, the grey area represents the combination of parameters where no halo formation is allowed, and white lines represent individual halo growth histories.

Figure C1 .
Figure C1.Effect of enhanced accretion in PBH-seeded haloes.Similar to Figure5, we show here the integrated (bolometric) energy density from BH accretion feedback of different origins vs. redshift: PBHs in the IGM (blue), in haloes (green), SRBHs (red), HMXBs (cyan), and SMBHs (yellow).For PBHs, we assume a monochromatic mass distribution with  c = 10 M ⊙ , and mass fractions of  PBH = 10 −3 (upper panel) and 10 −5 (lower panel).Here, we consider the same range for the abundance of BHs as shown in Figure4, when deriving possible radiation outputs, but include the effect of halo enhancement on PBH accretion.
Liu et al. (2021)13))), where one of the stars in the binary has an initial mass  i ≳  ★,c .To estimate the abundance of HMXBs, we consider the cases of Pop III and Pop I/II stars separately.For Pop III stars, we adopt the prescription from previous studies as describedin Jeon et al. (2014)andStacy & Bromm (2013);Liu et al. (2021), assuming that about ∼ 30% of Pop III stars form in binaries.On the other hand, observations suggest a larger binary fraction of approximately 0.5 in the local universe (e.g. Δ • , within a halo of mass  h at redshift .Furthermore,  ( • ) is the luminosity associated with BHs of mass  • , and ( • ) is the BH mass distribution function in a general form.To obtain the total number of BHs, we divide the total BH mass in a halo,  •,tot( h), by the average BH mass, as calculated in the denominator of the second equality above.When integrating,  •,min and  •,max represent the respective BH mass limits.For PBHs, we substitute the denominator with the characteristic BH mass  c , and integrate over the mass range of 0.1 − 10 c .For SRBHs and HMXBs, we employ Equ.(17- 19)with the mass limits discussed in Section 2.3.For SMBHs, we assume that each massive halo contains one central SMBH, and that the luminosity is given by the Eddington value, • ≃  Edd ( • ) ≃  0  Edd  2 .Once we have the aggregate power from one type of BH in a halo of mass  h at redshift , the (proper) average power density,  • (), is obtained by integrating over the halo mass function, as mentioned previously: •,i ( • ,  h ) is the number of BHs in mass bin  of width11For a typical value,  max ≲ 10 13 M ⊙ at redshift 6.
. As our results indicate, accreting PBHs located within haloes are possibly