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

Saiyang Zhang, 1,2 Boyuan Liu 3 and Volker Bromm2,4
1Department of Physics, University of Texas, Austin, TX 78712, USA
2Weinberg Institute for Theoretical Physics, University of Texas, Austin, TX
78712, USA
3Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge,
CB3 0HA, UK
4Department of Astronomy, University of Texas, Austin, TX 78712, USA E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

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.

###### keywords:

black hole physics – dark ages, reionization, first stars – dark matter –
diffuse radiation – early Universe

††pubyear: 2023††pagerange: Distinguishing the impact and signature of black
holes from different origins in early cosmic history–B

##  1 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, 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 $\sim$150
$\mathrm{M}_{\odot}$ (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 $M_{\mathrm{PBH}}\sim
1-100\ \mathrm{M}_{\odot}$, typical for stars and their remnants, provide all
of the cold dark matter (CDM) measured in the Universe, based on multiple
empirical constraints (e.g. Carr & Kühnel, 2020). Their contribution can be
expressed as $f_{\mathrm{PBH}}=\Omega_{\mathrm{PBH}}/\Omega_{\mathrm{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 $f_{\mathrm{PBH}}\lesssim 0.1-10^{-3}$. Furthermore, the non-
detection of the stochastic gravitational wave background (SGWB) by LVK gives
a similar but weaker constraint of $f_{\mathrm{PBH}}\lesssim 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 ($f_{\rm PBH}\lesssim 0.01$) as well (e.g.
Tristram 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, $f_{\mathrm{PBH}}\lesssim 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 near-infrared
(NIR) electromagnetic (EM) wave band [$0.7-5\ \mathrm{\mu m}$], astronomers
start to unveil an increasing number of galaxies at $z\gtrsim 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
$\Lambda$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; Liu & Bromm 2022).
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 $M_{\rm
PBH}\sim 1-100\ \rm M_{\odot}$ 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
$z\sim 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 ($\gtrsim 30\mathrm{\,M}_{\odot}$) in virialized haloes (e.g. Husain et
al., 2021), as well as SMBHs in more massive haloes with
$M_{\mathrm{h}}\gtrsim 10^{8}\mathrm{\,M}_{\odot}$ (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 $\sim 10^{5}-10^{13}\mathrm{\,M}_{\odot}$, whereas Ziparo et al.
(2022) only considered minihaloes without star formation ($\lesssim
10^{6}\mathrm{\,M}_{\odot}$). Also, we will consider the case of PBHs with an
extended mass function given by different physical formation mechanisms,
whereas Ziparo et al. (2022) assumed a monochromatic DM mass, and Cappelluti
et al. (2022) a mass function based on observational constraints (from
Hasinger, 2020). 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 $\Lambda$CDM cosmology with parameters: $\Omega_{\mathrm{m}}=0.3153$,
$\Omega_{\mathrm{b}}=0.04930$, $n_{\mathrm{s}}=0.96$, and $h=0.6736$ (Planck
Collaboration et al., 2020). The density fraction of CDM is then given by:
$\Omega_{\mathrm{DM}}=\Omega_{\mathrm{m}}-\Omega_{\mathrm{b}}$, with a
possible PBH sub-component of $\Omega_{\rm PBH}=f_{\rm PBH}\Omega_{\rm DM}$.

##  2 Black Holes in the Early Universe

Figure 1: 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 ($\sim
10^{5}-10^{6}\mathrm{\,M}_{\odot}$) 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 ($\sim 10-130\mathrm{\,M}_{\odot}$) 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 $\gtrsim
10^{8}\mathrm{\,M}_{\odot}$, 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 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 $\Lambda$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.

###  2.1 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 $z$, consisting of adiabatic and isocurvature terms,
$\delta_{\mathrm{ad}}(z)$ and $\delta_{\text{iso }}(z)$, respectively.
Adiabatic perturbations uniformly affect all mass-energy components, such that
$\delta_{\mathrm{ad}}(z)$ 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

$\delta_{\mathrm{iso}}(\mathbf{x})\equiv\frac{\delta\rho_{\rm
PBH}}{\bar{\rho}_{\rm DM}}=\frac{f_{\rm
PBH}}{\bar{n}_{\mathrm{PBH}}}\sum_{i}\delta_{D}\left(\mathbf{x}-\mathbf{x}_{i}\right)-1.$
(1)

Here, $\delta_{D}$ is the Dirac delta function, $\mathbf{x}_{i}$ marks the
position of the $i$-th PBH, and $\bar{n}_{\mathrm{PBH}}$ 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
$M_{\rm c}$, as stated in, e.g., Carr (1975); Choptuik (1993). In general, we
thus write the average density as

$\bar{n}_{\mathrm{PBH}}=\frac{f_{\mathrm{PBH}}}{M_{\rm
c}}\Omega_{\mathrm{DM}}\frac{3H_{0}^{2}}{8\pi G}\mbox{\ .}$ (2)

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 the PDM, such that the combined evolution of the
density perturbations in both DM (PDM) and PBHs is given as follows:

$\displaystyle\delta_{\mathrm{DM}}(z)$
$\displaystyle=\delta_{\mathrm{ad}}(z)+\delta_{\mathrm{iso}}(z)=\delta_{\mathrm{ad}}(z)+\left[T_{\mathrm{iso}}(z)-1\right]f_{\mathrm{PBH}}\delta_{\text{iso
}}^{0},$ (3) $\displaystyle\delta_{\mathrm{PBH}}(z)$
$\displaystyle=\delta_{\mathrm{DM}}(z)+\delta_{\text{iso
}}^{0}=\delta_{\mathrm{ad}}(z)+T_{\text{iso
}}^{\mathrm{PBH}}(z)\delta_{\text{iso }}^{0}.$

The evolution of adiabatic perturbations with respect to $z$ is expressed via
$\delta_{\mathrm{ad}}(z)=T_{\mathrm{ad}}(z)\delta_{\rm ad}^{0}$. Here,
$\delta_{\rm ad}^{0}$ and $\delta_{\rm iso}^{0}$ correspond to the primordial
adiabatic and PBH isocurvature perturbations, whereas $T_{\rm ad}(z)$ and
$T_{\rm iso}(z)$ are the linear transfer functions for the two modes, both
normalized to 1 as $z\rightarrow\infty$. The perturbation in PDM induced by
PBHs is $\delta_{\rm iso}(z)=\left[T_{\rm iso}(z)-1\right]f_{\rm
PBH}\delta_{\rm iso}^{0}$, and we can write the transfer function for PBH-
induced isocurvature fluctuations as $T^{\rm PBH}_{\rm iso}(z)=1+\left[T_{\rm
iso}(z)-1\right]f_{\rm PBH}$. We note that the transfer function is in turn
related to the structure growth factor, $T(z)=D(z)$. 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):

$\displaystyle D(z)$
$\displaystyle\approx\left(1+\frac{3\gamma}{2\alpha_{-}}s\right)^{\alpha_{-}},\quad
s=\frac{1+z_{\mathrm{eq}}}{1+z}$ (4) $\displaystyle\gamma$
$\displaystyle=\frac{\Omega_{m}-\Omega_{b}}{\Omega_{m}},\quad\alpha_{-}=\frac{1}{4}(\sqrt{1+24\gamma}-1),$

where $z_{\mathrm{eq}}\sim 3400$ represents the redshift at matter-radiation
equality. Since we restrict our treatment to the matter-dominated era with
$z\ll z_{\mathrm{eq}}$, we have $s\gg 1\ (z\rightarrow 0)$. Therefore,
$D_{\mathrm{iso}}\approx 3s/2$ (assuming $\gamma\approx 1$ and
$\alpha_{-}\approx 1$ to be constant).

Given the two-point correlation function for the density perturbations from
PBH DM:

$\left\langle\delta_{\mathrm{PBH}}(\mathbf{x})\delta_{\mathrm{PBH}}(\mathbf{0})\right\rangle=\frac{f_{\rm
PBH}^{2}}{\bar{n}_{\mathrm{PBH}}}\delta_{D}(\mathbf{x})+\xi_{\mathrm{PBH}}(x),$
(5)

where $\xi_{\mathrm{PBH}}(x)$ represents the reduced correlation function as a
function of separation distance $x=|\mathbf{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 $P_{\rm iso}(k)$ from PBHs approximates to $f_{\rm PBH}^{2}/\bar{n}_{\rm
PBH}$. We can then express the total power spectrum111In our calculation, we
have used the package Colossus (Diemer, 2018) to generate the $\Lambda$CDM
power spectrum, $P_{\rm\Lambda CDM}(k)$, as calibrated by Planck Collaboration
et al. (2020). at $z=0$ as follows:

$\displaystyle P_{\mathrm{tot}}(k)$ $\displaystyle=P_{\rm\Lambda
CDM}(k)+P_{\rm iso}(k),$ (6) $\displaystyle P_{\rm iso}(k)$
$\displaystyle=\left[f_{\rm PBH}(D_{0}-1)\right]^{2}/\bar{n}_{\rm
PBH}+T^{2}_{\rm mix}(k)P_{\rm\Lambda CDM}(k).$

Initially, the adiabatic and iso-curvature modes are not correlated. However,
at a later stage ($z\sim 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
$T_{\mathrm{mix}}$. Considering the inter-PBH separation scale,
$k_{\mathrm{PBH}}=\left(2\pi^{2}\bar{n}_{\mathrm{PBH}}\right)^{1/3}$, we use
the heuristic formula calibrated to simulation results in Liu et al. (2022):

$T_{\mathrm{mix}}^{2}(k)=\begin{cases}f_{\mathrm{PBH}}(D_{0}-1)^{2}D_{\mathrm{ad},0}^{-1}\left(k/k_{\mathrm{PBH}}\right)^{3},&k\leq
3k_{\mathrm{PBH}}\\\ 0,&k>3k_{\mathrm{PBH}}\end{cases},$ (7)

where $D_{0}=D(z=0)$ is the growth factor for the isocurvature mode, and
$D_{\rm ad,0}=D_{\rm ad}(z=0)/D_{\rm ad}(z=z_{\rm eq})$ for the adiabatic
mode222We denote the growth factor of the adiabatic mode as $D_{\rm ad}(z)$,
and adopt the semi-analytical fit from equ.(4.74)-(4.76) in Mo et al. (2010).,
both evaluated at redshift 0.

From the total power spectrum, we calculate the density contrast as a function
of the window mass $M$:

$\sigma(M)=\left[\frac{1}{2\pi^{2}}\int
P_{\mathrm{tot}}(k)W_{\mathrm{TH}}\left(kr_{\rm
M}\right)k^{2}dk\right]^{1/2},$ (8)

where $W_{\mathrm{TH}}\left(kr_{\rm M}\right)$ is the top-hat Fourier
transform of the spherical window function, and $r_{\rm M}$ corresponds to the
comoving radius associated with the mass window $M$. We note that for non-
linear scales, where $\sigma(M)\gtrsim 1$, this mass corresponds to a
virialized halo structure. Therefore, in the following, we will assume
$M_{h}=M$ for the halo mass. Employing the Press-Schechter formalism for the
halo distribution (Press & Schechter, 1974), the number density of virialized
structures (halo mass function) between $M_{\rm h}$ and $M_{\rm h}+dM_{\rm h}$
is:

$\displaystyle n(M_{\rm h},z)\mathrm{d}M_{\rm h}$
$\displaystyle=\sqrt{\frac{2}{\pi}}\frac{\bar{\rho}}{M_{\rm
h}^{2}}\frac{\delta_{\mathrm{c}}(z)}{\sigma(M_{\rm h})}$ (9)
$\displaystyle\times\exp\left(-\frac{\delta_{\mathrm{c}}^{2}(z)}{2\sigma^{2}(M_{\rm
h})}\right)\left|\frac{\mathrm{d}\ln\sigma(M_{\rm h})}{\mathrm{d}\ln M_{\rm
h}}\right|\mathrm{d}M_{\rm h},$

where $\bar{\rho}$ is the average comoving matter density, and the critical
overdensity factor, $\delta_{\mathrm{c}}(z)$, as a function of redshift is
given by

$\delta_{c}(z)\simeq\frac{1.686}{D_{\rm ad}(z)}.$ (10)

Integrating the halo number density over the mass range of interest and
dividing by the average mass density, the collapse fraction of the cosmic
matter is given by

$f_{\rm col}(>M_{\rm
h}|z)=\frac{1}{2}\text{erfc}\left(\frac{\delta_{c}(z)}{\sqrt{2}\sigma(R)}\right),$
(11)

and is plotted in Figure 2.

(a) Constant mass fraction $f_{\rm PBH}=0.001$, but considering changes in
mass scale: $M_{\rm c}=1\rm\ M_{\odot}$ (red) and $100\rm\ M_{\odot}$ (cyan).

(b) Varying PBH abundance, $f_{\rm PBH}=0.1$ (yellow) and $0.001$ (magenta),
keeping the mass scale constant at $M_{\rm c}=100\rm\ M_{\odot}$.

Figure 2: Collapse fraction in haloes larger than $10^{5}$ (pentagon),
$10^{6}$ (triangle), and $10^{8}M_{\odot}$ (circle), vs. redshift $z$. Here,
we consider the changes in the formation of structures in the presence of PBHs
(red, cyan and yellow, compared to the standard $\Lambda$CDM case (blue).

As can be seen, in the presence of PBHs, more minihaloes ($M_{\rm h}\sim
10^{5}\rm\ M_{\odot}$) form than in the standard $\Lambda$CDM case. Also, the
abundance of minihaloes will depend on the characteristic PBH mass, $M_{\rm
c}$, but more strongly on the PBH abundance. This behaviour can be explained
by the shift of total power spectrum $P_{\rm tot}$ in the presence of PBHs by
the isocurvature term $f_{\rm PBH}^{2}/\bar{n}_{\rm PBH}$ 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 $\gtrsim 100\rm\
M_{\odot}$, minihaloes of mass $\lesssim 10^{6}\rm\ M_{\odot}$ 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.

###  2.2 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
$\psi(M_{\mathrm{PBH}})$ with respect to PBH mass $M_{\mathrm{PBH}}$ is
related to the differential number density as follows:

$\frac{dn_{\mathrm{PBH}}}{d\ln
M_{\mathrm{PBH}}}=f_{\mathrm{PBH}}\Omega_{\mathrm{DM}}\frac{3H_{0}^{2}}{8\pi
G}\psi(M_{\mathrm{PBH}}).$ (12)

Here, $\psi(M_{\mathrm{PBH}})$ is the mass probability distribution function,
normalized to 1 when integrating over $\ln M_{\mathrm{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:

$\psi\left(M_{\mathrm{PBH}}\right)=\delta(M_{\mathrm{PBH}}-M_{\rm c}),$ (13)

where $M_{\rm 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:

$\psi\left(M_{\mathrm{PBH}}\right)=\frac{1}{\sqrt{2\pi}\sigma
M_{\mathrm{PBH}}}\exp\left(-\frac{\ln^{2}\left(M_{\mathrm{PBH}}/M_{\rm
c}\right)}{2\sigma^{2}}\right),$ (14)

where $\sigma$ denotes the width of the PBH distribution. Constraints on the
parameter space $(f_{\rm PBH},\sigma,M_{\mathrm{c}})$ for the lognormal
distribution have been calculated in Kühnel & Freese (2017), specifically
concluding that $0<\sigma\leq 1$. We note that as $\sigma\rightarrow 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). For
this paper, we restrict our attention to the two simple cases above. A
comprehensive study with respect to different shapes of the initial PBH mass
function will require another study.

###  2.3 Astrophysical Black Holes

####  2.3.1 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 optimistic333Recent JWST observations have discovered a
population of active galactic nuclei (AGN) at $z\gtrsim 5$ with a much higher
abundance than previously expected, which implies a surprising ubiquity of
SMBHs in high-$z$ galaxies (e.g. Furtak et al., 2023; Greene et al., 2023;
Harikane et al., 2023a; 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. assumption that SMBHs are located at the
center of all massive haloes with $M_{\rm h}\gtrsim 10^{8}\rm\ M_{\odot}$.
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:

$M_{\rm SMBH}=\eta(M_{\rm h},z)M_{\rm h},$ (15)

where the BH formation efficiency, $\eta(M_{\rm h},z)$, is a function of halo
mass and redshift. In Jeon et al. (2022), it is calculated using a global
energy balance argument as:

$\eta\simeq 2\times
10^{-7}\left(\frac{1+z}{10}\right)\left(\frac{f}{0.1}\right)^{-1}\left(\frac{M_{\rm
h}}{10^{10}\mathrm{M}_{\odot}}\right)^{2/3},$ (16)

where $f$ is the fraction of the BH rest mass energy deposited into the host
halo gas as heat. However, this idealized formalism scales as
$\eta\propto(1+z)$, and it does not converge to $0$ at high redshifts
($z\gtrsim 15$), where SMBHs begin to form. Therefore, we also adopt the
fitted value of $\eta(M_{\rm h},z)$ provided by the UNIVERSE MACHINE, which is
based on observational data (Zhang et al., 2023). This is discussed in more
detail in Appendix A.

####  2.3.2 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 star-formation 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
$M_{\mathrm{i}}$ can be summarized as follows (Heger et al., 2003):

  * •

$M_{\mathrm{i}}<10\mathrm{\,M}_{\odot}\rightarrow$ white dwarf

  * •

$10\mathrm{\,M}_{\odot}<M_{\mathrm{i}}<M_{\star,\mathrm{c}}\rightarrow$
neutron star

  * •

$M_{\star,\mathrm{c}}<M_{\mathrm{i}}<140\mathrm{\,M}_{\odot}\rightarrow$ black
hole

  * •

$140\mathrm{\,M}_{\odot}<M_{\mathrm{i}}<260\mathrm{\,M}_{\odot}\rightarrow$
pair-instability supernova (PISN) $\rightarrow$ no remnant

  * •

$M_{\mathrm{i}}>260\mathrm{\,M}_{\odot}\rightarrow$ massive black hole

Here, $M_{\star,\mathrm{c}}$ is the critical mass, before mass loss, above
which a single star will eventually produce a BH, with a value of
approximately $25\rm\ M_{\odot}$. 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 $Z\gtrsim 10^{-4}\ \rm Z_{\odot}$, 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, to maximize the production of
BHs, 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:

$\displaystyle\frac{dN}{dM_{\mathrm{i}}}=CM_{\mathrm{i}}^{\alpha},$ (17)

where $C$ is a normalization constant. For simplicity, we assume that any BH
with a progenitor mass larger than $260\mathrm{\,M}_{\odot}$ 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 $M_{\rm i,min}=0.1\mathrm{\,M}_{\odot}$
and $M_{\rm i,max}=260\mathrm{\,M}_{\odot}$ represent the minimum and maximum
stellar mass for our idealized IMF, respectively. Integrating, we obtain the
total stellar mass:

$M_{\star}=\int_{M_{\rm i,\min}}^{M_{\rm
i,\max}}\left(\frac{dN}{dM_{\mathrm{i}}}\right)M_{\mathrm{i}}dM_{\mathrm{i}}$
(18)

in turn determining the normalizing constant. Here, we will consider both top-
heavy ($\alpha=-1.35$; Stacy & Bromm 2013) and Salpeter ($\alpha=-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:

$\left(\frac{dN}{dM_{\rm
SRBH}}\right)=\left(\frac{dN}{dM_{\mathrm{i}}}\right)\frac{dM_{\mathrm{i}}}{dM_{\rm
SRBH}}.$ (19)

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\mathrm{\,M}_{\odot}$ or smaller
than the critical mass $M_{\star,\mathrm{c}}$, this function will return 0.

Overall, the total stellar BH mass in a halo of mass $M_{\rm h}$ at a given
redshift can be estimated by integrating over the halo star formation history
(SFH), as follows:

$\displaystyle M_{\bullet,\rm tot}\left(M_{\mathrm{h}},z\right)$
$\displaystyle=\int_{z_{\mathrm{ini}}}^{z}\mathrm{SFH}(M_{\mathrm{h}},z^{\prime})f_{\text{loss
}}(t(z)-t(z^{\prime})|M_{\mathrm{i}}\geq M_{\star,\mathrm{c}})$ (20)
$\displaystyle\times
f_{\mathrm{rem}}(t(z)-t(z^{\prime}))\left|\frac{dt}{dz^{\prime}}\right|\mathrm{d}z^{\prime}.$

In the above expression, $f_{\text{loss }}(t|M_{\mathrm{i}}\geq
M_{\star,\mathrm{c}})$ denotes the mass fraction of the stellar population
contained in stars with initial masses exceeding the critical mass
$M_{\star,\mathrm{c}}$ that have lifetimes shorter than $t$. This fraction is
a function of the time difference between the current redshift $z$ and the
redshift $z^{\prime}$ at which the stars initially form. For definiteness, we
assume that star formation begins at $z_{\rm ini}=50$.

We note that the lifetime of the massive stars ($\gtrsim 10\rm\ M_{\odot}$)
capable of forming BHs is less than $\sim\mathcal{O}(10\ \mathrm{Myr})$. This
is short in comparison to the time gap between the initial star formation and
the end of reionization, which is of order $\sim 900\mathrm{\,Myr}$. Thus, we
can remove the time dependence and approximate it with a nearly constant term
by integrating the IMF over the range of initial mass capable of forming
BHs444As noted earlier, the upper mass limit for this integration depends on
the stellar population. The maximum BH progenitor mass, $M_{\rm SRBH,\max}$,
is 140 $\rm\ M_{\odot}$ for Pop III stars and 260 $\rm\ M_{\odot}$ for Pop
I/II stars.:

$f_{\text{loss }}\simeq\frac{\int_{M_{\star,\mathrm{c}}}^{M_{\rm
i,\max}}\left(\frac{dN}{dM_{\mathrm{i}}}\right)M_{\mathrm{i}}dM_{\mathrm{i}}}{M_{\star}}$
(21)

In Equation (20), $f_{\mathrm{rem}}$ signifies the fraction of the progenitor
star’s mass converted into (astrophysical) BH remnants, given by the ratio of
the final BH mass, $M_{\mathrm{SRBH}}$, to the initial stellar mass,
$M_{\mathrm{i}}$. Upon integrating over the relevant mass interval, we arrive
at the remnant mass fraction, as follows:

$f_{\text{rem }}\simeq\frac{\int_{M_{\star,\mathrm{c}}}^{M_{\rm
i,\max}}\left(\frac{dN}{dM_{\mathrm{i}}}\right)M_{\rm
i}(M_{\mathrm{i}})dM_{\mathrm{i}}}{\int_{M_{\star,\mathrm{c}}}^{M_{\rm
i,\max}}\left(\frac{dN}{dM_{\mathrm{i}}}\right)M_{\mathrm{i}}dM_{\mathrm{i}}}.$
(22)

To model the star-formation history in a halo at redshift $z$, we separately
consider the formation history of Pop III and Pop I/II stars as a function of
redshift, represented as $\mathrm{SFH}(M_{\mathrm{h}},z)=\mathrm{SFH}_{\rm
III}(M_{\mathrm{h}},z)+\mathrm{SFH}_{\rm I/II}(M_{\mathrm{h}},z)$. In a first
order approximation, we can express this as the product of the comoving volume
corresponding to a given halo, and the cosmic star formation rate density
(SFRD):

$\mathrm{SFH}_{\rm i}(M_{\mathrm{h}},z)=\mathrm{SFRD}_{\rm i}(z)\times
V(M_{\rm h}).$ (23)

Here, $V(M_{\rm h})=8\pi GM_{h}/3H_{0}^{2}$ represents the comoving volume
corresponding to the halo mass $M_{\rm h}$, and the subscript $i$ 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 $z\lesssim 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
$\mathrm{SFH}_{\rm I/II}(M_{\mathrm{h}},z)$. The parametric fit, along with
the SMBH masses, are summarized in Appendix A. However, the formation rate at
$z\gtrsim 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 $z\sim 20$ of Pop III star formation:

$\frac{\mathrm{SFRD}_{\rm
III}(z)}{\mathrm{M}_{\odot}\mathrm{yr}^{-1}\mathrm{Mpc}^{-3}}=1.1\times
10^{-5}\frac{[(1+z)/21]^{2.63}}{0.054+[(1+z)/21]^{8.55}}.$ (24)

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 $\sim 10^{-4}M_{\rm h}$ at $z\simeq 6$.

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 ($z\gtrsim 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.

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
$z=9-16$ (Harikane et al., 2023b), with the work by Madau & Dickinson (2014)
serving as the upper bound. Therefore, at $z\gtrsim 10$, we can express the
upper bound and extrapolate it to even higher redshifts as follows:

$\frac{\mathrm{SFRD}_{\rm
I/II}(z)}{\mathrm{M}_{\odot}\mathrm{yr}^{-1}\mathrm{Mpc}^{-3}}\lesssim
0.007\frac{[(1+z)/10]^{2.7}}{0.001+[(1+z)/10]^{5.6}}$ (25)

####  2.3.3 X-ray Binaries

One important source of radiation feedback arises from the remnants of massive
stars that form in binaries. Specifically, we focus on high-mass X-ray
binaries (HMXBs, reviewed by e.g. Fornasini et al., 2023), where one of the
stars in the binary has an initial mass $M_{\mathrm{i}}\gtrsim
M_{\star,\mathrm{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 described in Jeon et al. (2014) and
Stacy & Bromm (2013); Liu et al. (2021), assuming that about $\sim 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. 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). For simplicity, we assume
that all binaries 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 $t_{\rm acc}$$\sim$ $1-10\rm\,Myr$ to all
HMXBs to allow them to accrete a substantial amount of mass
($\sim\mathcal{O}(\rm M_{\odot})$) from their companions (Mirabel et al.,
2011). We estimate the total mass of active HMXBs in a given halo of mass
$M_{h}$ as follows:

$\displaystyle M_{\rm HMXB}$
$\displaystyle\left(M_{\mathrm{h}},t(z)\right)/f_{\rm HMXB}$ (26)
$\displaystyle=$ $\displaystyle\int_{t(z)-t_{\rm
acc}}^{t(z)}\mathrm{SFH}(M_{\mathrm{h}},t^{\prime})f_{\text{loss
}}f_{\mathrm{rem}}\mathrm{d}t^{\prime}$ $\displaystyle\approx$ $\displaystyle\
\mathrm{SFH}(M_{\mathrm{h}},t(z))f_{\text{loss }}f_{\mathrm{rem}}t_{\rm acc}\
.$

Here, $f_{\rm HMXB}$ refers to the aforementioned fraction of massive stars
that form in binaries. The second equality is valid since $t_{\rm acc}$ is
much smaller than the cosmic time $t$. 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; Fragos et al., 2013b;
Sartorio et al., 2023; Liu et al., 2023), beyond the scope of the current
paper.

With the prescriptions above, 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 $M_{\rm c}=10\mathrm{\,M}_{\odot}$,
and a mass fraction in DM of $f_{\rm PBH}=10^{-3}$. We conclude that PBHs will
dominate the BH mass budget throughout the redshift range considered here,
unless $f_{\rm PBH}\lesssim 10^{-4}$. However, as more stars form and die,
around $z\sim 6$, the mass density in SRBHs (including HMXBs) becomes non-
negligible in dark matter haloes. Moreover, we notice that the SRBH fraction
in the form of active X-ray binaries is comparable to the abundance of SMBHs,
generally consistent with the result from Jeon et al. (2022).

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 $M_{\rm
c}=10\mathrm{\,M}_{\odot}$, and a mass fraction in DM of $f_{\rm
PBH}=10^{-3}$. For SRBHs, we assume instantaneous formation from the remnants
of massive stars. To comprehensively calculate the BH mass budget, 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, 2023b), 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 $\sim 10\rm\,Myr$.
For the other cases, we assign a duty cycle of 0.1 and an accretion time of
$\sim 1\rm\,Myr$. Mass densities are given in comoving units.

##  3 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):

$\dot{M}_{\bullet}=\frac{4\pi(GM_{\bullet})^{2}\rho_{\mathrm{gas}}}{v_{\rm
eff}^{3}}=\frac{4\pi(GM_{\bullet})^{2}\rho_{\mathrm{gas}}}{(c_{s}^{2}+v_{\mathrm{rel}}^{2})^{3/2}}.$
(27)

Here, $\rho_{\mathrm{gas}}=\mu nm_{\mathrm{p}}$ represents the average gas
density of the environment, given by the product of the average molecular
weight $\mu=1.22$, the gas number density $n$, and the proton mass
$m_{\mathrm{p}}$. Furthermore, $c_{\mathrm{s}}$ represents the sound speed
determined by the gas temperature, and $v_{\mathrm{rel}}$ is the relative
velocity of the BH with respect to the gas. Thus, the effective velocity is
given by $v_{\mathrm{eff}}=\sqrt{c_{\mathrm{s}}^{2}+v_{\mathrm{rel}}^{2}}$.

Next, we introduce a dimensionless accretion factor:

$\quad\dot{m}\equiv\frac{\dot{M}_{\bullet}}{\dot{M}_{\mathrm{Edd}}},$ (28)

where $\dot{M}_{\mathrm{Edd}}$ represents the Eddington accretion rate given
by:

$\dot{M}_{\mathrm{Edd}}=2.7\times
10^{-7}\mathrm{M}_{\odot}\mathrm{yr}^{-1}\left(\frac{M_{\bullet}}{10\mathrm{\,M}_{\odot}}\right)\left(\frac{\epsilon_{0}}{0.1}\right)^{-1}.$
(29)

Here, $\epsilon_{0}$ refers to the radiative efficiency. We assume that BH
accretion is Eddington-limited, such that $\dot{m}\leq 1$.

###  3.1 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). The accretion factor, normalized to
typical values, can be expressed as follows:

$\dot{m}=3.40\times
10^{-6}\lambda\left(\frac{1+z}{100}\right)^{3}\left(\frac{M_{\bullet}}{10\mathrm{\,M}_{\odot}}\right)\left(\frac{v_{\rm
eff}}{10\mathrm{~{}km}\mathrm{~{}s}^{-1}}\right)^{-3}.$ (30)

Here, $\lambda$ represents a scaling constant that represents the damping
effect, numerically fitted as follows:

$\lambda=\exp\left(\frac{9/2}{3+\hat{\beta}^{0.75}}\right)\frac{-1+(1+\hat{\beta})^{1/2}}{\hat{\beta}}.$
(31)

The effective viscosity $\hat{\beta}$ is in turn determined by the Compton
drag and Hubble flow, given by:

$\displaystyle\hat{\beta}=$ $\displaystyle 1.4\times
10^{-5}\left(\frac{M_{\bullet}}{10\mathrm{\,M}_{\odot}}\right)\left(\frac{1+z}{100}\right)^{3/2}\left(\frac{v_{\mathrm{eff}}}{10\mathrm{~{}km}\mathrm{~{}s}^{-1}}\right)$
(32) $\displaystyle\times\left[1+0.178\left(\frac{x_{\rm
e}}{10^{-3}}\right)\left(\frac{1+z}{100}\right)^{5/2}\right].$

Here, $x_{\rm e}$ represents the ionization fraction as a function of
redshift. For simplicity, we approximate $x_{\rm e}$ as a step function, in
line with the CMB constraint from Planck Collaboration et al. (2020). For
$6\lesssim z\lesssim 100$, which is long after recombination, we take its
value to be a constant555For the case of PBH accretion in the IGM, previous
studies have calculated the ionization fraction, e.g. Ziparo et al. (2022). In
their study, although $x_{\rm e}$ changes by an order of $\sim 100$ with
respect to $f_{\rm PBH}\sim 10^{-4}-1$, this factor has a negligible effect on
$\hat{\beta}$., specifically $x_{\rm e}\sim 10^{-3}$. For $z\lesssim 6$, after
reionization, this factor is set to 1. In all cases considered, we have
$\hat{\beta}\ll 1$. Then, the scaling constant can be approximated as
$\lambda\simeq\exp(3/2)/2\simeq 2.24$.

For the gas effective velocity $v_{\mathrm{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:

$\left\langle v_{\mathrm{eff}}\right\rangle\sim
c_{s}\left(\frac{16}{\sqrt{2\pi}}\mathcal{M}^{3}\right)^{\frac{1}{6}}\theta(\mathcal{M}-1)+c_{s}\left(1+\mathcal{M}^{2}\right)^{\frac{1}{2}}\theta(1-\mathcal{M}).$
(33)

Here, $\mathcal{M}=\left\langle v_{\mathrm{rel}}\right\rangle/c_{\rm s}$
represents the Mach number, and $\left\langle v_{\mathrm{rel}}\right\rangle$
is the root-mean-square of the relative velocity between dark matter (DM) and
baryons, approximated as666 We use a slightly different prescription for
$\left\langle v_{\mathrm{rel}}\right\rangle$ vs. $z$, as compared to Ricotti
et al. (2008), by considering the relative motion between DM and baryons in
non-linear perturbation theory (see e.g. Tseliakhovich & Hirata, 2010; Stacy
et al., 2011; Dvorkin et al., 2014).:

$\left\langle v_{\rm rel}\right\rangle\approx\min\left[1,z/10^{3}\right]\times
30\mathrm{~{}km}\mathrm{\,s}^{-1}.$ (34)

The gas sound speed, $c_{\rm 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):

$c_{s}\simeq
1.80\left(\frac{1+z}{100}\right)^{1/2}\left[\left(\frac{1+z_{\mathrm{dec}}}{1+z}\right)^{\beta}+1\right]^{-1/2\beta}\mathrm{km}\mathrm{\,s}^{-1}.$
(35)

Here, $\beta=1.72$ is a fitting parameter, and $z_{\mathrm{dec}}\simeq 130$
the decoupling redshift for baryons from the radiation field. Combining Equ.
(29) and (30), we can approximate the mass growth rate as:

$\dot{M}_{\bullet}\sim 2.7\times 10^{-7}\dot{m}\
\left(\frac{M_{\bullet}(t)}{10\mathrm{\,M}_{\odot}}\right)\left(\frac{\epsilon_{0}}{0.1}\right)^{-1}\mathrm{M}_{\odot}\
\mathrm{yr}^{-1}.$ (36)

Using the above equations to calculate the accretion rate and integrating over
cosmic time, we find that for PBH masses of interest, $M_{\rm PBH}\lesssim
100\rm\ M_{\odot}$, the change in mass is $\delta M_{\rm PBH}/M_{\rm PBH}\ll
0.01$, which is negligible. However, we also note that, although yet to be
observed, the accretion rate of PBHs could be greatly enhanced if DM haloes
are seeded by PBHs in isolation (e.g. Ricotti, 2007; Mack et al., 2007;
Ricotti et al., 2008; De Luca et al., 2020a). In this senario, the PBH mass is
associated with its halo mass through $M_{\rm PBH}\propto M_{\rm h}/(1+z)$,
and can accrete a significant amount of mass by $z\lesssim 6$, if $10\rm\
M_{\odot}$ is the characteristic PBH mass scale. To carefully address this
uncertainty, another simulation study would be required. Therefore, to make
conservative estimates and simplify our calculations, for PBH accretion in the
IGM, we assume a mass function $\psi(M)$ that is independent of redshift.

###  3.2 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:

$\displaystyle\dot{m}=3.2\times
10^{-6}\left(\frac{M_{{}_{\bullet}}}{10\mathrm{\,M}_{\odot}}\right)\left(\frac{n(r)}{1\mathrm{\,cm}^{-3}}\right)\left(\frac{\mu}{1.22}\right)\left(\frac{v_{\rm
eff}}{10\mathrm{\,km}/\mathrm{s}}\right)^{-3}\mbox{.}$ (37)

Here, for a BH at a given radius $r$ relative to the center of the halo,
$n(r)$ 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: $n(r)\propto r^{-2}$. Additionally, we assume
the same scaling for the PBH number density. Meanwhile, we use the virial
velocity as the effective velocity $v_{\rm eff}$ between BHs and the gas:

$v_{\rm eff}\sim\sqrt{\frac{GM_{\mathrm{h}}}{R_{\mathrm{vir}}}}\sim
5.4\mathrm{~{}km}\mathrm{~{}s}^{-1}\left(\frac{M_{\mathrm{h}}}{10^{6}\mathrm{\,M}_{\odot}}\right)^{\frac{1}{3}}\left(\frac{1+z}{21}\right)^{-\frac{1}{2}}.$
(38)

Here, $R_{\rm vir}$ is the virial radius of the halo, which is related to halo
mass and redshift by:

$R_{\mathrm{vir}}\sim 0.147\left(\frac{\Delta_{\rm
c}}{18\pi^{2}}\right)^{-1/3}\left(\frac{M_{\rm
h}}{10^{6}\mathrm{M}_{\odot}}\right)^{1/3}\left(\frac{1+z}{21}\right)^{-1}\mathrm{kpc},$
(39)

where $\Delta_{\rm 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 $\dot{m}\lesssim 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:

$\displaystyle n(R_{\rm vir})$
$\displaystyle\sim\Omega_{\mathrm{b}}(1+z)^{3}\frac{3H_{0}^{2}}{8\pi
G}\frac{\Delta_{\rm c}}{\mu m_{p}}$ (40) $\displaystyle\sim
0.33\left(\frac{\Delta_{\rm
c}}{18\pi^{2}}\right)\left(\frac{\mu}{1.22}\right)^{-1}\left(\frac{1+z}{21}\right)^{3}\rm
cm^{-3}.$

Similarly, since no previous simulations have considered the co-evolution 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 $\gtrsim 260\rm\ M_{\odot}$), 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 $\gtrsim 260\rm\ M_{\odot}$ and restrict
their radii to be outside of $0.001R_{\rm vir}$. Typically, for a $10^{9}\
\mathrm{M}_{\odot}$ halo at $z\simeq 10$, the virial radius is $R_{\rm
vir}\sim 3\ \rm kpc$. The inner region, $r<0.001R_{\rm vir}$, thus has a size
$\sim\mathcal{O}(\rm pc)$, and can only contain a SMBH. For PBHs following a
log-normal mass function, we numerically integrate from $0.1-10M_{\rm 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 $f_{\rm PBH}$:

$\displaystyle\langle\dot{M}_{\rm PBH,tot}\rangle\sim$ $\displaystyle 0.02\rm\
M_{\odot}\mathrm{yr}^{-1}\left(\frac{M_{\rm h}}{10^{12}\rm\
M_{\odot}}\right)^{1.094}\left(\frac{f_{\rm PBH}}{10^{-3}}\right)$ (41)
$\displaystyle\times(1+1.75z)\sqrt{\Omega_{m}(1+z)^{3}+(1-\Omega_{m})}\mbox{\
.}$

This is of the same order of magnitude as the accretion rate obtained by
integrating over all PBHs with a characteristic mass of $100\rm\ M_{\odot}$ in
a $10^{12}M_{\rm h}$ halo at $z=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. $\sim 0.01\
\mathrm{SFH}(M_{\rm h},z)$, as calculated using Equ.(20), which is also
significantly larger than the total accretion rate. 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 impose
Eddington-limited accretion throughout ($\dot{m}=1$).

##  4 Radiative signature

###  4.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 $u$ from BH feedback into two parts based on the BH habitat:

$u=u_{\mathrm{IGM}}+u_{\mathrm{halo}},$ (42)

where $u_{\rm IGM}$ represents the energy feedback from PBHs accreting in the
IGM, and $u_{\rm 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:

$u_{\rm halo}=\begin{cases}u_{\rm PBH},&\\\ \hskip 86.11084pt(M_{\min}(z)\leq
M_{\rm h}<M_{\rm H_{2}-cool})&\\\ u_{\rm SRBH}+u_{\rm HMXB}+u_{\rm PBH},&\\\
\hskip 86.11084pt(M_{\rm H_{2}-cool}\leq M_{\rm
h}<10^{8}\mathrm{M}_{\odot})&\\\ u_{\rm SRBH}+u_{\rm HMXB}+u_{\rm PBH}+u_{\rm
SMBH}.&\\\ \hskip 86.11084pt(10^{8}\mathrm{M}_{\odot}\leq M_{\rm
h}<M_{\max}(z))&\end{cases}$ (43)

Here, several DM halo mass limits are introduced. Specifically, $M_{\min}(z)$
refers to the smallest halo mass that could host a PBH, given by
$M_{\min}=\max(M_{c}\Omega_{\rm m}/f_{\rm PBH}\Omega_{\rm DM},M_{\rm vir})$,
where $M_{\rm vir}$ refers to the smallest mass that can virialize at a given
redshift into a halo structure with a virial temperature $\gtrsim
100\mathrm{\,K}$. We also notice some other pertinent mass scales, such as
$\sim M_{\rm c}(1+z_{\rm eq})/(1+z)$, 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. We defer a more complete understanding of halo
accretion at high redshifts to future study with simulations.

$M_{\rm H_{2}-cool}$ is the smallest mass of a DM halo capable of forming Pop
III stars via H2 cooling, and it is fitted by $M_{{\rm
H}_{2}-\mathrm{cool}}\simeq 1.54\times 10^{5}\rm\
M_{\odot}\left[(1+z)/31\right]^{-2.074}$ (Trenti & Stiavelli, 2009). $M_{\rm
max}(z)$ is the maximum allowed mass of DM haloes at redshift $z$, given by
the progenitor of $\sim 10^{16}\rm\ M_{\odot}$ haloes at $z=0$, following
their growth history777For a typical value, $M_{\rm max}\lesssim 10^{13}\rm\
M_{\odot}$ at redshift 6.. Haloes with mass larger than $M_{\rm max}(z)$ 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}\rm\ M_{\odot}$ 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 $\dot{m}$ and bolometric luminosity is given by:

$L=\epsilon_{\mathrm{EM}}\dot{M}_{\bullet}c^{2},\
\epsilon_{\mathrm{EM}}=\frac{\epsilon_{0}A\dot{m}}{1+A\dot{m}}.$ (44)

where $\epsilon_{\mathrm{EM}}$ is the electromagnetic efficiency, $A=100$, and
$\epsilon_{0}=0.057$ the radiative efficiency for optically thin accretion.
Note that $\epsilon_{\mathrm{EM}}$ reduces to $\epsilon_{0}$ for higher
accretion rates, $\dot{m}\gg 0.01$.

To evaluate the total luminosity from BH accretion in a halo of mass $M_{\rm
h}$ at redshift $z$, 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:

$\displaystyle L_{\bullet}(M_{\rm h},z)$ $\displaystyle=\sum_{i}\Delta
N_{\bullet,\mathrm{i}}(M_{\bullet},M_{\rm h})L(M_{\bullet})$ (45)
$\displaystyle=M_{\bullet,\mathrm{tot}}(M_{\rm
h})\frac{\int_{M_{\bullet,\min}}^{M_{\bullet,\max}}\psi(M_{\bullet})L(M_{\bullet})dM_{\bullet}}{\int_{M_{\bullet,\min}}^{M_{\bullet,\max}}\psi(M_{\bullet})M_{\bullet}dM_{\bullet}},$

where $\Delta N_{\bullet,\mathrm{i}}(M_{\bullet},M_{\rm h})$ is the number of
BHs in mass bin $i$ of width $\Delta M_{\bullet}$, within a halo of mass
$M_{\rm h}$ at redshift $z$. Furthermore, $L(M_{\bullet})$ is the luminosity
associated with BHs of mass $M_{\bullet}$, and $\psi(M_{\bullet})$ 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, $M_{\bullet\mathrm{,tot}}(M_{\rm
h})$, by the average BH mass, as calculated in the denominator of the second
equality above. When integrating, $M_{\bullet,\min}$ and $M_{\bullet,\max}$
represent the respective BH mass limits. For PBHs, we substitute the
denominator with the characteristic BH mass $M_{\rm c}$, and integrate over
the mass range of $0.1-10M_{\rm 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, $L_{\bullet}\simeq L_{\rm
Edd}(M_{\bullet})\simeq\epsilon_{0}\dot{M}_{\rm Edd}c^{2}$.

Once we have the aggregate power from one type of BH in a halo of mass $M_{\rm
h}$ at redshift $z$, the (proper) average power density, $j_{\bullet}(z)$, is
obtained by integrating over the halo mass function, as mentioned previously:

$j_{\bullet}(z)=\int_{M_{\rm min}(z)}^{M_{\rm max}(z)}(1+z)^{3}n\left(M_{\rm
h},z\right)L(M_{\rm h},z)dM_{\rm h},$ (46)

where $n\left(M_{\rm h},z\right)$ is the comoving number density of haloes per
unit halo mass, given by Equ.(9). The factor of $(1+z)^{3}$ is applied to
derive the physical number density at a given $z$. We note that, if the
maximum halo mass $M_{\rm max}$ is smaller than the lower limit $M_{\rm min}$
at some redshift $z$, there is no contribution from that particular type of
BH. The lower integration limit $M_{\rm min}$ depends on the BH type. For PBHs
located in a halo, $M_{\rm min}=\max(M_{c}\Omega_{\rm m}/f_{\rm
PBH}\Omega_{\rm DM},M_{\rm vir})$. For SRBHs, we choose $M_{\rm min}=M_{\rm
H_{2}-cool}$, whereas for SMBHs, we adopt $M_{\rm
min}=10^{8}\mathrm{\,M}_{\odot}$ as the minimum halo mass that can host an
AGN. The proper bolometric power density $u$ at redshift $z$ from BHs can be
obtained by integrating the power density at a given redshift $z^{\prime}$
from different sources over cosmic time:

$u(z)=\int_{z_{\rm
i}}^{z}\left(\frac{1+z}{1+z^{\prime}}\right)^{4}j_{\bullet}(z^{\prime})\left|\frac{dt}{dz^{\prime}}\right|dz^{\prime}\mbox{\
,}$ (47)

where $z_{\rm i}$ and $z$ are the initial and final redshift of consideration,
respectively. Here, we consider the most distant source emitting at $z_{\rm
i}\simeq 100$.

The results are shown in Figure 5 for different PBH mass fractions: $f_{\rm
PBH}=10^{-3}$ (upper panel) and $f_{\rm PBH}=10^{-1}$ (lower panel), both
assuming a monochromatic mass distribution with $M_{\rm
c}=10\mathrm{\,M}_{\odot}$. 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 $z\lesssim 30$. We also note that, since SMBHs
and HMXBs are accreting at a much higher rate than SRBHs, their emission power
at lower redshift ($z\lesssim 10$) starts to surpass that from SRBHs.

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-$\alpha$ photons scatter with H atoms, thus indirectly
modifying the hyperfine level populations (Madau et al., 1997). We can
estimate the required Lyman-$\alpha$ energy density, in proper units and
assuming a flat spectrum, as follows: $u_{\rm Ly\alpha}\simeq 4\pi
J_{\alpha}\nu_{\alpha}/c\gtrsim 3.5\times 10^{-14}[(1+z)/7]\ \rm erg\,cm^{-3}$
(Ciardi & Madau, 2003). Similarly, we can estimate the minimum energy density
in H-ionizing UV radiation required to reionize the IGM: $u_{\rm UV}\simeq
4\pi J_{\rm UV}\nu_{\rm UV}/c\gtrsim 2\times 10^{-15}[(1+z)/7]^{1.5}\ \rm
erg\,cm^{-3}$(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.

(a) PBH mass fraction: $f_{\rm PBH}=10^{-3}$

(b) PBH mass fraction: $f_{\rm PBH}=10^{-1}$

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 $M_{\rm c}=10\rm\ M_{\odot}$, and mass
fractions of $f_{\rm PBH}=10^{-3}$ (upper panel) and $10^{-1}$ (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-$\alpha$ 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).

The features of the emerging spectrum depend mostly on the normalized
accretion rate $\dot{m}$, 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 obtain the (angle-averaged) radiation intensity
at frequency $\nu$ and a given redshift $z$, we need to calculate the specific
emission coefficient, $j_{\bullet,\nu^{\prime}}(z^{\prime})$, from sources at
higher redshifts $z^{\prime}>z$ with $\nu^{\prime}=\nu(1+z^{\prime})/(1+z)$.
Applying the same procedure as for Equ. (45-46), we first sum over the
contributions from all BHs in a halo, and then integrate over the halo mass
function. Finally, integrating over all contributions from preceding
redshifts, the resulting specific intensity at $z$ is (e.g. Schauer et al.,
2019):

$J_{\nu}(z)=\frac{c}{4\pi}\int_{z_{\rm
i}}^{z}\left(\frac{1+z}{1+z^{\prime}}\right)^{3}j_{\bullet,\nu^{\prime}}(z^{\prime})\left|\frac{dt}{dz^{\prime}}\right|dz^{\prime}.$
(48)

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,
$z_{\rm 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 $J_{\nu}$ over the corresponding band, and $z_{\rm i}$ denotes the
highest redshift below which radiation produced by BHs can reach $z$ 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 $f_{\rm 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).

###  4.2 UV Intensity

Figure 6: Background intensity at the Lyman limit, $\nu_{\rm UV}={\rm
13.6\,eV}/h$, 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\ \rm M_{\odot}$ and a mass
fraction of $10^{-3}$. We show the minimum intensity, $J_{\rm UV}$, required
for reionization (Madau et al., 1999), as well as the minimum Lyman-$\alpha$
background, $J_{\alpha}$, needed for efficient Wouthuysen-Field (WF) coupling
(Ciardi & Madau, 2003). For the comparison to the latter, we for simplicity
assume $\nu_{\alpha}\simeq\nu_{\rm 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 $z\sim 15-20$ (purple
shaded area), inferred from the EDGES 21cm-cosmology measurement (Mirocha &
Furlanetto, 2019). Additionally, the brown region denotes the reionization
limit at $z\sim 6$.

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 $M_{\rm c}=10\ \rm M_{\odot}$ and a mass fraction
$f_{\rm PBH}=10^{-3}$. We evaluate the radiation intensity at $\nu_{\rm
UV}=13.6~{}\mathrm{eV}/h$, and assume for simplicity the same value at
$\nu_{\alpha}$, i.e. $J_{\alpha}\simeq J_{\rm 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 $z\lesssim 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 $\sim 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
$z\sim 6$ (reviewed in Robertson et al. 2010), with early evidence based on
the Gunn-Peterson absorption trough in high-$z$ 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 $f_{\rm PBH}\gtrsim
10^{-3}$. For our fiducial value of cosmic PBH density (0.1% of all DM), we
can see that at $z\lesssim 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 $z\sim 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 ($f_{\rm 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 $f_{\rm PBH}\lesssim
10^{-3}/f_{\rm esc}$.

Similarly, we assess whether the PBH contribution to the cosmic Lyman-$\alpha$
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 $z\simeq 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 $f_{\rm PBH}\lesssim 10^{-3}$ for $M_{\rm
c}=10\rm\ M_{\odot}$. Using the same characteristic mass, we find a similar
constraint on the PBH abundance, as $f_{\rm 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-$\alpha$ 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). As our results indicate, accreting PBHs located within
haloes are possibly potent sources for Lyman-$\alpha$ background radiation at
high redshifts, whereas the other BH sources considered here fall short in
producing a significant number of Lyman-$\alpha$ photons. For a discussion of
the contribution of stellar sources to the cosmic Lyman-$\alpha$ 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).

###  4.3 X-ray Intensity

Figure 7: Contribution of accreting high-$z$ BH sources to the present-day
($z=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
Figure 6. To show the (very small) radiation contribution from PBHs in the
IGM, using the same scaling, we artificially multiply by a factor of $10^{6}$,
as denoted in the legend. The black dash-dotted line represents the unresolved
X-ray 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-2\rm keV]$
range, while the lower panel depicts the integrated hard X-ray band with
photon energies of $[2-10\rm keV]$.

In Figure 7, we illustrate the contribution of BHs from various sources to the
unresolved portions of the present-day hard ($[2-10~{}\mathrm{keV}]$) and soft
($[0.5-2~{}\mathrm{keV}]$) X-ray backgrounds, again considering PBHs with
$M_{\rm c}=10\ \rm M_{\odot}$ for a monochromatic mass function and $f_{\rm
PBH}=10^{-3}$. We emphasize that any unresolved sources are preferentially
expected at $z\gtrsim 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 $z^{\prime}$ to the radiation observed at $z=0$, we use
$\nu^{\prime}=(1+z^{\prime})\nu$ 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 $z\lesssim 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~{}\rm keV]$), radiation from PBHs at $z\gtrsim
6$ almost reaches the intensity level of the unresolved X-ray background,
whereas in the hard X-ray band ($[2-10~{}\rm 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~{}\mathrm{keV}]$) more efficiently than hard X-ray radiation
($[2-10~{}\mathrm{keV}]$), due to the approximately $\nu^{-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 $f_{\rm PBH}\gtrsim 10^{-3}$.

##  5 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-$\Lambda$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 $\lesssim 10^{6}\rm\ M_{\odot}$.
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 ($f_{\rm PBH}\lesssim 10^{-3}$), we find that HMXBs and SMBHs start to
generate comparable accretion power as PBHs within halo environments only at
lower redshifts $z\lesssim 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 $z\lesssim
30$. Initially, the radiative output from PBHs dominates, but it becomes
comparable to that from SRBHs around $z\sim 7$. However, if the PBH abundance
were much greater than $f_{\rm PBH}\sim 10^{-3}/f_{\rm 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 $\sigma$, the overall energy density further increases.
This can be explained by the fact that the scaling of the accretion power $L$
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-$\alpha$ background
radiation to affect cosmic reionization or to produce detectable absorption in
the (redshifted) 21-cm line via efficient Wouthuysen-Field coupling. However,
for $f_{\rm PBH}\lesssim 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 $f_{\rm 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 are negligibly small (see Equ. 37), rendering the mass
function effectively time independent, within the characteristic mass range of
$M_{\rm c}\sim 1-100~{}\mathrm{M}_{\odot}$. However, for PBHs in a halo
environment, changes to the mass function become more complex. 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). The accretion
rate could be significantly enhanced by haloes seeded by PBHs, and thus
changing the mass function for $M_{\rm c}\gg 10~{}\mathrm{M}_{\odot}$.
However, there are additional factors that need to be considered, such as the
distribution of PBHs with an extended mass function in haloes, the case of
infalling PBHs through halo accretion from the IGM or through mergers with
other haloes, and the merger of PBHs with other BHs. In this work, for the
sake of simplicity, we assume a time-independent PBH mass function, for both
IGM and halo environments throughout cosmic history. 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. 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 near-infrared 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.

## Acknowledgements

BL is supported by the Royal Society University Research Fellowship.

## Data Availability

The data and codes used in this article will be shared upon reasonable request
to the corresponding authors.

## References

  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2020) Abbott R., et al., 2020, Phys. Rev. Lett., 125, 101102
  * Abbott et al. (2021) Abbott R., et al., 2021, ApJ, 913, L7
  * Afshordi et al. (2003) Afshordi N., McDonald P., Spergel D. N., 2003, ApJ, 594, L71
  * Agazie et al. (2023) Agazie G., et al., 2023, ApJ, 951, L8
  * Ali-Haïmoud & Kamionkowski (2017) Ali-Haïmoud Y., Kamionkowski M., 2017, Phys. Rev. D, 95, 043534
  * Ali-Haïmoud et al. (2017) Ali-Haïmoud Y., Kovetz E. D., Kamionkowski M., 2017, Phys. Rev. D, 96, 123523
  * Amendola et al. (2018) Amendola L., et al., 2018, Living Reviews in Relativity, 21, 2
  * Bagui & Clesse (2022) Bagui E., Clesse S., 2022, Physics of the Dark Universe, 38, 101115
  * Becker et al. (2001) Becker R. H., et al., 2001, AJ, 122, 2850
  * Behroozi & Silk (2018) Behroozi P., Silk J., 2018, MNRAS, 477, 5382
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57
  * Behroozi et al. (2019) Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488, 3143
  * Behroozi et al. (2020) Behroozi P., et al., 2020, MNRAS, 499, 5702
  * Bertone et al. (2012) Bertone G., Cumberbatch D., Ruiz de Austri R., Trotta R., 2012, J. Cosmology Astropart. Phys., 2012, 004
  * Bird et al. (2016) Bird S., Cholis I., Muñoz J. B., Ali-Haïmoud Y., Kamionkowski M., Kovetz E. D., Raccanelli A., Riess A. G., 2016, Phys. Rev. Lett., 116, 201301
  * Bondi (1952) Bondi H., 1952, MNRAS, 112, 195
  * Bowman et al. (2018) Bowman J. D., Rogers A. E. E., Monsalve R. A., Mozdzen T. J., Mahesh N., 2018, Nature, 555, 67
  * Boylan-Kolchin (2023) Boylan-Kolchin M., 2023, Nature Astronomy, 7, 731
  * Cappelluti et al. (2017) Cappelluti N., et al., 2017, ApJ, 837, 19
  * Cappelluti et al. (2022) Cappelluti N., Hasinger G., Natarajan P., 2022, ApJ, 926, 205
  * Carr (1975) Carr B. J., 1975, ApJ, 201, 1
  * Carr & Kühnel (2020) Carr B., Kühnel F., 2020, Annual Review of Nuclear and Particle Science, 70, 355
  * Carr & Silk (2018) Carr B., Silk J., 2018, MNRAS, 478, 3756
  * Carr et al. (2021) Carr B., Kohri K., Sendouda Y., Yokoyama J., 2021, Reports on Progress in Physics, 84, 116902
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chisholm (2006) Chisholm J. R., 2006, Phys. Rev. D, 73, 083504
  * Choptuik (1993) Choptuik M. W., 1993, Phys. Rev. Lett., 70, 9
  * Ciardi & Madau (2003) Ciardi B., Madau P., 2003, ApJ, 596, 1
  * Clark et al. (2023) Clark J. S., Lohr M. E., Najarro F., Patrick L. R., Ritchie B. W., 2023, MNRAS,
  * Clesse & García-Bellido (2017) Clesse S., García-Bellido J., 2017, Physics of the Dark Universe, 15, 142
  * De Luca et al. (2020a) De Luca V., Franciolini G., Pani P., Riotto A., 2020a, Phys. Rev. D, 102, 043505
  * De Luca et al. (2020b) De Luca V., Franciolini G., Pani P., Riotto A., 2020b, J. Cosmology Astropart. Phys., 2020, 052
  * De Luca et al. (2020c) De Luca V., Franciolini G., Pani P., Riotto A., 2020c, J. Cosmology Astropart. Phys., 2020, 044
  * De Luca et al. (2020d) De Luca V., Desjacques V., Franciolini G., Riotto A., 2020d, J. Cosmology Astropart. Phys., 2020, 028
  * De Luca et al. (2022a) De Luca V., Franciolini G., Riotto A., Veermäe H., 2022a, Phys. Rev. Lett., 129, 191302
  * De Luca et al. (2022b) De Luca V., Franciolini G., Kehagias A., Pani P., Riotto A., 2022b, Physics Letters B, 832, 137265
  * Decarli et al. (2018) Decarli R., et al., 2018, ApJ, 854, 97
  * Desjacques & Riotto (2018) Desjacques V., Riotto A., 2018, Phys. Rev. D, 98, 123533
  * Diemer (2018) Diemer B., 2018, ApJS, 239, 35
  * Dizgah et al. (2019) Dizgah A. M., Franciolini G., Riotto A., 2019, Journal of Cosmology and Astroparticle Physics, 2019, 001
  * Dolgov & Silk (1993) Dolgov A., Silk J., 1993, Phys. Rev. D, 47, 4244
  * Dvorkin et al. (2014) Dvorkin C., Blum K., Kamionkowski M., 2014, Phys. Rev. D, 89, 023519
  * Escrivà et al. (2023) Escrivà A., Bagui E., Clesse S., 2023, J. Cosmology Astropart. Phys., 2023, 004
  * Feng (2010) Feng J. L., 2010, ARA&A, 48, 495
  * Fialkov & Barkana (2014) Fialkov A., Barkana R., 2014, MNRAS, 445, 213
  * Field (1959) Field G. B., 1959, ApJ, 129, 536
  * Finkelstein et al. (2022) Finkelstein S. L., et al., 2022, ApJ, 940, L55
  * Fornasini et al. (2023) Fornasini F. M., Antoniou V., Dubus G., 2023, arXiv e-prints, p. arXiv:2308.02645
  * Fragos et al. (2013a) Fragos T., et al., 2013a, ApJ, 764, 41
  * Fragos et al. (2013b) Fragos T., Lehmer B. D., Naoz S., Zezas A., Basu-Zych A., 2013b, ApJ, 776, L31
  * Furtak et al. (2023) Furtak L. J., et al., 2023, arXiv e-prints, p. arXiv:2308.05735
  * Gessey-Jones et al. (2022) Gessey-Jones T., et al., 2022, MNRAS, 516, 841
  * Gnedin et al. (2008) Gnedin N. Y., Kravtsov A. V., Chen H.-W., 2008, ApJ, 672, 765
  * Gouttenoire et al. (2023) Gouttenoire Y., Trifinopoulos S., Valogiannis G., Vanvlasselaer M., 2023, arXiv e-prints, p. arXiv:2307.01457
  * Greene et al. (2023) Greene J. E., et al., 2023, arXiv e-prints, p. arXiv:2309.05714
  * Harikane et al. (2022) Harikane Y., et al., 2022, ApJS, 259, 20
  * Harikane et al. (2023a) Harikane Y., et al., 2023a, arXiv e-prints, p. arXiv:2303.11946
  * Harikane et al. (2023b) Harikane Y., et al., 2023b, ApJS, 265, 5
  * Hasinger (2020) Hasinger G., 2020, J. Cosmology Astropart. Phys., 2020, 022
  * Hawking (1971) Hawking S., 1971, Monthly Notices of the Royal Astronomical Society, 152, 75
  * Heger et al. (2003) Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H., 2003, ApJ, 591, 288
  * Helgason et al. (2016) Helgason K., Ricotti M., Kashlinsky A., Bromm V., 2016, MNRAS, 455, 282
  * Husain et al. (2021) Husain R., Liu B., Bromm V., 2021, MNRAS, 508, 2169
  * Inayoshi et al. (2020) Inayoshi K., Visbal E., Haiman Z., 2020, ARA&A, 58, 27
  * Inman & Ali-Haïmoud (2019) Inman D., Ali-Haïmoud Y., 2019, Phys. Rev. D, 100, 083528
  * Inoue & Kusenko (2017) Inoue Y., Kusenko A., 2017, J. Cosmology Astropart. Phys., 2017, 034
  * Jaacks et al. (2018) Jaacks J., Finkelstein S. L., Bromm V., 2018, MNRAS, 475, 3883
  * Jeon et al. (2014) Jeon M., Pawlik A. H., Bromm V., Milosavljević M., 2014, MNRAS, 440, 3778
  * Jeon et al. (2022) Jeon J., Bromm V., Finkelstein S. L., 2022, MNRAS, 515, 5568
  * Kashlinsky (2016) Kashlinsky A., 2016, ApJ, 823, L25
  * Kashlinsky et al. (2015) Kashlinsky A., Mather J. C., Helgason K., Arendt R. G., Bromm V., Moseley S. H., 2015, ApJ, 804, 99
  * Kashlinsky et al. (2018) Kashlinsky A., Arendt R. G., Atrio-Barandela F., Cappelluti N., Ferrara A., Hasinger G., 2018, Reviews of Modern Physics, 90, 025006
  * Khaire et al. (2016) Khaire V., Srianand R., Choudhury T. R., Gaikwad P., 2016, MNRAS, 457, 4051
  * Klessen & Glover (2023) Klessen R. S., Glover S. C. O., 2023, ARA&A, 61, 65
  * Kocevski et al. (2023) Kocevski D. D., et al., 2023, ApJ, 954, L4
  * Kühnel & Freese (2017) Kühnel F., Freese K., 2017, Phys. Rev. D, 95, 083508
  * Labbé et al. (2023) Labbé I., et al., 2023, Nature, 616, 266
  * Larson et al. (2023) Larson R. L., et al., 2023, ApJ, 953, L29
  * Liu & Bromm (2020a) Liu B., Bromm V., 2020a, MNRAS, 497, 2839
  * Liu & Bromm (2020b) Liu B., Bromm V., 2020b, ApJ, 903, L40
  * Liu & Bromm (2022) Liu B., Bromm V., 2022, ApJ, 937, L30
  * Liu et al. (2021) Liu B., Meynet G., Bromm V., 2021, MNRAS, 501, 643
  * Liu et al. (2022) Liu B., Zhang S., Bromm V., 2022, MNRAS, 514, 2376
  * Liu et al. (2023) Liu B., Sartorio N. S., Izzard R. G., Fialkov A., 2023, arXiv e-prints, p. arXiv:2308.06154
  * Mack et al. (2007) Mack K. J., Ostriker J. P., Ricotti M., 2007, ApJ, 665, 1277
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  * Madau et al. (1997) Madau P., Meiksin A., Rees M. J., 1997, ApJ, 475, 429
  * Madau et al. (1999) Madau P., Haardt F., Rees M. J., 1999, ApJ, 514, 648
  * Maiolino et al. (2023) Maiolino R., et al., 2023, arXiv e-prints, p. arXiv:2308.01230
  * Matthee et al. (2023) Matthee J., et al., 2023, arXiv e-prints, p. arXiv:2306.05448
  * McBride et al. (2009) McBride J., Fakhouri O., Ma C.-P., 2009, MNRAS, 398, 1858
  * Mena et al. (2019) Mena O., Palomares-Ruiz S., Villanueva-Domingo P., Witte S. J., 2019, Phys. Rev. D, 100, 043540
  * Meszaros (1974) Meszaros P., 1974, A&A, 37, 225
  * Mirabel et al. (2011) Mirabel I. F., Dijkstra M., Laurent P., Loeb A., Pritchard J. R., 2011, A&A, 528, A149
  * Mirocha & Furlanetto (2019) Mirocha J., Furlanetto S. R., 2019, MNRAS, 483, 1980
  * Mo et al. (2010) Mo H., Van den Bosch F., White S., 2010, Galaxy formation and evolution. Cambridge University Press
  * Narayan & Yi (1995) Narayan R., Yi I., 1995, ApJ, 452, 710
  * Natarajan et al. (2017) Natarajan P., Pacucci F., Ferrara A., Agarwal B., Ricarte A., Zackrisson E., Cappelluti N., 2017, ApJ, 838, 117
  * Negri & Volonteri (2017) Negri A., Volonteri M., 2017, MNRAS, 467, 3475
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Poulin et al. (2017) Poulin V., Serpico P. D., Calore F., Clesse S., Kohri K., 2017, Phys. Rev. D, 96, 083524
  * Power et al. (2009) Power C., Wynn G. A., Combet C., Wilkinson M. I., 2009, Monthly Notices of the Royal Astronomical Society, 395, 1146
  * Press & Schechter (1974) Press W. H., Schechter P., 1974, ApJ, 187, 425
  * Pringle (1981) Pringle J. E., 1981, ARA&A, 19, 137
  * Pritchard & Loeb (2012) Pritchard J. R., Loeb A., 2012, Reports on Progress in Physics, 75, 086901
  * Ricotti (2007) Ricotti M., 2007, ApJ, 662, 53
  * Ricotti et al. (2008) Ricotti M., Ostriker J. P., Mack K. J., 2008, ApJ, 680, 829
  * Robertson et al. (2010) Robertson B. E., Ellis R. S., Dunlop J. S., McLure R. J., Stark D. P., 2010, Nature, 468, 49
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Sana et al. (2013) Sana H., et al., 2013, A&A, 550, A107
  * Sartorio et al. (2023) Sartorio N. S., et al., 2023, MNRAS,
  * Schauer et al. (2019) Schauer A. T. P., Liu B., Bromm V., 2019, ApJ, 877, L5
  * Serpico et al. (2020) Serpico P. D., Poulin V., Inman D., Kohri K., 2020, Physical Review Research, 2, 023204
  * Sidoli & Paizis (2018) Sidoli L., Paizis A., 2018, MNRAS, 481, 2779
  * Singh et al. (2022) Singh S., et al., 2022, Nature Astronomy, 6, 607
  * Smith & Bromm (2019) Smith A., Bromm V., 2019, Contemporary Physics, 60, 111
  * Smith et al. (2019) Smith A., Ma X., Bromm V., Finkelstein S. L., Hopkins P. F., Faucher-Giguère C.-A., Kereš D., 2019, MNRAS, 484, 39
  * Stacy & Bromm (2013) Stacy A., Bromm V., 2013, MNRAS, 433, 1094
  * Stacy et al. (2011) Stacy A., Bromm V., Loeb A., 2011, ApJ, 730, L1
  * Sureda et al. (2021) Sureda J., Magaña J., Araya I. J., Padilla N. D., 2021, MNRAS, 507, 4804
  * Takhistov et al. (2022) Takhistov V., Lu P., Gelmini G. B., Hayashi K., Inoue Y., Kusenko A., 2022, J. Cosmology Astropart. Phys., 2022, 017
  * Tremmel et al. (2017) Tremmel M., Karcher M., Governato F., Volonteri M., Quinn T. R., Pontzen A., Anderson L., Bellovary J., 2017, MNRAS, 470, 1121
  * Trenti & Stiavelli (2009) Trenti M., Stiavelli M., 2009, ApJ, 694, 879
  * Tristram Acuña & Tseng (2023) Tristram Acuña J., Tseng P.-Y., 2023, arXiv e-prints, p. arXiv:2304.10084
  * Tseliakhovich & Hirata (2010) Tseliakhovich D., Hirata C., 2010, Phys. Rev. D, 82, 083520
  * Vaskonen & Veermäe (2020) Vaskonen V., Veermäe H., 2020, Phys. Rev. D, 101, 043015
  * Wang et al. (2021) Wang F., et al., 2021, ApJ, 907, L1
  * Weltman et al. (2020) Weltman A., et al., 2020, Publ. Astron. Soc. Australia, 37, e002
  * Wilms et al. (2000) Wilms J., Allen A., McCray R., 2000, ApJ, 542, 914
  * Wouthuysen (1952) Wouthuysen S. A., 1952, AJ, 57, 31
  * Yuan et al. (2023) Yuan G.-W., et al., 2023, arXiv e-prints, p. arXiv:2303.09391
  * Zel’dovich & Novikov (1967) Zel’dovich Y. B., Novikov I. D., 1967, Soviet Ast., 10, 602
  * Zhang et al. (2023) Zhang H., Behroozi P., Volonteri M., Silk J., Fan X., Hopkins P. F., Yang J., Aird J., 2023, MNRAS, 518, 2123
  * Ziparo et al. (2022) Ziparo F., Gallerani S., Ferrara A., Vito F., 2022, MNRAS, 517, 1086

##  Appendix A Halo Fitting Formulae

The total stellar mass, as a function of redshift and halo mass, is provided
by Behroozi et al. (2013, 2019); Behroozi et al. (2020); Zhang et al. (2023),
based on a parametric fit to observational data. In this section, we outline
how we utilize their model to calculate the total SRBH and SMBH mass within a
halo.

Initially, we apply the parametric fit for a single halo accretion history,
following the model detailed in the appendix of Behroozi et al. (2013). This
model is initially based on the growth trajectory of the progenitor of a
$10^{13}\rm\ M_{\odot}$ halo, $M_{13}(z)$, and scaled to the haloes of
different masses with the parametric fit. For a halo of mass $M_{0}$ (in units
of $\rm\ M_{\odot}$) at redshift $z=0$, its progenitor mass $M_{h}$ at
redshift $z$ is given by:

$\displaystyle M_{h}\left(M_{0},z\right)$
$\displaystyle=M_{13}(z)10^{f\left(M_{0},z\right)},$ (49) $\displaystyle
M_{13}(z)$ $\displaystyle=10^{13.276}\rm
M_{\odot}(1+z)^{3.00}\left(1+\frac{z}{2}\right)^{-6.11}\exp(-0.503z),$
$\displaystyle f\left(M_{0},z\right)$
$\displaystyle=\log_{10}\left(\frac{M_{0}}{M_{13}(0)}\right)\frac{g\left(M_{0},1\right)}{g\left(M_{0},\frac{1}{1+z}\right)},$
$\displaystyle g\left(M_{0},z\right)$
$\displaystyle=1+\exp\left(-4.651\left(\frac{1}{1+z}-a_{0}\left(M_{0}\right)\right)\right),$
$\displaystyle a_{0}\left(M_{0}\right)$
$\displaystyle=0.205-\log_{10}\left[\left(\frac{10^{9.649}\rm
M_{\odot}}{M_{0}}\right)^{0.18}+1\right].$

For the remainder of the parametric model, we adopt the fitted parameters from
Zhang et al. (2023) to ensure consistency. The star formation rate (SFR, in
units of $\mathrm{\,M}_{\odot}\mathrm{yr}^{-1}$) as a function of halo mass
and redshift is fitted by a double power law

$\displaystyle\mathrm{SFR}_{\mathrm{SF}}\left(M_{h},a,V_{h},\epsilon_{h},\alpha_{h},\beta_{h}\right)$
$\displaystyle=\frac{\epsilon_{h}}{v_{h}^{\alpha_{h}}+v_{h}^{\beta_{h}}}\rm
M_{\odot}\mathrm{yr}^{-1},$ (50) $\displaystyle v_{h}$
$\displaystyle=\frac{v_{\mathrm{Mpeak}}}{V_{h}\cdot\mathrm{km}\mathrm{s}^{-1}},$

where $v_{\mathrm{Mpeak}}$ is the maximum velocity within a halo of mass
$M_{h}$ at redshift $z$:

$\displaystyle v_{\mathrm{Mpeak}}\left(M_{h},a\right)$
$\displaystyle=200\mathrm{~{}km}\mathrm{~{}s}^{-1}\left[\frac{M_{h}}{M_{200\mathrm{kms}}(a)}\right]^{1/3},$
(51) $\displaystyle M_{200\mathrm{kms}}(a)$ $\displaystyle=\frac{1.64\times
10^{12}\mathrm{M}_{\odot}}{\left(\frac{a}{0.378}\right)^{-0.142}+\left(\frac{a}{0.378}\right)^{-1.79}}.$

Given that we have also assumed a halo growth history, this $M_{h}$ will
change with respect to the redshift and eventually evolve into a halo of mass
$M_{0}$ at $z=0$. The remaining parameters
($V_{h},\epsilon_{h},\alpha_{h},\beta_{h}$) can be fitted with the following
values:

$\displaystyle\log_{10}(V_{h})=$ $\displaystyle
2.289+1.548(a-1)+1.218\ln(1+z)-0.087z,$ (52)
$\displaystyle\log_{10}(\epsilon_{h})=$ $\displaystyle
0.556-0.944(a-1)-0.042\ln(1+z)+0.418z,$ $\displaystyle\alpha_{h}=$
$\displaystyle-3.907+32.223(a-1)+20.241\ln(1+z)-2.193z,$
$\displaystyle\beta_{h}=$ $\displaystyle 0.329+2.342(a-1)+0.492z.$

Here, $a=1/(1+z)$ is the scale factor.

Figure 8: Star formation rate (SFR), in units of $\rm\
M_{\odot}\mathrm{\,yr}^{-1}$, as a function of redshift and host halo mass,
from the parametric fit generated by the UNIVERSEMACHINE (Zhang et al., 2023).
Here, in addition to the overall star formation rate, we add a Pop III
component as the dominant star formation process at high redshifts ($z\gtrsim
13$)(Liu & Bromm, 2020a). Similar to Figure 3, the grey region represents the
part of parameter space for the growth trajectories of $\gtrsim 10^{16}\ \rm
M_{\odot}$ haloes at $z=0$, that are not expected to exist in the observable
Universe due to low number densities. White lines represent halo growth
trajectories for a single halo.

The resulting SFR as a function of redshift and current halo mass at redshift
$z$ is plotted as a color map in Figure 8. Since this parametric fit only
provides the SFR at low redshift ($z=0-10$), which corresponds to the SFR of
higher metallicity stars, we also incorporate the SFR for Pop III stars from
Liu & Bromm (2020a), using Equ. (23-24). We find that it generally aligns with
the JWST observation of star-forming galaxies at redshift $z\simeq 10-16$ in
terms of SFR (Harikane et al., 2023b). For higher redshifts $z\gtrsim 10$, due
to the lack of data, we extrapolate the fitted star formation rate from the
UNIVERSEMACHINE to approximate the stellar mass.

To calculate the SMBH mass in a halo of mass $M_{\rm h}$ at redshift
$z=1/a-1$, we require a few additional steps. Firstly, we calculate the total
stellar mass $M_{\star}(M_{h},a)$ as

$\displaystyle\log_{10}\left(M_{\star}\right)$
$\displaystyle=\log_{10}\left(\epsilon_{\star}M_{1}\right)+f\left(\log_{10}\left(\frac{M_{h}}{M_{1}}\right)\right)-f(0),$
(53) $\displaystyle f(x)$
$\displaystyle=-\log_{10}\left(10^{\alpha_{\star}x}+1\right)+\delta_{\star}\frac{\left(\log_{10}(1+\exp(x))\right)^{\gamma_{\star}}}{1+\exp\left(10^{-x}\right)},$

with fitted parameters
($v_{\star},\epsilon_{\star},M_{1},\alpha_{\star},\delta_{\star},\gamma_{\star}$):

$\displaystyle v_{\star}=$ $\displaystyle\exp\left(-4a^{2}\right),$ (54)
$\displaystyle\log_{10}(\epsilon_{\star})=$
$\displaystyle-1.777-0.006(a-1)v_{\star}-0.119(a-1),$
$\displaystyle\log_{10}\left(M_{1}\right)=$ $\displaystyle
11.514+\left(-1.793(a-1)-0.251z\right)v_{\star},$
$\displaystyle\alpha_{\star}=$ $\displaystyle-1.412+0.731(a-1)v_{\star},$
$\displaystyle\delta_{\star}=$ $\displaystyle
3.508+\left(2.608(a-1)-0.043z\right)v_{\star},$ $\displaystyle\gamma_{\star}=$
$\displaystyle 0.316+\left(1.319(a-1)+0.279z\right)v_{\star}.$

Once we have the stellar mass for a halo at redshift $z$, the corresponding
bulge mass $M_{\text{bulge }}$ can be calculated using the following
equations:

$\displaystyle M_{\text{bulge }}$
$\displaystyle=\frac{f_{z}(z)M_{\star}}{1+\exp\left\\{-1.13\left[\log_{10}\left(M_{\star}/10^{10.2}\rm
M_{\odot}\right)\right]\right\\}},$ (55) $\displaystyle f_{z}(z)$
$\displaystyle=\frac{z+2}{2z+2}.$ (56)

Finally, $\tilde{M}_{\mathrm{SMBH}}\left(M_{\text{bulge
}},\gamma_{\bullet},\beta_{\bullet}\right)$, the average SMBH mass
corresponding to the bulge mass $M_{\text{bulge }}$ is determined by the
following fitting formula:

$\displaystyle\log_{10}\tilde{M}_{\mathrm{SMBH}}=\beta_{\bullet}+\gamma_{\bullet}\log_{10}\left(\frac{M_{\text{bulge
}}}{10^{11}\rm\ M_{\odot}}\right),$ (57)

where

$\displaystyle\gamma_{\bullet}=$ $\displaystyle 1.028+0.036(a-1)+0.052z,$ (58)
$\displaystyle\beta_{\bullet}=$ $\displaystyle 8.343-0.173(a-1)+0.044z.$

The resulting plot of the average SMBH mass in a given host halo of mass
$M_{h}$ at redshift $z$ is given in Figure 9. When compared with the SRBH mass
abundance within a halo, we find that it is smaller by at least one order of
magnitude in higher mass haloes at $z\gtrsim 6$. However, the higher accretion
rate for SMBHs will compensate for this difference by generating more energy
feedback while growing in mass.

Figure 9: 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^{7}\rm\ M_{\odot}$.
Similar to Figure 3, the grey area represents the combination of parameters
where no halo formation is allowed, and white lines represent individual halo
growth histories.

##  Appendix B Accreting BH Spectra

Next, we summarize the prescription for the radiation feedback spectrum from
BH accretion, as described in Takhistov et al. (2022). To streamline the
analysis, we categorize the accretion process into two scenarios based on the
normalized accretion factor $\dot{m}$: the thin-disk regime and the advection-
dominated accretion flow (ADAF) regime.

For accretion in the IGM, where we have calculated a low accretion rate for
the BHs of interest, we exclusively apply the ADAF model. Conversely, for
SMBHs and HMXBs accreting at the Eddington rate, we assume a thin disk regime,
where the blackbody spectrum is dominant. For other cases, we employ a
combination of the two regimes depending on the BH accretion rate and the
viscosity parameter $\alpha=0.1$.

Thin disk ($\dot{m}\geq 0.07\alpha$)

In the thin disk regime, the temperature profile of the accretion disk $T(r)$
can be modeled as a function of the radius (See also Pringle, 1981):

$T(r)=T_{i}\left(\frac{r_{i}}{r}\right)^{3/4}\left[1-\left(\frac{r_{i}}{r}\right)^{\frac{1}{2}}\right]^{\frac{1}{4}}.$
(59)

Here, $r_{i}$ is the inner radius of the disk, which is equal to three times
the Schwarzschild radius $r_{s}$. The temperature $T_{i}$ at the inner radius
is derived from the balance between gravitational energy release and radiative
cooling:

$T_{\rm i}=\left(\frac{3GM_{\bullet}\dot{M}_{\bullet}}{8\pi
r_{i}^{3}\sigma_{B}}\right)^{1/4}=1.3\times 10^{3}\
\mathrm{eV}/k\,\dot{m}^{\frac{1}{4}}\left(\frac{M_{\bullet}}{10\mathrm{M}_{\odot}}\right)^{-\frac{1}{4}}.$
(60)

The maximum temperature $T_{\rm max}=T(r_{\rm max})$ of the disk is attained
at $r_{\rm max}=1.36r_{i}$, as given by the temperature profile. The minimum
temperature $T_{\rm o}=T(r_{\rm out})$ is reached at the outer radius of the
disk:

$r_{\text{out }}\simeq 1.16\times 10^{10}r_{\rm
s}\left(\frac{M_{\bullet}}{10\mathrm{M}_{\odot}}\right)^{\frac{2}{3}}\left(\frac{\tilde{v}}{10\mathrm{~{}km}/\mathrm{s}}\right)^{-\frac{10}{3}}.$
(61)

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.

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:

$L_{\nu}=\begin{cases}c_{\alpha}\left(\frac{T_{\max}}{T_{o}}\right)^{\frac{5}{3}}\left(\frac{h\nu}{kT_{\max}}\right)^{2}&(\nu<kT_{o}/h)\\\
c_{\alpha}\left(\frac{h\nu}{kT_{\max}}\right)^{\frac{1}{3}}&(kT_{o}/h<\nu<kT_{\max}/h)\\\
c_{\alpha}\left(\frac{h\nu}{kT_{\max}}\right)^{2}e^{1-\frac{h\nu}{kT_{\max}}}&(\nu>kT_{\max}/h)\\\
\end{cases}.$ (62)

Here, $c_{\alpha}$ is the normalizing flux constant such that the integration
over all frequencies gives a total power of $0.057\dot{M}_{\bullet}c^{2}$,
which recovers the subgrid model (Tremmel et al., 2017).

ADAF ($\dot{m}<0.07\alpha$):

The ADAF disk dominates when the accretion rate $\dot{m}$ is significantly
sub-Eddington888For a complete derivation, see e.g. Narayan & Yi (1995).. In
this regime, we primarily consider the spectrum from the sychrotron radiation
and Inverse-Compton (IC) scattering process:

$L_{\nu}=\begin{cases}L_{\nu_{p}}\left(\frac{\nu}{\nu_{p}}\right)^{5/2}&(\text{Sychrotron},\
\nu<\nu_{p})\\\
L_{\nu_{p}}\left(\frac{\nu}{\nu_{p}}\right)^{-\alpha_{c}}&(\text{IC},\
\nu_{p}\leq\nu\lesssim 3kT_{e}/h)\end{cases}.$ (63)

where the peak frequency $\nu_{p}$ is

$\displaystyle\nu_{p}=$ $\displaystyle 1.83\times 10^{3}\mathrm{eV}/h$ (64)
$\displaystyle\dot{m}^{\frac{3}{4}}\left(\frac{\alpha}{0.1}\right)^{-\frac{1}{2}}\left(\frac{1-\beta}{1/11}\right)^{\frac{1}{2}}\theta_{e}^{2}\left(\frac{r_{i}}{3r_{s}}\right)^{-\frac{5}{4}}\left(\frac{M_{\bullet}}{10\mathrm{M}_{\odot}}\right)^{-\frac{1}{2}},$

and the corresponding peak luminosity:

$L_{\nu_{p}}=5.06\times
10^{37}\frac{\mathrm{erg}}{\mathrm{s}\cdot\mathrm{eV}/h}\dot{m}^{\frac{3}{2}}\alpha^{-1}(1-\beta)\theta_{e}^{5}\left(\frac{r_{i}}{r_{s}}\right)^{-\frac{1}{2}}\left(\frac{M_{\bullet}}{10\mathrm{M}_{\odot}}\right).$
(65)

To ensure overall consistency, we normalize the product of the peak frequency
and luminosity, $\nu_{p}L_{\nu_{p}}$, with the results from Equ. (44).

In these equations, $\beta$ is a partition parameter for total pressure
originating from the gas and magnetic field. with values equaling $10/11$.
$\theta_{e}$ is the dimensionless electron temperature written in units of
electron rest mass $m_{e}$:

$\displaystyle\theta_{e}=$
$\displaystyle\frac{kT_{e}}{m_{e}c^{2}}=\frac{T_{e}}{5.93\times 10^{9}K}$ (66)
$\displaystyle\simeq$ $\displaystyle
0.17A_{c}^{-\frac{1}{7}}\delta^{\frac{1}{7}}\alpha^{\frac{3}{14}}(1-\beta)^{-\frac{1}{14}}\left(\frac{r_{i}}{r_{s}}\right)^{\frac{3}{28}}\left(\frac{M_{\bullet}}{10\mathrm{M}_{\odot}}\right)^{\frac{1}{14}}\dot{m}^{-\frac{5}{28}},$

where $A_{c}^{-\frac{1}{7}}$ is a pre-factor related to the relative
contribution of IC to the total power depending on the accretion rate. Since
we do not further divide the ADAF case, we will assume a uniform value of 1.1.
$\delta$ is the electron heating efficiency approximated by a value of $0.3$.

The exponent $\alpha_{c}$ for the IC part of the spectrum in Equ. (63), is
determined by the following expression:

$\alpha_{c}=-\frac{\ln\tau_{\rm es}}{\ln A},$ (67)

where $\tau_{\rm es}=12.4\dot{m}\alpha^{-1}(r_{i}/r_{s})^{-1/2}$ is the
optical depth for electron scattering, and $A=1+4\theta_{e}+16\theta_{e}^{2}$
is the amplification factor. Typically, for a $10\rm\,M_{\odot}$ BH that
accretes at a rate of $\dot{m}\sim 10^{-6}$, this exponent will yield a value
of $\alpha_{c}\sim 2$, where the slope becomes asymptotically flatter with
increasing accretion rate $\dot{m}$ and BH mass $M_{\bullet}$. When the
accretion is efficient, reaching the critical limit for the thin disk of
regime $\dot{m}\sim 10^{-2}$, for the same BH mass of $10\rm\,M_{\odot}$, the
power-law exponent approaches $\alpha_{c}\sim 1$.