© The Authors 2025

1. Introduction

Cosmic radiation backgrounds serve as integral probes of the energetic and dynamic processes that shaped the early Universe. These diffuse emissions, spanning from radio to high-energy X-rays, encapsulate the integrated contributions of astrophysical sources over cosmic time. Notably, the cosmic X-ray background (CXB), cosmic radio background (CRB), and the Lyman-Werner background (LWB) provide critical insights into energetic astrophysical processes, the interstellar medium (ISM), and the regulation of star formation in the early Universe. Observations have reported significant excesses in the CXB and CRB beyond the expected contributions from known astrophysical sources (Fixsen et al. 2011; Cappelluti et al. 2017). At the same time, the LWB plays a crucial role in regulating star formation by dissociating molecular hydrogen (H₂), the primary coolant in metal-poor gas during the early universe (Stecher & Williams 1967; Saslaw & Zipoy 1967; Peebles & Dicke 1968).

The CXB, primarily measured by satellite missions such as the Chandra X-ray Observatory, the Röntgen Satellite (ROSAT), and the X-ray Multi-Mirror Mission (XMM-Newton), traces high-energy processes associated with active galactic nuclei (AGNs), star-forming galaxies, and high-redshift black holes (Gilli et al. 2007). Although most of the CXB can be attributed to known populations, observations reveal an unresolved excess in the 0.5–10 keV band. This excess is quantified as non-source contributions (nsCXB), defined as the component of the CXB that remains after masking all sources detected in X-ray surveys and those with counterparts identified in other wavelengths. nsCXB accounts for $9.7_{-1.8}^{+1.6}$ per cent of the total CXB in the 0.5–2 keV range and ${17}_{-7.0}^{+5.9}$ per cent in the 2–10 keV range (Cappelluti et al. 2017)¹. Potential contributors include heavily obscured, Compton-thick AGNs, which still avoid detection with traditional X-ray surveys, and star-forming galaxies at high redshifts. However, their contribution is estimated to be minor. Population synthesis models, which combine AGN luminosity functions and spectral properties, suggest that these sources can only partially explain the nsCXB, leaving a significant fraction unresolved (Ananna et al. 2020). Ongoing and upcoming missions, such as the Nuclear Spectroscopic Telescope Array (NuSTAR; Harrison et al. 2013) and the High-Energy X-ray Probe (HEX-P; Civano et al. 2024), are expected to provide sensitive broad-band X-ray observations, helping to refine these models and uncover the nature of the remaining unresolved sources.

The CRB observed using instruments such as the Absolute Radiometer for Cosmology, Astrophysics, and Diffuse Emission (ARCADE 2) is primarily composed of contributions from well-known astrophysical sources, including the cosmic microwave background (CMB), radio-loud AGNs, and star-forming galaxies. However, ARCADE 2 measurements reveal a puzzling excess at frequencies between 3 and 8 GHz, with a reported temperature rise of 54 ± 6 mK at 3.3 GHz (Fixsen et al. 2011). This excess remains unexplained by known sources, leading to exploring alternative explanations. Proposed hypotheses include synchrotron radiation from highly efficient, accreting supermassive black holes (SMBHs; Ewall-Wice et al. 2018) and extreme populations of supernova remnants (Mirocha & Furlanetto 2019). More exotic processes have also been suggested, such as synchrotron radiation emitted by relativistic electrons through late decays of a metastable particle (Cline & Vincent 2013), decays of dark matter (DM) particles resulting in dark photons oscillating into ordinary photons (Caputo et al. 2023), and radiative decays of relic neutrinos into sterile neutrinos (Bhupal Dev et al. 2024). However, these scenarios require extreme efficiencies or unconventional physics that remain challenging to reconcile with current observations. Further investigations into the origins of the CRB excess are necessary to fully understand its implications for cosmology and the physics of diffuse backgrounds.

The LWB, consisting of ultraviolet photons with energies between 11.2 eV and 13.6 eV, plays a pivotal role in shaping the early Universe by regulating star formation. These photons, mainly emitted by massive stars, dissociate molecular hydrogen (H₂) through the Solomon process (Stecher & Williams 1967), suppressing the primary cooling mechanism in metal-poor gas (Saslaw & Zipoy 1967; Peebles & Dicke 1968). Molecular cooling becomes inefficient above temperatures of ∼10⁴ K, where atomic cooling processes dominate, particularly in halos with sufficient gravitational potential to sustain such mechanisms (Haiman et al. 1996; Machacek et al. 2001). The intensity of the LWB and its interaction with self-shielding effects in dense regions determine the extent of this suppression (e.g. Johnson et al. 2008). In regions where self-shielding is effective, H₂ molecules can survive, allowing localised cooling and star formation to proceed (Wolcott-Green et al. 2011; Hartwig et al. 2015). This typically occurs in halos with masses above ∼10⁶ − 10⁷ M_⊙, where higher gas densities enable efficient self-shielding against the dissociating LW radiation. However, in diffuse environments, such as lower-mass halos below this threshold, the LWB efficiently dissociates molecular hydrogen, delaying or suppressing star formation. This feedback mechanism increases the critical halo mass threshold for star formation, governing the chemothermal evolution of the primordial gas (Latif & Khochfar 2019; Schauer et al. 2021). Beyond its role in star formation, the LWB sets the conditions necessary for the direct collapse of gas into black hole seeds, which may serve as precursors to SMBHs observed at high redshifts (Agarwal et al. 2019; Wise et al. 2019).

Primordial black holes (PBHs) have emerged as a compelling hypothesis to explain the observed anomalies in the CXB and CRB, as well as potential contributors to the LWB. These objects, theorised to form from density fluctuations in the early Universe (Carr & Hawking 1974), have been proposed as candidates for DM and could constitute a significant fraction of its total density (Carr et al. 2016; Khlopov 2010). Unlike stellar black holes, PBHs are not limited by the processes of stellar collapse and can span a wide range of masses (Carr 1975), offering unique insights into early universe physics. Their accretion processes produce radiation across the electromagnetic spectrum, including the X-ray, radio, and ultraviolet bands bands (Carr 1979; Ricotti et al. 2008; Mack et al. 2007). Specifically, accreting PBHs could generate X-rays and synchrotron emission from relativistic electrons, partially accounting for the observed excesses in the CXB and CRB (Hasinger 2020; Cappelluti et al. 2022; Mittal & Kulkarni 2022). In addition, PBHs could influence the LWB by emitting ultraviolet radiation during gas accretion, further regulating molecular hydrogen dissociation and impacting star formation in low-mass halos.

Recent James Webb Space Telescope (JWST) discoveries have significantly expanded our understanding of the early Universe, revealing numerous massive galaxies and supermassive black holes (SMBHs) at high redshifts. Many galaxy candidates at z ≳ 10 have been identified using NIRCam imaging through programmes such as Early Release Observations (ERO; Pontoppidan et al. 2022) and Director’s Discretionary Early Release Science (ERS). Examples include GN-z11, identified at z = 10.6 (Tacchella et al. 2023), which exhibits a higher-than-expected stellar mass. Such discoveries, if correct, challenge the standard ΛCDM cosmology by highlighting the rapid formation of massive structures in the early universe (e.g. Atek et al. 2022; Castellano et al. 2022; Finkelstein et al. 2022; Naidu et al. 2022; Yan et al. 2022; Bradley et al. 2023; Robertson et al. 2023; Yan et al. 2023). Similarly, SMBHs with masses exceeding 10⁷M_⊙ have been observed at z = 7 − 10, suggesting growth rates that are difficult to reconcile with traditional black hole formation models (Goulding et al. 2023; Larson et al. 2023; Bogdán et al. 2024; Greene et al. 2024; Natarajan et al. 2023). These findings suggest that PBHs could act as seeds for SMBHs, accelerating their growth and contributing to the rapid formation of massive galaxies in the early universe (Carr & Kühnel 2020). An extreme example of such a BH-dominated high-redshift system discovered by JWST is Abell 2744-QSO1, which is overmassive even with respect to the host halo dynamical mass (as opposed to its stellar mass), challenging standard formation channels (Ji et al. 2025).

Simulations further support this hypothesis. Colazo et al. (2024) demonstrated that even a small abundance of PBHs, accounting for only 0.5 per cent of the DM, with extended mass functions from Sureda et al. (2021), can provide the seeds for SMBHs and enhance the abundance of high-redshift galaxies observed by JWST. Their study, employing the SWIFT simulation code, highlights how PBHs can drive structure formation in the early Universe with realistic star formation efficiencies. Additionally, Liu & Bromm (2022) showed, through an analytical model based on linear perturbation theory and the Press-Schechter (PS) formalism, that massive PBHs (≳10⁹ M_⊙), assuming a monochromatic mass function and comprising ∼10⁻⁶ − 10⁻³ of the DM, can significantly accelerate the formation of massive galaxies at z ≳ 10. Their analysis demonstrated that the gravitational potential of such PBHs enhances structure formation by reducing the star formation efficiency required to explain the observed abundance of massive galaxies at high redshifts, aligning theoretical predictions with JWST observations. In addition, Liu et al. (2022), using cosmological hydrodynamic zoom-in simulations and semi-analytical models, investigated the impact of stellar-mass PBHs on the first generation of stars, showing that PBHs amplify small-scale density perturbations and influence early star formation through accretion feedback. These models not only align with observational constraints on PBH parameters but also predict potential signatures that JWST could confirm in future studies, offering a promising avenue for testing the role of PBHs in cosmic evolution. Collectively, these findings underscore the ability of PBHs, even in minimal fractions, to serve as seeds for SMBHs and play a pivotal role in shaping the formation and evolution of massive structures in the early Universe. While PBHs are a compelling candidate to explain the observed abundance of massive galaxies and SMBHs at high redshifts, other hypotheses have also been proposed. For instance, Koehler et al. (2024) explored the role of cosmic string loops in seeding galaxy formation, demonstrating through cosmological simulations that they could also contribute to explaining JWST observations.

Focusing on cosmic radiation backgrounds, Ziparo et al. (2022) investigated emissions produced by PBHs, particularly in the CXB and CRB. Using the Bondi–Hoyle–Lyttelton accretion model, they calculated the accretion rate based on gas density, sound speed, and the relative velocity between the PBH and the surrounding gas. Their findings indicate that PBHs can produce significant emissions in X-ray and radio bands through accretion. However, their results suggest that while PBHs could account for a fraction of the CXB, they cannot fully explain the CRB excess observed by ARCADE 2.

Similarly, PBHs have been explored as a possible explanation for the anomalous 21 cm absorption signal reported by the Experiment to Detect the Global Epoch of Reionisation Signature (EDGES; Bowman et al. 2018; Ewall-Wice et al. 2018). However, this detection remains debated, and it is important to acknowledge that the Shaped Antenna Measurement of the Background Radio Spectrum (SARAS) experiment found no evidence supporting the signal and suggested that systematic effects might be responsible for the observed feature (Singh et al. 2022). Despite this controversy, the EDGES signal is not definitively ruled out and, therefore, continues to be investigated in theoretical studies, particularly in relation to PBH-driven heating. In this context, Ziparo et al. (2022) showed that PBHs could impact the 21 cm signal by heating the intergalactic medium (IGM) through X-ray emission and by contributing to a potential excess radio background; however, their results indicate that the conditions required for PBHs to reproduce the observed EDGES signal are tightly constrained by existing observational limits.

Our previous work, Casanueva-Villarreal et al. (2024), hereafter CV2024, investigated the role of PBHs in the early Universe by analysing their impact on gas properties at very high redshifts (z ∼ 23). This study employed a combination of hydrodynamical simulations, semi-analytical models, and analytical approaches to explore various accretion regimes, including standard advection-dominated accretion flow (ADAF), electron ADAF (eADAF), luminous hot accretion flow (LHAF), and thin disc models. This comprehensive analysis revealed how the temperature and hydrogen abundance near PBH-hosting halos depend on different DM scenarios, offering a more nuanced understanding of PBH accretion processes and their influence on early galaxy formation.

Despite their theoretical potential, PBHs remain heavily constrained by multiple observations. For example, Ziparo et al. (2022) concluded that while PBHs could contribute to the CXB, their model could not fully account for the observed CRB. In addition, various multi-messenger observations have placed stringent bounds on the allowed PBH abundance. These include microlensing constraints from surveys such as Subaru (Niikura et al. 2019a; Kusenko et al. 2020), EROS-2/MACHO (Tisserand et al. 2007), and OGLE (Niikura et al. 2019b; Mróz et al. 2024), gravitational wave limits from LIGO/Virgo (Abbott et al. 2019; Wong et al. 2021), constraints from cmB spectral distortions due to accretion (Poulin et al. 2017), and limits from radio and X-ray observations of the Galactic Centre (Manshanden et al. 2019). Together, these independent constraints span a broad mass range and significantly restrict the parameter space in which PBHs could constitute a substantial fraction of DM (see Carr et al. 2021 for a review). Additionally, while PBHs might play a role in SMBH formation and galaxy assembly, further observational and theoretical refinement is needed to evaluate their overall impact comprehensively; for a detailed discussion, see Liu et al. (2024).

This study aims to constrain the fraction of DM that can be composed of PBHs by examining the radiation backgrounds they produce. Specifically, we compared the results considering different gas density profiles within halos to understand how these profiles affect the predicted radiation and the resulting constraints on the PBH contribution to DM. Building on CV2024, we analysed PBH accretion across various sub-regimes, such as eADAF, ADAF, LHAF, and thin disc models. This approach not only refines constraints on PBH parameters but also sheds light on their potential role as DM candidates and their broader influence on early universe cosmology.

This paper is structured as follows. Section 2 outlines the theoretical framework of this study. In Section 2.1.1, we describe the cosmological distribution of PBHs, including their interplay with the large-scale structure of the Universe. Section 2.2 focuses on modelling accretion processes, highlighting the different regimes and their implications for radiation production. Section 3 presents the main results, evaluating the contributions of PBHs to the X-ray, Lyman-Werner, and radio backgrounds and comparing these to current observational constraints. In Section 3.2, we explore the impact of varying key modelling assumptions on our results to assess their robustness. Finally, Section 4 synthesises the findings and outlines prospects for future research.

2. Model

In this section, we elaborate on the theoretical framework employed in our study to examine the effects of PBHs on cosmic radiation backgrounds. Our model incorporates the cosmological distribution of PBHs as well as the physical processes associated with their accretion and the resulting radiation.

We adopted a ΛCDM cosmology consistent with the Planck 2018 results (Planck Collaboration VI 2020): Ω_m = 0.315, Ω_Λ = 0.685, Ω_b = 0.049, σ₈ = 0.811, n_s = 0.965, and H₀ = 67.4 km s⁻¹ Mpc⁻¹.

2.1. Cosmological distribution and structure formation

2.1.1. Cosmological distribution of PBHs

We assume a monochromatic mass distribution of PBHs which constitute a fraction f_PBH of DM, defined as f_PBH = Ω_PBH/Ω_DM. The PBH density distribution follows that of DM, which is organised into virialised objects (halos) and a diffuse component in the IGM. The number density of PBHs at redshift z is given by

$$\begin{array}{c}\hfill n_{\mathrm{PBH}}(z)=\frac{{\mathrm{\Omega }}_{\mathrm{DM}}{\rho }_{\mathrm{c}}{f}_{\mathrm{PBH}}{(1+z)}^3}{M_{\mathrm{PBH}}}=n_{\mathrm{PBH},\mathrm{IGM}}(z)+n_{\mathrm{PBH},\mathrm{halos}}(z),\end{array}$$(1)

where ρ_c = 3H₀²/(8πG) is the critical density of the Universe at z = 0, with G = 6.674 × 10⁻⁸ cm³ g⁻¹ s⁻² being the gravitational constant. M_PBH denotes the monochromatic mass of PBHs, while n_PBH, IGM and n_PBH, halos represent their number densities in the IGM and within DM halos, respectively.

• PBHs in the IGM: The number density of PBHs in IGM can be expressed as

  $$\begin{array}{c}\hfill n_{\mathrm{PBH}}^{\mathrm{IGM}}(z)=\frac{{\mathrm{\Omega }}_{\mathrm{DM}}{\rho }_{\mathrm{c}}{(1+z)}^3(1-f_{\mathrm{coll}}(z))f_{\mathrm{PBH}}}{M_{\mathrm{PBH}}},\end{array}$$(2)

  where f_coll(z) is the fraction of DM that has collapsed into virialised halos by redshift z, calculated using the PS formalism

  $$\begin{array}{c}\hfill f_{\mathrm{coll}}(M_{\min },z)=\mathrm{erfc}\left[\frac{{\delta }_{\mathrm{c}}(z)}{\sqrt{2{\sigma }_M^2}}\right],\end{array}$$(3)

  with δ_c(z) = 1.686/D(z) being the critical density for collapse, D(z)∝(1 + z)⁻¹ the growth factor, and σ_M the standard deviation of the linearly extrapolated matter power spectrum. M_min is the minimum halo mass for efficient cooling processes to occur. For M_min, we followed Liu et al. (2022) and defined it as the maximum between (i) the virial mass corresponding to a virial temperature of 100 K, and (ii) the minimum mass required for a halo to host at least one PBH, given its mass and fractional abundance. The 100 K threshold is motivated by the characteristic relative velocity between baryons and DM at z ≈ 30, which corresponds to a baryonic kinetic temperature of approximately 100 K. This velocity suppresses early gas collapse and sets a natural temperature scale for gas capture into halos. Additionally, PBH accretion can heat the surrounding gas above this temperature, reinforcing its use as a conservative lower limit for the onset of significant PBH feedback in halos. Eq. (2) indicates that $n_{\mathrm{PBH}}^{\mathrm{IGM}}(z)$ decreases with decreasing redshift. Besides cosmic expansion, this decrease is attributed to the ongoing structure formation, where more DM collapses into virialised halos, leading to a higher fraction of PBHs being locked into these structures rather than remaining in the IGM.

• PBHs in DM halos: We assume that the DM distribution within halos follows a Navarro-Frenk-White (NFW) density profile, which is widely used to describe the structure of DM halos. This profile is characterised by a density that decreases with increasing radial distance from the halo centre, described as

  $$\begin{array}{c}\hfill {\rho }_{\mathrm{DM}}(r)=\frac{{\rho }_{\mathrm{c}}{\delta }_{\mathrm{c}}}{cx{(1+cx)}^2},\end{array}$$(4)

  where x = r/r_vir is the radial distance normalised by the virial radius r_vir. The virial radius is expressed as (Barkana & Loeb 2001)

  $$\begin{array}{c}\hfill r_{\mathrm{vir}}=0.784{\left(\frac{M_{\mathrm{vir}}}{{10}^8{h}^{-1}{\mathrm{M}}_{\odot }}\right)}^{1/3}{\left[\frac{{\mathrm{\Omega }}_m}{{\mathrm{\Omega }}_m^z}\frac{{\mathrm{\Delta }}_{\mathrm{c}}}{18{\pi }^2}\right]}^{-1/3}{\left(\frac{1+z}{10}\right)}^{-1}{h}^{-1}\mathrm{kpc},\end{array}$$(5)

  where the overdensity relative to ρ_c at the collapse redshift is given by

  $$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{c}}=18{\pi }^2+82d-39d^2,\end{array}$$(6)

  with d = Ω_m^z − 1 and Ω_m^z = Ω_m(1 + z)³/(Ω_m(1 + z)³ + Ω_Λ). The parameter Δ_c is related to δ_c through

  $$\begin{array}{c}\hfill {\delta }_{\mathrm{c}}=\frac{{\mathrm{\Delta }}_{\mathrm{c}}}{3}\frac{c^3}{ln(1+c)-c/(1+c)}=\frac{{\mathrm{\Delta }}_{\mathrm{c}}}{3}\frac{c^3}{F(c)},\end{array}$$(7)

  where c is the concentration parameter, which encapsulates the relationship between the halo’s characteristic density and its virial mass. This parameter reflects how tightly matter is bound within the halo, varying with both the virial mass M_vir and the redshift. It is expressed as

  $$\begin{array}{c}\hfill logc=1.071-0.098[log\left(\frac{M_{\mathrm{vir}}}{{\mathrm{M}}_{\odot }}\right)-12].\end{array}$$(8)

  This relation is based on the N-body simulations by Macciò et al. (2007), which are calibrated for halos at z = 0. Following Ziparo et al. (2022), we implemented the redshift evolution of the concentration parameter as c ∝ (1 + z)⁻¹, as suggested by previous works (Barkana & Loeb 2002; Duffy et al. 2008; Ricotti 2009). Within the virial radius of the halo, we can determine the number of PBHs enclosed within any given radial distance r, assuming their distribution follows that of the DM. This is expressed as

  $$\begin{array}{c}\hfill N_{\mathrm{PBH}}(r)=f_{\mathrm{PBH}}\frac{4\pi }{M_{\mathrm{PBH}}}{\displaystyle {\int }_0^r}{\rho }_{\mathrm{DM}}(r\prime )r{\prime }^2\mathrm{d}r\prime .\end{array}$$(9)

  This approach allows us to analyse how PBHs are distributed within a DM halo and how their distribution might influence structure formation and radiation emission.

2.1.2. Halo mass function in structure formation with PBHs

To model the Halo Mass Function (HMF) in cosmologies that include PBHs, we employ the formalism developed by Zhang et al. (2024a) (see also Padilla et al. 2021), which extends the standard PS approach (Press & Schechter 1974) to incorporate the effects of PBHs on the power spectrum of density perturbations. This approach modifies the linear power spectrum to include adiabatic perturbations, isocurvature contributions induced by PBHs, and a correlation term that captures the mode mixing between these components.

The total power spectrum, extrapolated to z = 0, is expressed as

$$\begin{array}{c}\hfill P_{\mathrm{tot}}(k)=P_{\mathrm{ad}}(k)+P_{\mathrm{iso}}(k)+P_{\mathrm{corr}}(k),\end{array}$$(10)

where P_ad(k) denotes the standard adiabatic power spectrum in ΛCDM cosmology, computed using the Planck 2018 cosmological parameters (Planck Collaboration VI 2020). The term P_iso(k) represents the contribution from isocurvature perturbations arising due to the discreteness of PBHs. Finally, the correlation term P_corr(k) captures the interaction between adiabatic and isocurvature modes, which is significant at intermediate scales where the influence of PBHs becomes dominant.

The isocurvature contribution P_iso(k) is given by

$$\begin{array}{c}\hfill P_{\mathrm{iso}}(k)=\frac{{[f_{\mathrm{PBH}}{D}_0]}^2}{{\overline{n}}_{\mathrm{PBH}}}.\end{array}$$(11)

Here, ${\overline{n}}_{\mathrm{PBH}}$ represents the average comoving number density of PBHs, while D₀ = D_PBH(a = 1) denotes the growth factor of PBH-induced perturbations at z = 0. This growth factor evolves as a function of the scale factor a and can be well approximated using the analytical formalism presented in Inman & Ali-Haïmoud (2019)

$$\begin{array}{cc}\hfill D_{\mathrm{PBH}}(a)& \approx {(1+\frac{3\gamma a}{2a_{-}-a_{\mathrm{eq}}})}^{a_{-}},\gamma =\frac{{\mathrm{\Omega }}_m-{\mathrm{\Omega }}_{\mathrm{b}}}{{\mathrm{\Omega }}_{\mathrm{m}}},\hfill \\ \hfill a_{-}& =\frac{1}{4}(\sqrt{1+24\gamma }-1),\hfill \end{array}$$(12)

where a_eq ∼ 1/3400 corresponds to the scale factor at the epoch of matter-radiation equality.

The mode-mixing term P_corr(k) is expressed as

$$\begin{array}{c}\hfill P_{\mathrm{corr}}(k)=f_{\mathrm{PBH}}{D}_0^2{\left(\frac{k}{k_{\mathrm{PBH}}}\right)}^3{P}_{\mathrm{ad}}(k),\mathrm{for}k\le 3k_{\mathrm{PBH}},\end{array}$$(13)

and vanishes for k > 3k_PBH. Here, $k_{\mathrm{PBH}}={(2{\pi }^2{\overline{n}}_{\mathrm{PBH}})}^{1/3}$ is the characteristic scale below which isocurvature effects dominate.

Using this modified power spectrum, the HMF is computed via the PBH-modified PS formalism, incorporating corrections for ellipsoidal collapse as described in Sheth & Tormen (1999). While this formalism provides a practical framework for studying PBH cosmologies, it has well-known limitations. Specifically, it underestimates the abundance of small halos due to the suppression of their formation and growth by nearby PBH-seeded halos (Zhang et al. 2024b). At the high-mass end, the PS formalism tends to overpredict the abundance of massive halos by up to a factor of 10 for f_PBH ≳ 0.01, as it does not fully account for the nonlinear suppression of halo growth due to PBH-induced structures and their interactions (Liu et al. 2024). Furthermore, this formalism fails to capture the complex nonlinear effects surrounding individual PBHs, as well as the characteristic bimodal distribution of halos, which emerges from the interplay between isocurvature perturbations and PBH-driven structure formation (Zhang et al. 2024b; Liu et al. 2024). As highlighted in Liu et al. (2024) and Zhang et al. (2024b), addressing these limitations is crucial for improving our understanding of PBH-induced structure formation. While the PS formalism is strictly valid for modelling the Poisson effect in the linear regime, the seed effect is not explicitly included in the HMFs derived from this formalism. Future work could explore high-resolution simulations and refined analytical frameworks to better capture nonlinear dynamics and PBH-environment interactions, which remain an open challenge in PBH cosmologies.

In the top panel of Fig. 1, we show the total power spectrum P_tot(k) for different values of f_PBH, assuming a PBH monochromatic mass of M_PBH = 100 M_⊙. Compared to higher masses, such as M_PBH = 10⁶ M_⊙, which are explored in Zhang et al. (2024b), the enhancement at small scales (k ≳ 10² h Mpc⁻¹) due to the isocurvature component is moderate but still present. This enhancement reflects the gravitational influence of PBHs on small-scale density perturbations, with stronger effects for larger values of f_PBH.

Fig. 1. Impact of PBHs on the power spectrum and HMF (assuming MPBH = 100 M⊙). Top panel: total power spectrum Ptot(k) as a function of wavenumber k for various fractions of DM in PBHs (fPBH = 0.0001, 0.001, 0.01, 0.1, and 1), compared to the standard ΛCDM case (fPBH = 0). Bottom panel: HMF at z = 10, showing the number density of halos as a function of halo mass M, also assuming MPBH = 100 M⊙.

The bottom panel of Fig. 1 presents the resulting HMF at z = 10 for the same set of f_PBH and assuming M_PBH = 100 M_⊙. For this PBH mass, the differences relative to the standard ΛCDM case are most prominent for low-mass halos (M ≲ 10⁹ M_⊙), where an increase in their abundance is observed with larger f_PBH, driven by the Poisson effect. In contrast, for halos with M ≳ 10⁹ M_⊙, the HMF for M_PBH = 100 M_⊙ converges to the ΛCDM prediction, indicating negligible influence from PBHs of this mass range on larger halos. As discussed in Zhang et al. (2024b), the analytical results presented here differ from their simulation findings, which exhibit a significant enhancement in the abundance of massive halos and a suppression of low-mass halos. The latter is driven by nonlinear dynamics, such as the engulfment of smaller halos by PBH-hosting massive halos. The differences arise because the analytical framework used in this study accounts only for linear Poisson-induced modifications to the HMF and does not incorporate nonlinear processes, including the seed effect and the redistribution of matter around PBHs.

In cases of extreme f_PBH values, such as M_PBH = 10 M_⊙ and 33 M_⊙ with f_PBH = 1, and 100 M_⊙ with f_PBH = 0.1 and f_PBH = 1, the integration of the modified HMF over all halos resulted in a total DM density in halos (${\rho }_{\mathrm{DM}}^{\mathrm{halos}}$) exceeding the total cosmological DM density (${\rho }_{\mathrm{DM}}^{\mathrm{total}}$). This outcome arises because the HMF used here incorporates modifications due to PBHs, which significantly alter DM distribution compared to the standard ΛCDM HMF (Zhang et al. 2024b). These modifications account for the enhanced abundance of low-mass halos seeded by PBHs, especially in scenarios where PBHs dominate the DM composition (f_PBH ∼ 1).

To address this inconsistency, a correction factor was applied to scale the HMF, ensuring that the total DM density in halos matches the cosmological value (${\rho }_{\mathrm{DM}}^{\mathrm{halos}}={\rho }_{\mathrm{DM}}^{\mathrm{total}}$). This adjustment effectively assigns all DM to halos, setting the DM density in the IGM to zero (${\rho }_{\mathrm{DM}}^{\mathrm{IGM}}=0$) for these specific configurations, i.e. in a Universe where massive PBHs dominate the DM budget, and halos encapsulate the entirety of the available DM.

The correction factors applied to the HMF are close to unity across all configurations and were only necessary for f_PBH = 1 and for M_PBH = 100 M_⊙, also for f_PBH = 0.1. In most cases, the factors deviate by less than 10 per cent, with the most significant deviation observed in the extreme scenario of M_PBH = 100 M_⊙ and f_PBH = 1, where the factor remains below 30 per cent. As discussed in Section 3.2, we tested the impact of using the unmodified HMF, defined here as the standard ΛCDM HMF without accounting for the effects of PBHs, instead of the modified one. These tests revealed that the derived constraints on f_PBH remain largely unaffected, as the corrections are applied primarily in extreme configurations with large fractions of DM composed of PBHs. For lower f_PBH values, where the constraints are typically established, the differences between the modified and unmodified HMFs are negligible.

It is worth noting that the dominant contribution of PBHs to observable cosmic backgrounds arises from accretion within halos, where the densities are sufficient to sustain significant accretion rates. In contrast, the contribution from PBHs in the IGM is negligible. Consequently, ρ_DM, IGM = 0 in these extreme cases does not impact the main results, as the dominant contributions to cosmic backgrounds are governed by the accretion processes occurring within halos.

2.2. Accretion and emission from PBHs

We considered the same PBH accretion model as in CV2024, outlined below.

2.2.1. Bondi–Hoyle–Lyttleton model

We assumed the Bondi–Hoyle–Lyttleton accretion model characterises the accretion processes of PBHs in both the IGM and within halos. The accretion rate is given by

$$\begin{array}{c}\hfill \dot{M}=4\pi r_{\mathrm{B}}^2\stackrel{\sim }{v}\rho =\frac{4\pi G^2{M}^{2}n\mu m_{\mathrm{p}}}{{\stackrel{\sim }{v}}^3},\end{array}$$(14)

where $r_{\mathrm{B}}=GM/{\stackrel{\sim }{v}}^2$ is the Bondi radius, n is the gas number density, μ is the mean molecular weight, for which we adopted μ = 1.22, corresponding to that of a primordial neutral gas. m_p = 1.673 × 10⁻²⁴ g is the proton mass, and $\stackrel{\sim }{v}\equiv {(v^2+c_{\mathrm{s}}^2)}^{1/2}$ represents the characteristic velocity. Here, v is the relative velocity of the PBH with respect to the gas, and c_s is the sound speed of the gas.

It is worth noting that, following the prescription of Takhistov et al. (2022), our formulation does not include a Bondi-Eigenvalue factor λ, which is sometimes introduced in the literature as an ad hoc correction to account for uncertainties in the accretion process. In this approach, the accretion and emission regimes are determined directly from the local physical conditions of the PBH and the surrounding medium, without the need for additional free parameters.

2.2.2. Condition for accretion disc formation

To determine if an accretion disc forms around a PBH of mass M, we examine whether the outer edge of the PBH accretion disc, r_out, lies within the radius of the innermost stable circular orbit (ISCO). Efficient disc formation is unlikely if the outer edge is within the ISCO radius. The radius of the outer disc edge is approximately given by Agol & Kamionkowski (2002)

$$\begin{array}{c}\hfill r_{\mathrm{out}}\simeq 5.4\times {10}^9{r}_{\mathrm{s}}{\left(\frac{M}{100{\mathrm{M}}_{\odot }}\right)}^{\frac{2}{3}}{\left(\frac{\stackrel{\sim }{v}}{10\mathrm{km}/\mathrm{s}}\right)}^{-\frac{10}{3}},\end{array}$$(15)

where r_s = 2GM/c² is the Schwarzschild radius of the BH, and c = 2.998 × 10¹⁰ cm/s is the speed of light.

For a non-rotating, spherically symmetric PBH (Schwarzschild PBH), the ISCO radius is defined as

$$\begin{array}{c}\hfill r_{\mathrm{ISCO}}=3r_{\mathrm{s}}.\end{array}$$(16)

This radius represents the minimum distance at which an object can maintain a stable circular orbit around the PBH.

2.2.3. Accretion regimes

Following Takhistov et al. (2022), we consider different accretion regimes characterised by the dimensionless accretion rate ṁ. This rate is defined as

$$\begin{array}{c}\hfill \dot{m}=\frac{\dot{M}}{{\dot{M}}_{\mathrm{Edd}}}=2.64\times {10}^{-7}\left(\frac{M}{{\mathrm{M}}_{\odot }}\right)\left(\frac{n}{1{\mathrm{cm}}^{-3}}\right)\left(\frac{\mu }{1}\right){\left(\frac{\stackrel{\sim }{v}}{10\mathrm{km}/\mathrm{s}}\right)}^{-3},\end{array}$$(17)

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

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{Edd}}=\frac{L_{\mathrm{Edd}}}{{ϵ}_0{c}^2}=6.7\times {10}^{-16}\left(\frac{M}{{\mathrm{M}}_{\odot }}\right){\mathrm{M}}_{\odot }/\mathrm{s},\end{array}$$(18)

with ϵ₀ = 0.1 being the radiative efficiency characteristic of thin discs. The radiative efficiency ϵ₀ can vary from 0.057 for a non-rotating Schwarzschild BH to 0.42 for an extremal Kerr BH (see e.g. Kato et al. 2008). As in Takhistov et al. (2022), we assume ϵ₀ = 0.1 as a typical value for thin discs.

In this study, we include both thin discs and ADAFs. We consider the formation of a thin disc when ṁ ≥ 0.07 α, where α is the viscosity parameter. For a more robust analysis, we further divide ADAFs into three distinct sub-regimes based on the accretion rate, following Takhistov et al. (2022). Specifically, we consider LHAF for accretion rates $0.1{\alpha }^2<\dot{m}<0.07\alpha$, standard ADAF in the range ${10}^{-3}{\alpha }^2<\dot{m}<0.1{\alpha }^2$, and eADAF for ṁ < 10⁻³α². Below, we provide a description of each regime, utilising the model detailed by Takhistov et al. (2022), which has also been employed in CV2024.

The thin disc formed around PBHs is optically thick, allowing it to efficiently radiate blackbody emission (Shakura & Sunyaev 1973), which facilitates a comprehensive analytical model. The temperature profile of the disc, beyond the inner region, is given by (Pringle 1981)

$$\begin{array}{c}\hfill T(r)=T_{\mathrm{i}}{\left(\frac{r_{\mathrm{i}}}{r}\right)}^{3/4}{[1-{\left(\frac{r_{\mathrm{i}}}{r}\right)}^{\frac{1}{2}}]}^{\frac{1}{4}},\end{array}$$(19)

where r_i represents the inner radius of the disc, assumed to be the ISCO radius, and

$$\begin{array}{c}\hfill T_{\mathrm{i}}=\left(\frac{3GM\dot{M}}{8\pi r_{\mathrm{i}}^3{\sigma }_{\mathrm{B}}}\right)=53.3\mathrm{eV}{\left(\frac{n}{1{\mathrm{cm}}^{-3}}\right)}^{1/4}{\left(\frac{\stackrel{\sim }{v}}{10\mathrm{km}/\mathrm{s}}\right)}^{-3/4},\end{array}$$(20)

where σ_B = 5.67 × 10⁻⁵ g/s³/K⁴ is the Stefan-Boltzmann constant.

Beyond the inner edge, the disc reaches a maximum temperature T_max = 0.488T_i at a radius r = 1.36r_i.

The thin disc spectrum is a combination of blackbody emissions from different radii. Using the scaling relations from Pringle (1981) and ensuring continuity, the resulting spectrum can be approximated as

$$\begin{array}{cc}& \nu <T_{\mathrm{o}}:L_{\nu }=c_{\alpha }{\left(\frac{T_{max}}{T_{\mathrm{o}}}\right)}^{5/3}{\left(\frac{\nu }{T_{max}}\right)}^2,\hfill \\ \hfill & T_{\mathrm{o}}<\nu <T_{max}:L_{\nu }=c_{\alpha }{\left(\frac{\nu }{T_{max}}\right)}^{1/3},\hfill \\ \hfill & T_{max}<\nu :L_{\nu }=c_{\alpha }{\left(\frac{\nu }{T_{max}}\right)}^2{e}^{1-\nu /T_{max}},\hfill \end{array}$$(21)

where ν is the photon energy, T_o = T(r_out) is the temperature at the outer edge of the disc, and

$$\begin{array}{cc}\hfill c_{\alpha }& =1.27\times {10}^{29}\mathrm{erg}\mathrm{eV}{}^{-1}{\mathrm{s}}^{-1}{\left(\frac{M}{{\mathrm{M}}_{\odot }}\right)}^2\hfill \\ \hfill & {\left(\frac{n}{1{\mathrm{cm}}^{-3}}\right)}^{3/4}{\left(\frac{\stackrel{\sim }{v}}{10\mathrm{km}/\mathrm{s}}\right)}^{-9/4}.\hfill \end{array}$$(22)

c_α is normalised such that the emission achieves the maximum possible efficiency for a Schwarzschild BH (${\displaystyle \int L_{\nu }}\mathrm{d}\nu =0.057\dot{M}{c}^2$).

In contrast to the thin disc, when an ADAF forms, the heat generated by viscosity is not efficiently radiated away, and a significant amount of energy is advected into the BH event horizon along with the gas. This results in a complex emission spectrum with contributions from electron cooling, synchrotron radiation, inverse-Compton (IC) scattering, and bremsstrahlung processes.

To characterise the flow, we used the following parameters: the fraction of viscously dissipated energy that heats electrons directly, δ = 0.3; the ratio of gas pressure to total pressure, β = 10/11; the minimum flow radius, equal to the ISCO radius, r_min = 3r_s; and the viscosity parameter, α = 0.1. These values are consistent with recent simulations and observations (Yuan & Narayan 2014), and are also used in the model by Takhistov et al. (2022).

The synchrotron emission is self-absorbed and peaks at a photon energy of

$$\begin{array}{cc}\hfill {\nu }_{\mathrm{p}}& =1.83\times {10}^{-2}\mathrm{eV}{\left(\frac{\alpha }{0.1}\right)}^{-\frac{1}{2}}{\left(\frac{1-\beta }{1/11}\right)}^{\frac{1}{2}}\hfill \\ \hfill & {\theta }_{\mathrm{e}}^2{\left(\frac{r_{min}}{3r_{\mathrm{s}}}\right)}^{-\frac{5}{4}}{\left(\frac{M}{{\mathrm{M}}_{\odot }}\right)}^{-\frac{1}{2}}{\left(\frac{\dot{m}}{{10}^{-8}}\right)}^{\frac{3}{4}},\hfill \end{array}$$(23)

where θ_e represents the temperature in units of the electron mass m_e

$$\begin{array}{c}\hfill {\theta }_{\mathrm{e}}=\frac{k{T}_{\mathrm{e}}}{m_{\mathrm{e}}{c}^2}=\frac{T_{\mathrm{e}}}{5.93\times {10}^9\mathrm{K}}.\end{array}$$(24)

where k_B = 1.3807 × 10⁻¹⁶ cm² g s⁻² K⁻¹ is the Boltzmann constant and T_e the electron temperature. The peak luminosity is given by

$$\begin{array}{c}\hfill L_{{\nu }_{\mathrm{p}}}=5.06\times {10}^{38}\mathrm{erg}/\mathrm{s}/\mathrm{eV}{\alpha }^{-1}(1-\beta )\left(\frac{M}{{\mathrm{M}}_{\odot }}\right){\dot{m}}^{3/2}{\theta }_{\mathrm{e}}^5{\left(\frac{r_{\min }}{r_{\mathrm{s}}}\right)}^{-1/2}.\end{array}$$(25)

In this simple description, the synchrotron spectrum is assumed to terminate at ν_p, resulting in L_ν, syn = 0 for ν > ν_p (Mahadevan 1997).

The synchrotron photons undergo IC scattering with the surrounding electron plasma. The resulting IC spectrum is described by

$$\begin{array}{c}\hfill L_{\nu ,\mathrm{IC}}=L_{{\nu }_{\mathrm{p}}}{\left(\frac{\nu }{{\nu }_{\mathrm{p}}}\right)}^{-{\alpha }_{\mathrm{c}}},\end{array}$$(26)

where

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{c}}=-\frac{ln{\tau }_{\mathrm{e}\mathrm{s}}}{lnA},\end{array}$$(27)

with ${\tau }_{\mathrm{e}\mathrm{s}}=12.4\dot{m}{\alpha }^{-1}{\left(\frac{r_{\min }}{r_{\mathrm{s}}}\right)}^{-1/2}$ as the electron scattering optical depth, and A = 1 + 4θ_e + 16θ_e². This expression for L_ν, IC is valid only in the frequency range ν_p ≤ ν ≲ 3k_BT_e, which corresponds to the energy range where thermal IC scattering significantly modifies the photon spectrum. Beyond this range, photons either lack sufficient energy for further upscattering or are dominated by other emission processes.

Finally, the bremsstrahlung emission from the thermal spectrum is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill L_{\nu ,\text{brems}}& =1.83\times {10}^{17}{\left(\frac{\alpha }{0.1}\right)}^{-2}{\left(\frac{c_1}{0.5}\right)}^{-2}ln\left(\frac{r_{max}}{r_{min}}\right)F\left({\theta }_{\mathrm{e}}\right)\hfill \\ \hfill & {\left(\frac{T_{\mathrm{e}}}{5\times {10}^9\mathrm{K}}\right)}^{-1}{e}^{-(h\nu /k{T}_{\mathrm{e}})}\left(\frac{M}{{\mathrm{M}}_{\odot }}\right){\dot{m}}^2{\mathrm{ergs}\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1},\hfill \end{array}\end{array}$$(28)

where c₁ = 0.5 and F(θ_e) is given by Mahadevan (1997)

$$\begin{array}{c}\hfill F({\theta }_{\mathrm{e}})=\{\begin{array}{cc}\begin{array}{cc}& 4{\left(\frac{2{\theta }_{\mathrm{e}}}{{\pi }^3}\right)}^{1/2}(1+1.78{\theta }_{\mathrm{e}}^{1.34})\hfill \\ \hfill & +1.73{\theta }_{\mathrm{e}}^{3/2}(1+1.1{\theta }_{\mathrm{e}}+{\theta }_{\mathrm{e}}^2-1.25{\theta }_{\mathrm{e}}^{5/2}),\hfill \end{array}\hfill & {\theta }_{\mathrm{e}}\le 1,\hfill \\ \begin{array}{cc}& \left(\frac{9{\theta }_{\mathrm{e}}}{2\pi }\right)(ln[1.12{\theta }_{\mathrm{e}}+0.48]+1.5)\hfill \\ \hfill & +2.30{\theta }_{\mathrm{e}}(ln\left[1.12{\theta }_{\mathrm{e}}\right]+1.28),\hfill \end{array}\hfill & {\theta }_{\mathrm{e}}\ge 1.\hfill \end{array}\end{array}$$(29)

Regarding frequency ranges, synchrotron emission dominates the radio to sub-millimetre spectrum due to its origin in the thermal electron population of the ADAF. The IC scattering process primarily contributes to the X-ray region, with photon energies typically extending up to ∼3k_BT_e. This upper limit corresponds to the maximum energy gain for photons undergoing thermal IC scattering before other emission mechanisms become dominant. Finally, bremsstrahlung emission is significant in the hard X-ray region, where it provides a relatively constant contribution across a wide frequency range before an exponential cutoff dictated by the high-energy tail of the electron distribution.

Each ADAF subregime has different temperature dependences because θ_e is determined by balancing heating and radiative processes (Mahadevan 1997). Direct viscous electron heating dominates at low accretion rates, while ion-electron collisional heating becomes significant at higher rates. The methodology for calculating θ_e in each subregime can be found in Section 3.1.3 of Takhistov et al. (2022).

We acknowledge that the emission model adopted in this work, following Takhistov et al. (2022), involves several simplifying assumptions that may limit its physical completeness. These include the use of fixed microphysical parameters–such as the viscosity parameter α, the plasma β (which reflects the relative importance of gas pressure over magnetic pressure), and the electron heating fraction δ–which are kept constant across all accretion regimes where they apply. The model also considers only thermal electron populations and does not include non-thermal particles, which could affect the emitted spectrum in some ADAF scenarios. Furthermore, the emission is computed assuming spherical symmetry in the geometry of the accretion flow, neglecting potential effects of angular momentum and anisotropies in the gas surrounding PBHs. While the slim disc regime, relevant at high Eddington ratios, is not included in this model, we explicitly verified that its omission does not affect our results: the super-Eddington regime is only marginally reached in the densest central regions of the most massive halos, which contribute negligibly to the total background. These assumptions are standard in the ADAF literature and enable a tractable and physically motivated framework, but they may introduce systematic uncertainties that cannot be robustly quantified within the scope of this work. We highlight these caveats as limitations of the current implementation and note that they warrant further investigation.

2.3. PBH accretion in halos

As discussed earlier, the accretion rate of PBHs depends not only on the PBH mass but also on the characteristic velocity of the surrounding medium and the gas density. Within halos, we assume that the virial velocity of the halos determines the characteristic velocity

$$\begin{array}{c}\hfill \stackrel{\sim }{v}\approx \sqrt{\frac{G{M}_{\mathrm{h}}}{r_{\mathrm{vir}}}}\sim 5.4{\mathrm{km}\mathrm{s}}^{-1}{\left(\frac{M_{\mathrm{h}}}{{10}^6{\mathrm{M}}_{\odot }}\right)}^{-1/2}{\left(\frac{1+z}{21}\right)}^{-1/2}.\end{array}$$(30)

Regarding the gas density, we consider three distinct halo density profiles to assess their impact on the calculated radiation:

• Simple isothermal density profile. The first density profile considered is a simple isothermal profile, where the density distribution follows an inverse-square law given by

  $$\begin{array}{c}\hfill n(r)=\frac{n_0}{3}{\left(\frac{r_{\mathrm{vir}}}{r}\right)}^2,\end{array}$$(31)

  where n₀ is the average gas density, defined as

  $$\begin{array}{c}\hfill n_0=200{\rho }_{\mathrm{m}}(a)\frac{{\mathrm{\Omega }}_{\mathrm{b}}}{{\mathrm{\Omega }}_{\mathrm{m}}}\frac{1}{\mu m_{\mathrm{p}}},\end{array}$$(32)

  with ρ_m(a) the matter density of the Universe as a function of the scale factor a, given by

  $$\begin{array}{c}\hfill {\rho }_{\mathrm{m}}(z)=\frac{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}{8\pi G}{(1+z)}^3,\end{array}$$(33)

  This profile assumes a spherically symmetric distribution where the density falls off as r⁻², providing a simple and widely used model in astrophysical studies. Simulations indicate that the gas distribution in high-z atomic-cooling halos approximately follows ρ ∝ r⁻² at r ≳ 0.003 pc (Safarzadeh & Haiman 2020). Similar density profiles are observed in the simulations by Liu et al. (2022) for molecular-cooling minihalos. Consequently, Liu et al. (2022) employs this profile for their calculation of the radiation background produced by PBHs.

• Makino+1998 density profile. The second density profile is the one derived in Makino et al. (1998). The gas density profile in hydrostatic equilibrium within the potential of a DM halo, assuming an isothermal distribution, can be approximated by the β-model

  $$\begin{array}{c}\hfill {\rho }_{\mathrm{gas}}={\rho }_{0}exp\{-\frac{\mu m_{\mathrm{p}}}{2k_{\mathrm{B}}{T}_{\mathrm{vir}}}[v_{\mathrm{e}}^2(0)-v_{\mathrm{e}}^2(r)]\},\end{array}$$(34)

  where ρ₀ is the central gas density

  $$\begin{array}{c}\hfill {\rho }_0(z)=\frac{({\mathrm{\Delta }}_{\mathrm{c}}/3)c^3{e}^A}{{\displaystyle {\int }_0^c}{(1+t)}^{A/t}{t}^2\mathrm{d}t}\left(\frac{{\mathrm{\Omega }}_{\mathrm{b}}}{{\mathrm{\Omega }}_{\mathrm{m}}}\right){\mathrm{\Omega }}_{\mathrm{m}}^z{\rho }_{\mathrm{c}}(z),\end{array}$$(35)

  where A = 2c/F(c), with F(c) defined in Eq. (7), and ρ_c(z) = 3H(z)²/(8πG) represents the critical density of the Universe at a redshift z. The viral temperature of the halo, T_vir, is given by

  $$\begin{array}{c}\hfill T_{\mathrm{vir}}=1.98\times {10}^4\left(\frac{\mu }{0.6}\right){\left(\frac{M_{\mathrm{vir}}}{{10}^8{h}^{-1}{\mathrm{M}}_{\odot }}\right)}^{2/3}{\left[\frac{{\mathrm{\Omega }}_{\mathrm{m}}{\mathrm{\Delta }}_{\mathrm{c}}}{{\mathrm{\Omega }}_{\mathrm{m}}^{z}18{\pi }^2}\right]}^{1/3}(1+z)\mathrm{K},\end{array}$$(36)

  and v_e is the escape velocity

  $$\begin{array}{c}\hfill v_{\mathrm{e}}^2(r)=2{\displaystyle {\int }_r^{\infty }}\frac{GM(r\prime )}{r{\prime }^2}\mathrm{d}r\prime =2v_{\mathrm{c},\mathrm{vir}}^2\frac{F(cx)+\frac{\mathit{cx}}{1+cx}}{xF(c)},\end{array}$$(37)

  where $v_{c,\mathrm{vir}}=\sqrt{G{M}_{\mathrm{vir}}/r_{\mathrm{vir}}}$ is the virial velocity. These equations collectively describe how the gas density profile is determined in the context of a universal DM halo, accounting for the effects of hydrostatic equilibrium and the influence of the DM distribution. Ziparo et al. (2022) use this profile to calculate the radiation backgrounds produced by PBHs.

• Rescaled halo from Liu+2022. The third density profile is derived by rescaling a halo from the study by Liu et al. (2022). In this model, the density profile at the onset of star formation at z ≈ 30 within a halo of 2 × 10⁵ M_⊙ is utilised. The profile is rescaled to match the virial mass of the halos in our study by multiplying the reference profile by the ratio of the enclosed mass within our halos at their respective virial radii to the enclosed mass of the reference halo at its virial radius. This method was also employed in CV2024 to address resolution issues and to better account for the density distribution within halos at z ≈ 30. This profile is only employed up to z = 20, beyond which its applicability diminishes. At lower redshifts, hierarchical mergers, baryonic feedback processes (e.g. supernova-driven outflows and radiation pressure), and longer dynamical timescales significantly alter halo density structures. These effects deviate from the assumptions underlying the rescaled profile derived from a halo at z ∼ 30. Limiting its use to z = 20 ensures the validity of the model within the appropriate dynamical and physical regimes.

As illustrated in Fig. 2, the hydrogen number density profiles exhibit distinct behaviours depending on the chosen model and halo mass, with significant variations observed across different redshifts. These differences increase with halo mass as the density profiles deviate more substantially from one another at higher masses. However, we calculated that the contribution of halos with M_h > 10⁶ M_⊙ to the total emission is negligible, as halos with M_h < 10⁶ M_⊙ dominate the total emission, contributing the overwhelming majority of the radiation.

Fig. 2. Comparison of hydrogen number density profiles nH(r) within DM halos of different masses (Mh = 104 M⊙, 106 M⊙, and 108 M⊙) at redshifts z = 6 (solid lines), z = 10 (dashed-dotted lines), and z = 40 (dotted lines). The three density profiles examined include a simple isothermal profile (magenta), the profile derived by Makino et al. (1998) (blue), and the rescaled profile from Liu et al. (2022) (green). The vertical grey lines denote the virial radii of the corresponding halos at each redshift.

2.4. PBH accretion in the IGM

Following Ricotti et al. (2008), we assume that PBHs in the IGM are surrounded by a uniformly distributed gas with a density given by

$$\begin{array}{c}\hfill {\rho }_{\mathrm{IGM}}=250\mu m_{\mathrm{p}}{\left(\frac{1+z}{1000}\right)}^3{\mathrm{g}\mathrm{cm}}^{-3}.\end{array}$$(38)

This relation provides the mean gas density in the IGM under the assumption of a homogeneous and isotropic Universe. As discussed by Ricotti et al. (2008), the scaling ρ_IGM ∝ (1 + z)³ reflects the cosmic expansion, where the gas density evolves inversely with the comoving volume. The prefactor 250μ corresponds to the baryon density at z = 1000, with μ accounting for the primordial composition of the gas, dominated by hydrogen and helium. This density serves as a crucial input for modelling gas dynamics and accretion processes around PBHs.

While we follow the assumption of a uniform gas distribution in the IGM, as adopted in Ziparo et al. (2022), it is important to acknowledge that the presence of DM halos seeded by PBHs could concentrate the surrounding gas. This effect may enhance accretion rates and luminosity, potentially influencing feedback processes and the thermal evolution of the IGM. However, our constraints indicate that the fraction of DM composed of PBHs is quite small in order to remain consistent with the observational and theoretical limits explored in this work, suggesting that this effect is probably not important.

The sound speed of the gas, c_s, is determined using the parametric relation provided by Luca et al. (2020), which incorporates the evolution of the cmB temperature with redshift. This expression is given by

$$\begin{array}{c}\hfill c_{\mathrm{s}}\simeq 5.70{\left(\frac{1+z}{1000}\right)}^{1/2}{[{\left(\frac{1+z_{\mathrm{dec}}}{1+z}\right)}^{\beta }+1]}^{-1/(2\beta )}{\mathrm{km}\mathrm{s}}^{-1},\end{array}$$(39)

where β = 1.72 is a fitting parameter, and z_dec ≃ 130 corresponds to the redshift at which baryons decouple from the radiation field. This formulation captures the transition from tightly coupled baryons and radiation at early times to a thermally evolving IGM. The inclusion of this relation ensures an accurate representation of the baryonic sound speed across cosmic epochs, which is critical for assessing gas dynamics and the efficiency of accretion onto PBHs.

To model the relative velocity between DM and baryons, we adopted the relation proposed by Ali-Haïmoud & Kamionkowski (2017)

$$\begin{array}{c}\hfill v_{\mathrm{rel}}(z)\simeq 30min[1,\left(\frac{1+z}{1000}\right)]{\mathrm{km}\mathrm{s}}^{-1}.\end{array}$$(40)

This residual streaming velocity originates after the recombination epoch, when baryons kinetically decouple from photons at z ∼ 1000. The root-mean-square (RMS) value of this velocity is approximately 30 km/s at recombination, as calculated by Ali-Haïmoud & Kamionkowski (2017). The velocity decreases as (1 + z) at earlier times, reflecting the expansion of the Universe, and remains constant for z ≲ 1000. This streaming motion is crucial for understanding the thermal and dynamical evolution of the IGM, as it affects small-scale structure formation and the interactions between baryons and DM (see also Kashlinsky 2021; Atrio-Barandela 2022).

2.5. Radiation backgrounds

The radiation background produced by PBHs accreting within halos and the IGM is analysed across three frequency regimes: X-ray, LW, and radio. This section outlines the equations used to characterise the background intensity in each regime, while a detailed derivation is provided in Appendix A.

We acknowledge that the Cosmic Infrared Background (CIB) has also been studied as a potential probe of PBHs. Some works have examined whether accreting PBHs could account for the observed CIB–CXB cross-correlation (e.g. Kashlinsky 2016; Cappelluti et al. 2017, 2022). More recent analyses incorporating detailed accretion physics and thermal feedback suggest that PBHs are subdominant contributors. For example, Hasinger (2020), assuming a broad PBH mass function extending from 10⁻⁸ to 10¹⁰ M_⊙, Bondi-like accretion from homogeneous background gas, and a PBH abundance equal to the total DM, estimated a contribution of only ∼0.5% to the CIB fluctuations. Manzoni et al. (2023), considering both halo and background-gas accretion with a self-consistent treatment of thermal evolution, evaluated PBH contributions to the CIB for both monochromatic and lognormal mass functions in the range 1–10³ M_⊙. In the limiting case where PBHs constitute all of the DM, they found a contribution below 1% to the observed NIRB fluctuations. This value drops to below 0.1% when current observational constraints on PBH abundance are taken into account, with a maximal contribution at M_PBH ∼ 50 M_⊙ within the allowed parameter space. Given this, and our focus on radiation backgrounds that place more direct constraints on the PBH contribution to DM, we do not consider the CIB further in this work.

The comoving emissivity is calculated differently for halos and the IGM. In the case of halos, integration over the HMF is required, whereas for the IGM, the emissivity is determined directly from the PBH number density.

For halos, the comoving emissivity within a specific frequency band ν ∈ [ν₁, ν₂] is given by

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{cm},[{\nu }_1,{\nu }_2]}(z)={\displaystyle \int L_{[{\nu }_1,{\nu }_2]}}(M_{\mathrm{h}},z)\frac{\mathrm{d}{n}_{\mathrm{h}}}{\mathrm{d}{M}_{\mathrm{h}}}\mathrm{d}{M}_{\mathrm{h}},\end{array}$$(41)

where L_[ν1, ν2](M_h, z) represents the total luminosity of a halo with mass M_h at redshift z. This luminosity is computed using the accretion model described in Section 2.2.3, incorporating the gas density profiles and characteristic velocities discussed in Section 2.3. The term dn_h/dM_h corresponds to the HMF, which provides the comoving number density of halos per unit halo mass, further detailed in Section 2.1.2.

For PBHs in the IGM, the comoving emissivity is expressed as

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{cm},[{\nu }_1,{\nu }_2]}(z)=L_{[{\nu }_1,{\nu }_2]}(z)n_{\mathrm{IGM}}(z),\end{array}$$(42)

where L_[ν1, ν2](z) denotes the total luminosity per PBH within the band [ν₁, ν₂] in the IGM, and n_IGM(z) represents the comoving PBH number density, as defined in Section 2.1.1. As with halos, the luminosity is computed using the formalism in Section 2.2.3, with gas density and velocity profiles specific to IGM conditions, discussed in Section 2.4.

The background intensity in the X-ray band, observed at z = 0, is obtained by discretising the redshift integration to account for the contributions ΔJ_X, i of individual redshift bins

$$\begin{array}{c}\hfill J_{\mathrm{X}}={\displaystyle \sum _i}\mathrm{\Delta }{J}_{\mathrm{X},i}={\displaystyle \sum _i}\frac{1}{4\pi }\frac{{ϵ}_{\mathrm{cm},\mathrm{X}}(z_i)}{{(1+z_i)}^2}\mathrm{\Delta }{r}_i,\end{array}$$(43)

where ϵ_cm, X(z_i) is the comoving X-ray emissivity at redshift z_i, and Δr_i represents the comoving distance step.

For the LW band, the integrated physical intensity observed at z is first calculated as

$$\begin{array}{c}\hfill J_{\mathrm{LW},\mathrm{int}}(z)=\frac{{(1+z)}^4}{4\pi }{\displaystyle {\int }_z^{z_{max}(z)}}\frac{{ϵ}_{\mathrm{cm},\mathrm{LW}}(z,z\prime )}{{(1+z\prime )}^3}c{\frac{\mathrm{d}{t}_{\mathrm{U}}}{\mathrm{d}a\prime }|}_{z\prime }\mathrm{d}z\prime ,\end{array}$$(44)

where t_U is the cosmic age and z_max(z) = (13.6/11.2)(1 + z)−1, which corresponds to the maximum redshift from which photons originally emitted at 13.6 eV are redshifted to 11.2 eV at z, marking the relevant range for the LW band. The corresponding dimensionless specific intensity is given by

$$\begin{array}{c}\hfill J_{\mathrm{LW}}(z)=\frac{{10}^{21}{J}_{\mathrm{LW},\mathrm{int}}(z)/\mathrm{\Delta }{\nu }_{\mathrm{LW}}}{\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}{\mathrm{sr}}^{-1}{\mathrm{Hz}}^{-1}},\end{array}$$(45)

where Δν_LW = (13.6 − 11.2) eV/h_P, with h_P = 4.1357 × 10⁻¹⁵ eVs being the Planck constant.

For the radio band, the physical background intensity at an observed frequency ν is computed as

$$\begin{array}{c}\hfill J_{\nu }=\frac{1}{4\pi }{\displaystyle {\int }_z^{\infty }}\frac{{(1+z)}^3}{{(1+z\prime )}^2}{\dot{\rho }}_{{\nu }^{\prime }}(z^{\prime })c{\frac{\mathrm{d}{t}_{\mathrm{U}}}{\mathrm{d}{a}^{\prime }}|}_{z^{\prime }}\mathrm{d}z\prime .\end{array}$$(46)

where ${\dot{\rho }}_{\nu \prime }(z\prime )\equiv {ϵ}_{\nu \prime }(z\prime )/{(1+z\prime )}^3$ represents the comoving emissivity at source-frame frequency ν′=ν[(1 + z′)/(1 + z)]. The corresponding brightness temperature, which is a crucial observable for 21 cm cosmology, is calculated as

$$\begin{array}{c}\hfill T_{\mathrm{b}}(\nu ,z)=\frac{J_{\nu }(z)c^2}{2k_{\mathrm{B}}{\nu }^2}.\end{array}$$(47)

3. Results and discussion

In this section, we analyse the contribution of PBHs to the CXB, LWB, and CRB, exploring their potential role as significant sources of these cosmic backgrounds and as candidates for DM. Using theoretical predictions for different PBH masses and fractions, we compare their impact with current observational constraints. The analysis considers contributions from both halos and the IGM, highlighting the interplay between these components in shaping the radiation backgrounds. Finally, we discuss the implications of these results for the thermal and ionisation history of the Universe, the suppression of star formation, and the constraints they impose on the PBH parameter space.

3.1. Background intensity

3.1.1. X-ray background

Approximately 70 per cent of the CXB is attributed to resolved sources, predominantly AGNs, across both soft and hard X-ray bands (Cappelluti et al. 2017, 2022). The unresolved fraction of the CXB is likely due to a yet undetected population of BHs, such as those accreting in heavily obscured environments (Gilli et al. 2007; Treister et al. 2009), as well as extended X-ray emissions from galaxy clusters (Gilli et al. 1999). In addition, there could be a contribution from high-z AGNs, linked to the recent JWST discovery of such abundant sources, but with significant uncertainties linked to the physics of their X-ray emission (e.g. Jeon et al. 2022; Maiolino et al. 2025).

To investigate the potential contribution of PBHs to this unresolved background, we quantify their impact on the X-ray background in the soft (0.5–2 keV) and hard (2–10 keV) X-ray bands. Our analysis specifically aims to determine whether PBHs can account for a significant portion of the unresolved CXB, thereby offering insights into the nature of DM.

The results in Fig. 3 illustrate the cumulative contribution of accreting PBHs to the present-day soft X-ray background intensity, J_{0.5 − 2 keV}, for various PBH masses and fractions f_PBH. The shaded grey region indicates the observed excess X-ray background reported by Cappelluti et al. (2017), which serves as a constraint on unresolved contributions to the CXB.

Fig. 3. Contribution of accreting high-z BH sources to the present-day (z = 0) soft X-ray background (0.5–2 keV). The top row shows the cumulative evolution of the integrated background intensity, JX, as a function of redshift z, for PBH masses of MPBH = 1 M⊙, 10 M⊙, 33 M⊙, and 100 M⊙ in each column. Different values of fPBH are represented by distinct colours, as indicated in the legend, ranging from 10−4 to 1. The grey shaded region corresponds to the observed excess X-ray background reported in Cappelluti et al. (2017). Three halo density profiles are considered: simple isothermal (dashed lines), the profile from Makino et al. (1998) (dashed-dotted lines), and the profile from Liu et al. (2022) (dotted lines). The bottom row illustrates the relative contributions of halos ($J_{\mathrm{X}}^{\mathrm{h}}$) and the IGM ($J_{\mathrm{X}}^{\mathrm{IGM}}$) to the total X-ray background intensity ($J_{\mathrm{X}}^{\mathrm{tot}}$). These are represented as $J_{\mathrm{X}}^{\mathrm{h}}/J_{\mathrm{X}}^{\mathrm{tot}}$ (dashed-dotted lines) and $J_{\mathrm{X}}^{\mathrm{IGM}}/J_{\mathrm{X}}^{\mathrm{tot}}$ (solid lines), respectively. The relative contributions are computed using the same PBH masses and density profile configurations as in the top row.

The analysis shows that the PBH contribution depends significantly on both the mass and the fraction f_PBH. For the simple isothermal and Makino+1998 density profiles, PBHs with M_PBH = 1 M_⊙ are ruled out for f_PBH ≥ 0.01, while for M_PBH = 10 M_⊙, 33 M_⊙, and 100 M_⊙, fractions f_PBH ≥ 0.001 are excluded due to their overproduction of soft X-ray background intensity. The results for the simple isothermal and Makino+1998 profiles are fully consistent across all cases.

The results based on the rescaled halo profile from Liu et al. (2022) are shown only up to z = 20, due to the redshift limitations of the model, as discussed in Section 2.3. Within this range, the constraints are less restrictive compared to the simple isothermal and Makino+1998 profiles. For M_PBH = 1 M_⊙, fractions f_PBH = 1 are excluded, while for M_PBH = 10 M_⊙, 33 M_⊙, and 100 M_⊙, fractions f_PBH ≥ 0.1, f_PBH ≥ 0.01, and f_PBH ≥ 0.01, respectively, are ruled out.

The differences in constraints between the Liu+2022 halo profile and the other two profiles stem from variations in the predicted gas densities across different halo masses. For halos with M_h < 5 × 10⁵ M_⊙, the Liu+2022 profile predicts lower gas densities compared to the simple isothermal and Makino+1998 profiles, leading to reduced X-ray emissions from these low-mass halos. Since halos in this mass range dominate the total X-ray luminosity due to their abundance, as discussed in Section 2.3, this generally results in lower overall X-ray output for the Liu+2022 profile.

However, the Liu+2022 profile predicts higher gas densities for halos with M_h ≥ 5 × 10⁵ M_⊙, and this difference becomes relevant for certain configurations. Specifically, for higher PBH masses (M_PBH = 33 M_⊙ and M_PBH = 100 M_⊙) and very low fractions (f_PBH = 10⁻⁴), the increased contribution from halos in the intermediate mass range (5 × 10⁵ M_⊙ ≤ M_h ≤ 1 × 10⁶ M_⊙) leads to slightly higher X-ray emissions for the Liu+2022 profile compared to the other two profiles.

This behaviour reflects a balance between the higher gas densities of intermediate-mass halos and the greater abundance of low-mass halos. For small fractions of PBHs (f_PBH = 10⁻⁴), the number of PBHs increases significantly, amplifying the contribution from the intermediate-mass halos where the Liu+2022 profile predicts higher densities. This interplay results in slightly higher emissions for the Liu+2022 profile in these specific cases, despite its generally lower predictions for the total X-ray output. This distinction highlights the sensitivity of PBH constraints to the choice of halo density profile and the complex role of halo mass distribution in shaping X-ray emissions.

The bottom panels of Fig. 3 illustrate the relative contributions of PBHs in halos and the IGM to the total soft X-ray intensity. In general, the contribution from the IGM decreases with decreasing redshift, while the contribution from halos becomes dominant at lower redshifts, as expected. Across the redshift range studied (z = 6 to z = 40), the IGM contribution is subdominant in most cases, with the only exception being for f_PBH = 10⁻⁴ and M_PBH ≥ 10 M_⊙, where the IGM briefly dominates the emission at higher redshifts. Notably, for higher PBH masses, the redshift of the transition from IGM-dominated to halo-dominated emission shifts to lower values, reflecting the delayed onset of efficient accretion within halos for more massive PBHs. Nevertheless, for such low fractions, the total emission from both the IGM and halos is negligible in comparison to observational constraints.

Interestingly, for certain configurations with f_PBH > 10⁻⁴, the IGM contribution increases rather than decreases with decreasing redshift. This behaviour is observed, for example, for M_PBH = 1 M_⊙ with f_PBH ≥ 10⁻¹, M_PBH = 10 M_⊙ with f_PBH ≥ 10⁻², and M_PBH = 33 M_⊙ with f_PBH ≥ 10⁻². However, even in these cases, the IGM contribution remains negligible compared to the dominant emission from halos.

This counterintuitive increase in IGM emission at lower redshifts can be understood by examining the interplay between accretion dynamics and the evolving density profiles of individual halos. As illustrated in Fig. 2, the central densities of halos decrease with decreasing redshift for a given halo mass. This trend reflects the physical evolution of gas within halos, as the interplay between cooling, accretion, and internal feedback processes shapes the density structure over time. In contrast, the IGM remains less affected by such processes, leading to its relative contribution to X-ray emission becoming more noticeable under specific conditions.

Additionally, this behaviour is closely tied to the balance between the number of PBHs and their individual luminosities. For higher f_PBH, the overall contribution of the IGM is naturally enhanced due to the larger population of PBHs accreting in the diffuse medium. At the same time, for M_PBH ≥ 10 M_⊙, the relative importance of IGM accretion increases further as halo accretion becomes less efficient due to the decreasing gas densities at lower redshifts. This balance is particularly relevant at high fractions f_PBH because the combined luminosity from a larger number of accreting PBHs in the IGM can partially counteract the suppressed accretion rates within halos.

However, even in these scenarios, the total X-ray emission from the IGM remains significantly lower than that from halos, as the latter dominates the total intensity across all configurations.

In addition to the soft X-ray background, we also investigated the contribution of PBHs to the hard X-ray background (2–10 keV). The constraints for M_PBH = 33 M_⊙ and 100 M_⊙ are identical to those derived from the soft X-ray background, with f_PBH ≤ 0.001. However, for lower masses, the constraints are less restrictive in the hard X-ray regime. Specifically, M_PBH = 1 M_⊙ is allowed for fractions f_PBH ≤ 0.1, while M_PBH = 10 M_⊙ is consistent with fractions f_PBH ≤ 0.01. These findings suggest that the hard X-ray background provides complementary but generally weaker constraints compared to the soft X-ray background for low-mass PBHs.

We refine our calculations of f_PBH by deriving tighter limits that satisfy the observed unresolved soft X-ray background constraints and hence provide more precise results. Using the Makino+1998 and simple isothermal profiles, which yield consistent results, we find that for M_PBH = 1 M_⊙, the refined constraint is f_PBH ≤ 0.007, while for M_PBH = 10 M_⊙, M_PBH = 33 M_⊙, and M_PBH = 100 M_⊙, the constraints are f_PBH ≤ 0.0008, f_PBH ≤ 0.0006, and f_PBH ≤ 0.0007, respectively. These refined fractions allow PBHs to explain up to 99, 93, 80, and 91 per cent of the observed excess in the soft X-ray background, while contributing approximately 33, 37, 33, and 39 per cent to the observed excess in the hard X-ray background.

Interestingly, the results reveal that for the same fraction of PBHs (f_PBH), smaller masses (M_PBH) can produce slightly higher emissions compared to larger masses. This behaviour arises because smaller PBHs, for a fixed f_PBH, correspond to a larger number density, which enhances the total accretion luminosity. However, this effect is not universally dominant across all PBH fractions but becomes significant only for very small fractions of DM composition (f_PBH ≪ 0.01), where the increased number of PBHs at lower masses compensates for their reduced individual luminosities. For larger fractions (f_PBH ≳ 0.01), the luminosity per PBH begins to dominate, leading to higher contributions from more massive PBHs. This highlights the complex interplay between PBH mass, abundance, and accretion efficiency in determining their contribution to cosmic backgrounds.

3.1.2. Lyman-Werner background

The LWB plays a critical role in the early Universe by dissociating H₂, the primary coolant in low-mass halos, thereby regulating star formation and the collapse of structures. The critical threshold J_LW ∼ 1, where the H₂ destruction rate equals its formation rate, marks a tipping point for suppressing small-scale structure formation. PBHs contributing significantly to the LWB could therefore have profound implications for cosmic evolution by delaying or preventing the formation of stars and galaxies in high-redshift environments. A possible example for such efficient LWB delay of early star formation is the massive Pop III candidate GLIMPSE-16043 at z = 6.5, recently discovered by JWST (Fujimoto et al. 2025).

Our analysis, shown in Fig. 4, quantifies the PBH contribution to J_LW for various masses and fractions. For the simple isothermal and Makino profiles, PBHs with M_PBH = 1 M_⊙ are ruled out for f_PBH ≥ 10⁻², while those with M_PBH = 10 M_⊙, 33 M_⊙, and 100 M_⊙ are excluded for f_PBH ≥ 10⁻³. Importantly, these constraints are reached at z ≳ 25, before molecular cooling ceases to dominate as the primary cooling mechanism in minihalos. The Liu profile provides less restrictive constraints for M_PBH = 1 M_⊙ and 10 M_⊙, ruling out f_PBH ≥ 10⁻¹ and f_PBH ≥ 10⁻², respectively, while remaining consistent with the tighter limits for M_PBH = 33 M_⊙ and 100 M_⊙ (f_PBH > 10⁻³). This relaxation of constraints for the Liu profile stems from the same characteristics discussed in the context of the X-ray background, where its lower gas densities for minihalos at z ≳ 20 lead to reduced accretion luminosities compared to the other profiles.

Fig. 4. Background intensity of LW radiation (JLW) produced by accreting high-z PBH sources. The format of panels, colours, and line styles is identical to Fig. 3: each column corresponds to a different PBH mass, colours represent fPBH values, and line styles denote the density profiles considered (simple isothermal, Makino+1998, and Liu+2022). The grey dashed line represents the critical value JLW ∼ 1, where the destruction rate of H2 equals the formation rate, while the orange dashed-dotted line shows the contribution from stellar sources for reference. The bottom row illustrates the relative contributions of halos ($J_{\mathrm{LW}}^{\mathrm{h}}$) and the IGM ($J_{\mathrm{LW}}^{\mathrm{IGM}}$) to the total LW background intensity ($J_{\mathrm{LW}}^{\mathrm{tot}}$). These are represented as $J_{\mathrm{LW}}^{\mathrm{h}}/J_{\mathrm{LW}}^{\mathrm{tot}}$ (dashed-dotted lines) and $J_{\mathrm{LW}}^{\mathrm{IGM}}/J_{\mathrm{LW}}^{\mathrm{tot}}$ (solid lines), respectively. The relative contributions are computed using the same PBH masses and density profile configurations as in the top row. For fPBH > 0.0001, the halo contribution dominates at all redshifts, while for fPBH = 0.0001, the IGM contribution becomes significant at smaller redshifts, especially for larger PBH masses.

The bottom panels of Fig. 4 illustrate the relative contributions of halos ($J_{\mathrm{LW}}^{\mathrm{h}}$) and the IGM ($J_{\mathrm{LW}}^{\mathrm{IGM}}$) to the total LWB intensity for various f_PBH and M_PBH. For f_PBH = 10⁻⁴, the IGM contribution dominates the LWB over a range of redshifts z < 40 across all PBH masses. However, it transitions to being halo-dominated at lower redshifts in all cases except for M_PBH = 100 M_⊙, where the IGM contribution remains dominant even at z = 6. In this specific scenario, the IGM exceeds the halo contribution by approximately 42 per cent, accounting for 71 per cent of the total LWB intensity compared to 29 per cent from halos. Despite this, the IGM contribution does not influence the derived constraints, as it remains well below the theoretical LWB limit.

For M_PBH = 100 M_⊙ with f_PBH = 10⁻³, there is also a range at z < 40 where the IGM briefly dominates the total intensity. Nonetheless, by z = 6, the halo contribution becomes dominant again, as occurs in all other configurations except for M_PBH = 100 M_⊙ with f_PBH = 10⁻⁴.

The relative increase in the IGM contribution with decreasing redshift, observed for specific PBH masses and fractions, follows the same explanation outlined in the X-ray section. It results from the interplay between the evolving density profiles of halos and accretion dynamics. While this effect is notable, the IGM contribution remains negligible compared to that of halos in all cases that are relevant for setting constraints. Additionally, as in the X-ray background case, for higher PBH masses, the transition from IGM dominance to halo dominance occurs at progressively lower redshifts. This effect can be explained by the interplay between PBH mass, their number density for a given f_PBH, and the evolving density structures of the IGM and halos. For more massive PBHs, the number of such objects is lower due to the fixed fraction of DM they represent, which reduces their cumulative emission in the diffuse IGM. Simultaneously, their higher individual accretion rates in halos compensate for their lower abundance, resulting in a stronger contribution from halos at lower redshifts.

Moreover, while the density profiles of halos evolve and decrease with redshift, the central densities of the most massive halos remain sufficiently high to sustain efficient accretion at lower redshifts. In contrast, the IGM density decreases more rapidly with time, diminishing its relative contribution. This dynamic results in the transition to halo-dominated emission occurring at progressively lower redshifts for higher PBH masses, reflecting the combined effects of hierarchical structure formation, the balance between the number and efficiency of accreting PBHs, and the evolving density profiles of the IGM and halos.

Exceeding the J_LW ∼ 1 threshold at z ≳ 25 would disrupt molecular cooling in minihalos, critically delaying or suppressing the formation of Population III stars. This makes the constraints on f_PBH particularly significant, as molecular cooling remains the dominant mechanism for gas condensation in minihalos until z ∼ 10 − 15, when atomic hydrogen cooling begins to take over in more massive halos (M_h ≳ 10⁸ M_⊙). Any excess contribution to the LWB beyond the allowed limits at high redshifts would drastically alter the timeline of early star formation, delaying the onset of stellar populations and impacting the formation of early galaxies.

By comparing these findings with the constraints derived from the CXB, we observe that the limits in terms of orders of magnitude remain consistent. However, when considering the more refined constraints, an exception arises: for M_PBH = 10 M_⊙ and f_PBH = 0.0008, the LWB intensity exceeds the theoretical limit of J₂₁ ∼ 1 at z ∼ 25, reaching J₂₁ ∼ 1.15. Consequently, the constraint for this mass is tightened to f_PBH ≤ 0.0006.

This refinement highlights the complementarity between the CXB and LWB in probing PBH scenarios, as the general consistency between these independent backgrounds strengthens the robustness of the derived constraints while also providing additional precision in specific cases.

Earlier work by Liu et al. (2022) concluded that the LWB does not impose significant constraints on PBHs. However, we find that this conclusion arises from a coding implementation issue in their modelling, which underestimated the PBH contributions to J_LW. Our analysis indicates that the LWB provides robust constraints on PBH scenarios, particularly in regulating the thermal and chemical history of the early IGM and the formation of primordial stars.

While we have referred to the upper limits on LW background fluxes as strict constraints, a super-critical LW background (J_LW > 1) would, in principle, be possible, shifting the onset of Pop III star formation from minihalos to atomic cooling halos and delaying metal enrichment (Schauer et al. 2021). However, JWST observations suggest that metal enrichment was already in place for some of the highest-z galaxies (Curti et al. 2022), making this scenario less likely. In contrast, observational constraints from X-ray and radio backgrounds impose strict upper limits that cannot be exceeded, as they rely on direct measurements. Meanwhile, the LW background constraints are inferred indirectly from their impact on structure formation. Nonetheless, the constraints derived from the X-ray and LW backgrounds are mutually consistent in order of magnitude.

3.1.3. Radio background

The radio background brightness temperature (T_b) provides crucial insights into the contributions of accreting PBHs to the CRB. Observations from ARCADE 2 reveal a significant excess at low frequencies that cannot be explained by known astrophysical sources (Fixsen et al. 2011; Condon et al. 2012). Specifically, any Bremsstrahlung contributions from high-z HII regions, powered by star formation in the pre-reionisation universe, would be limited to less than 10 per cent of the signal measured by ARCADE 2 (Liu et al. 2019).

At z = 0, Fig. 5 shows the predicted T_b contributions from PBHs accreting in halos. The results indicate that, even in the extreme case of f_PBH = 1, the observed excess is not exceeded for both the Makino+1998 and simple isothermal profiles.

Fig. 5. Brightness temperature of the radio background (Tb) as a function of frequency for PBHs accreting in halos, evaluated at z = 0 (top row) and z = 18 (bottom row). Columns correspond to different PBH masses (MPBH = 1, 10, 33, 100 M⊙), while the line styles (solid and dashed-dotted) represent the simple isothermal and Makino+1998 density profiles, respectively. The colour scheme indicates fPBH: red for 1, green for 10−2, purple for 10−4, and gold for fmax. The maximum fraction of DM allowed for each PBH mass, fmax, is determined from CXB observational and LWB theoretical constraints: 7 × 10−3 for MPBH = 1 M⊙, 6 × 10−4 for MPBH = 10 M⊙ and 33 M⊙, and 7 × 10−4 for MPBH = 100 M⊙. In the top panels, the grey continuous line corresponds to the fit for the observed radio background excess at z = 0 from Condon et al. (2012).

The maximum fractions of PBHs allowed by the CXB (observational) and LWB (theoretical) constraints (f_PBH, max = 7 × 10⁻³ for M_PBH = 1 M_⊙, 6 × 10⁻⁴ for M_PBH = 10 M_⊙ and 33 M_⊙, and 7 × 10⁻⁴ for M_PBH = 100 M_⊙) remain within the observational CRB limits at z = 0. However, their contributions to the CRB excess are negligible, with T_b at 1420 MHz being at most of the order of 1 per cent of the observed excess for the Makino+1998 profile and 10⁻¹ per cent for the simple isothermal profile. These results suggest that PBHs, even if they constitute the entirety of DM, cannot significantly contribute to the CRB at the present epoch.

In the second row of Fig. 5, predictions for T_b at z = 18 are shown. This redshift is of particular interest due to the anomalous 21 cm absorption signal detected by EDGES (Bowman et al. 2018), corresponding to an observed frequency of ν_obs = 74.8 MHz. Explaining the depth of this absorption feature requires an enhanced radio background relative to the cmB, quantified by the parameter A_r. This parameter characterises the ratio of the radio background intensity to the cmB intensity at the relevant frequencies and must satisfy 1.9 < A_r < 418 (Fialkov & Barkana 2019). These constraints are based on the EDGES low-band detection and limits from the Long Wavelength Array 1 (LWA1), which observes the extragalactic radio background at low frequencies (Dowell & Taylor 2018).

As seen in Fig. 5, the T_b values at z = 18 are systematically higher than those at z = 0 across all configurations, reflecting the enhanced accretion rates and denser halo environments at earlier times, where PBHs have a more pronounced impact on the thermal and radiative properties of the gas. Additionally, at z = 18, the differences between the simple isothermal and Makino+1998 profiles become negligible

Table 1 summarises the computed A_r values at z = 18. Configurations that satisfy the EDGES-required range of 1.9 < A_r < 418 (Fialkov & Barkana 2019) correspond only to the extreme cases with f_PBH = 1 for M_PBH > 1. Although these configurations appear viable within the CRB limits alone, they are excluded by the tighter CXB and LWB constraints. Furthermore, the maximum PBH fractions permitted by these multi-frequency constraints fall short of producing the required A_r by 2–4 orders of magnitude. This suggests that, if the EDGES signal is indeed astrophysical in origin, PBHs alone would not be sufficient to account for the enhanced radio background. We also studied the CRB produced by PBHs in the IGM and found that the derived constraints remain unchanged as the contribution within halos is greater.

Ziparo et al. (2022) also studied the CRB and concluded that PBHs are not significant contributors to the observed radio excess. However, their results differ slightly from ours due to differences in the modelling of radio luminosity (L_R). While we calculate L_R based on detailed accretion regimes, including eADAF, standard ADAF, LHAF, and thin discs (as explained in Section 2.2), Ziparo et al. (2022) adopt the fundamental plane relation, which empirically connects radio luminosity, X-ray luminosity, and black hole mass. This approach simplifies the dependency of L_R on the accretion physics, focusing instead on an observationally derived correlation. Despite these differences, both studies suggest that, even with improved modelling, the contribution of PBHs to the CRB is likely to be small.

These results emphasise the complementarity of multi-frequency analyses in constraining PBH parameter space. While T_b predictions at z = 18 offer valuable insights into the early contributions of PBHs to radio backgrounds, constraints from CXB and LWB analyses remain the most stringent, ruling out the configurations that could explain the CRB excess.

3.2. Effects of model variations on inferred backgrounds and PBH abundances

To evaluate the influence of different modelling assumptions on the derived constraints for PBHs, we investigate four distinct variations of the baseline model outlined in Section 2. These variations are designed to emphasise the importance of key assumptions and assess their impact on the predicted cosmic backgrounds. The analysis and results presented in Fig. 6 show both the soft X-ray and Lyman-Werner background intensities generated by PBHs accreting within halos. Each column corresponds to a specific PBH mass–1, 10, 33, and 100 M_⊙–and the f_PBH value derived from the baseline model at which the observational threshold, discussed in Section 3.1.1, or the theoretical threshold, described in Section 3.1.2, is exceeded. The first row shows the cumulative soft X-ray background intensity, while the second row displays the LW background intensity. For all results, accretion in halos is modelled using two gas density profiles: the simple isothermal and Makino+1998 formulations, as indicated in the figure. The third density profile, derived from the rescaled halo in Liu+2022, is not included in this analysis as it does not yield the strongest constraints on f_PBH.

Fig. 6. Cumulative soft X-ray and Lyman-Werner background intensities from PBHs accreting within halos under different model assumptions. Each column corresponds to a specific PBH mass (MPBH = 1 M⊙, 10 M⊙, 33 M⊙, and 100 M⊙) and the fPBH value derived from the baseline model where the observational (nsCXB) or theoretical (JLW ∼ 1) thresholds are exceeded. The top row shows the soft X-ray background intensity (J0.5 − 2 keV), while the bottom row displays the LW background intensity (JLW). The coloured lines represent different model variations: baseline model (blue), Ziparo+2022 velocity (green), standard ADAF without subregimes (yellow), CDM power spectrum (red), and isothermal PBH distribution (purple). Dashed and dashed-dotted lines correspond to results using the simple isothermal and Makino+1998 density profiles, respectively. The grey shaded region in the top row indicates the observed excess X-ray background from Cappelluti et al. (2017), and the dashed grey line in the bottom row marks the critical JLW ∼ 1, where molecular hydrogen dissociation balances its formation.

Below, we outline the key features of each variation:

• Baseline model: As described in Section 2, this model assumes that PBHs constitute DM distributed according to an NFW density profile, with the HMF modified to account for the presence of PBHs. The characteristic velocity is taken as the virial velocity, which reflects the gravitational potential and velocity dispersion within the halo (blue lines in Fig. 6).

• Ziparo+2022 velocities: Instead of the virial velocity, the sound speed is adopted, and the relative velocity between PBHs and baryons is set to zero. The sound speed is expressed as

  $$\begin{array}{c}\hfill c_{\mathrm{s}}={\left(\frac{k_{\mathrm{B}}{T}_{\mathrm{vir}}}{\mu m_{\mathrm{p}}}\right)}^{1/2}\sim 8.3T_{\mathrm{vir},4}^{1/2}\mathrm{km}/\mathrm{s},\end{array}$$(48)

  where T_vir, 4 is the virial temperature in units of 10⁴ K. This assumes hydrodynamical equilibrium, which enhances accretion rates (green lines in Fig. 6).

• Standard ADAF (no sub-regimes): This variation only considers the standard ADAF and thin disc models, excluding the eADAF and LHAF sub-regimes. This leads to substantially weaker X-ray backgrounds, as the exclusion of LHAF contributions removes a key component of accretion luminosity, which, in this case, plays a dominant role in high-density regions within halos (yellow lines in Fig. 6).

• CDM power spectrum: The unmodified CDM HMF is adopted, neglecting the effects of PBHs on the mass function. This reduces sensitivity to f_PBH, though the impact on constraints is within an order of magnitude compared to the baseline model (red lines in Fig. 6).

• Isothermal PBH distribution: Following Zhang et al. (2024a), this variation assumes an isothermal density profile for PBHs instead of the NFW profile used in the baseline model. This results in enhanced accretion rates on smaller scales due to the smoother spatial distribution of PBHs (purple lines in Fig. 6).

The analysis of model variations highlights several key trends and sensitivities in the derived constraints for f_PBH, emphasising the importance of exploring these variations to understand the dependence of our results on different modelling considerations. When adopting the CDM power spectrum and assuming the Makino+1998 density profile, the constraints become moderately less restrictive across all PBH masses, consistently increasing by a factor of approximately 2 compared to the baseline model while remaining within the same order of magnitude. For instance, for M_PBH = 1 M_⊙, the combined X-ray and LW constraint shifts from f_PBH < 10⁻² to f_PBH < 2 × 10⁻². Similarly, for M_PBH = 10 M_⊙, 33 M_⊙, and 100 M_⊙, the constraints change from f_PBH < 10⁻³ to f_PBH < 2 × 10⁻³. Importantly, when using the Simple isothermal profile, the constraints derived under the CDM power spectrum are fully consistent with those obtained with the Makino+1998 profile. This uniform relaxation across all masses reflects the relatively limited impact of the unmodified CDM power spectrum on the HMF at lower f_PBH. As discussed in Section 2.1.2, the HMF becomes more significantly modified only at higher PBH fractions (f_PBH ∼ 1), where isocurvature perturbations dominate. These findings suggest that the PBH-induced changes to the HMF are minimal for the parameter space considered, further supporting the robustness of the baseline model’s assumptions for the HMF modification.

The ADAF without the sub-regimes model produces the most pronounced deviations, particularly for low-mass PBHs, under the Makino+1998 density profile. For M_PBH = 1 M_⊙, the X-ray constraint is relaxed by an entire order of magnitude, shifting from f_PBH < 10⁻² to f_PBH < 3 × 10⁻¹. This result implies that the fraction of DM that could consist of PBHs increases from at most 1–30 per cent. For higher masses (M_PBH ≥ 10 M_⊙), the constraints remain within the same order of magnitude, with the largest variations observed in the LW background. For example, for M_PBH = 33 M_⊙, the LW constraint shifts from f_PBH < 10⁻³ to f_PBH < 5 × 10⁻³, while for M_PBH = 100 M_⊙, the LW constraint changes slightly from f_PBH < 10⁻³ to f_PBH < 4 × 10⁻³. Notably, when the Simple isothermal profile is adopted, the relaxation observed for M_PBH = 1 M_⊙ under the ADAF without the sub-regimes model no longer spans an order of magnitude. Instead, the constraint shifts more modestly from f_PBH < 10⁻² to f_PBH < 2 × 10⁻², consistent with the order of magnitude of the baseline model. This reduced impact is attributable to the higher central gas densities in the Simple isothermal profile, which enhance accretion and mitigate the exclusion of the LHAF sub-regime. For higher PBH masses (M_PBH ≥ 10 M_⊙), the constraints derived with the Simple isothermal profile remain consistent with those obtained using the Makino+1998 profile, showing no significant deviations.

In contrast, the Ziparo+2022 velocity model and the Isothermal PBH distribution yield constraints that are very similar to those of the baseline model, regardless of the chosen density profile. Both models introduce slight tightening of the constraints but do not affect the order of magnitude of the allowed f_PBH. This stability suggests that the relative velocity assumptions and PBH spatial distributions do not significantly influence the inferred constraints, underscoring the robustness of these parameters in shaping the results.

Overall, these findings highlight the importance of modelling assumptions in shaping the constraints on f_PBH, particularly for low-mass PBHs. Notable relaxations arise under the ADAF without the sub-regimes model, especially for M_PBH = 1 M_⊙ with the Makino+1998 profile, where the exclusion of LHAF significantly weakens the constraints. However, the Simple isothermal profile mitigates these effects due to its higher central gas densities, which enhance accretion. Variations in the CDM power spectrum consistently produce moderate relaxations (factor of 2 across all masses) without changing the order of magnitude, reflecting the limited impact of HMF modifications at low f_PBH. By contrast, variations in velocity assumptions (Ziparo+2022 velocity model) and PBH spatial distributions (Isothermal PBH distribution) show minimal deviations, indicating that these specific assumptions have less influence within the explored parameter space. This analysis validates our results independently of these assumptions and underscores the need to consider such variations to ensure the robustness of the derived constraints.

As stated earlier, all results presented in this work assume a monochromatic PBH mass function. This choice facilitates comparison with previous studies and enables a detailed treatment of accretion processes without requiring additional assumptions about the mass distribution. While such an approximation is well motivated in scenarios where PBHs form from narrow features in the primordial power spectrum, broader mass distributions may weaken constraints by spreading the emission across multiple scales and modifying PBH clustering. Although a full treatment of extended mass functions lies beyond the scope of this paper, it remains an important avenue for future exploration. Constraints on PBH abundance derived under the monochromatic assumption can be reinterpreted for extended mass functions using the formalism of Bellomo et al. (2018).

Figure 7 summarises the main results of this work, showing the maximum allowed DM fraction in the form of PBHs, f_PBH, max, as a function of monochromatic PBH mass. The figure includes our constraints for different model variations and two halo gas density profiles, along with a comparison to existing observational limits and to the independent theoretical bound from Ziparo et al. (2022). Results are shown for a simple isothermal gas density profile (left panel) and for the Makino+1998 profile (right panel). In general, the allowed fraction decreases with all model variations converging to the same f_PBH, max for each individual mass at M_PBH ≥ 10⁴ M_⊙, where the IGM contribution dominates due to the decreasing number of halos massive enough to host at least one PBH. Notably, between 33 M_⊙ and 100 M_⊙, f_PBH, max remains roughly constant, which can be attributed to a reduction in the total number of PBHs at these higher masses, compensating for the increased mass of individual PBHs in this specific mass range.

Fig. 7. Constraints on the maximum allowed PBH DM fraction for different monochromatic PBH masses, including results derived in this work under various modelling assumptions, as well as observational limits and previous theoretical results. We consider two gas density profiles within halos: a simple isothermal (left panel) and the Makino+1998 profile (right panel). For each case, we compute the allowed PBH abundance under different modelling scenarios: baseline (blue squares), standard ADAF without subregimes (yellow circles), isothermal PBH spatial distribution (purple hexagons), Ziparo+2022 velocity prescription (green downward triangles), and CDM power spectrum (red upward triangles). Maximum fractions are computed at one-decade intervals in mass across the range 10−8 M⊙ to 106 M⊙. Each coloured marker with a black edge corresponds to the maximum allowed PBH fraction obtained by applying the constraint methodology developed in this work, under a specific model variation. Solid lines connecting the markers and shaded regions, all shown in the same colour as their corresponding markers, are included to aid visual interpretation: although intermediate values were not explicitly computed, the region above each evaluated point can be reliably considered excluded under the corresponding assumptions. In both panels, our constraints are compared with observational limits, with the results from (solid black line; Ziparo et al. 2022), which represent independent theoretical bounds derived from modelling PBH emission contributing to the cosmic X-ray background using a different formalism, and with the results from CV2024 (black stars). While all other constraints correspond to upper bounds on fPBH, the CV2024 markers indicate fractions for which significant PBH-induced feedback effects are already present in the gas within halos at z ∼ 23. Observational constraints, all derived assuming a monochromatic PBH mass function, include microlensing limits from the Subaru M31 survey (Niikura et al. 2019a) updated by Kusenko et al. (2020), green dashed line), from the Expérience pour la Recherche d’Objets Sombres/MAssive Compact Halo Objects (EROS-2/MACHO) (Tisserand et al. 2007, green dot-dashed line), from the OGLE 5 yr survey (Niikura et al. 2019b, dotted green line), and from the OGLE 20 yr survey (Mróz et al. 2024, solid green line). A dotted black line shows the 95% confidence region assuming the six ultrashort-timescale OGLE events are due to planetary-mass PBHs (Niikura et al. 2019b). Additional constraints include the PBH fraction inferred under the assumption that all BH mergers detected during the third observing run (O3) of the LIGO-Virgo-KAGRA Collaboration originate from PBHs (Wong et al. 2021, blue point with error bars), subsolar-mass PBH binary merger limits from LIGO O2 (Abbott et al. 2019, solid blue line), cmB accretion limits from Poulin et al. (2017, orange dashed line), and Galactic Centre constraints from X-ray (magenta solid) and radio (purple dashed) observations (Manshanden et al. 2019).

In the isothermal gas profile case, all model variants yield f_PBH, max = 1 for M_PBH ≲ 10⁻² M_⊙, indicating that our methodology does not constrain PBHs as DM at these scales. However, several observational limits are already more stringent in this mass range. At M_PBH ∼ 10⁻¹ M_⊙, the baseline and CDM power spectrum models provide constraints comparable to EROS-2/MACHO, while the standard ADAF model (excluding LHAF and eADAF sub-regimes) yields results consistent with the OGLE 5 yr limit. The models in which PBHs follow an isothermal distribution, or adopt the velocity prescription from Ziparo et al. (2022), yield f_PBH, max(M_PBH = 10⁻¹ M_⊙)∼0.005 − 0.01, below the EROS-2/MACHO and OGLE 5 yr limits.

For the Makino+1998 profile, the constraints are slightly weaker, mainly due to lower central gas densities. Both the standard ADAF (no sub-regimes) and CDM power spectrum variants allow f_PBH = 1 for M_PBH ≲ 10⁻¹ M_⊙, while the baseline, isothermal PBH distribution, and Ziparo+2022 velocity models yield f_PBH, max(M_PBH = 10⁻¹ M_⊙)∼0.1, remaining above the OGLE 20 yr and EROS-2/MACHO bounds. At M_PBH ≳ 10³ M_⊙, all models converge to f_PBH, max ≲ 10⁻⁵, independently of the modelling assumptions.

We also compare our results with the theoretical constraint from Ziparo et al. (2022), shown as a solid black line. Although the overall shape is similar, the underlying physical assumptions differ substantially. Ziparo et al. (2022) compute the PBH contribution to the unresolved cosmic X-ray background using a semi-analytic approach based on integrating over the halo population. They adopt the Makino+1998 gas profile and estimate accretion rates for each halo mass bin, assuming a fixed spectral shape for all PBHs: a power law with exponential cutoff at high energies, inspired by ADAF models. The bolometric-to-X-ray conversion is performed using a constant correction factor f_X = 0.1 in the 2 − 10  keV band, and no distinction is made between accretion regimes or halo-specific physical conditions. In contrast, our methodology includes regime-dependent spectral modelling, incorporating transitions between ADAF, LHAF, and eADAF, and constructs halo-by-halo spectral energy distributions that depend explicitly on gas density, velocity dispersion, and accretion rate. Despite these differences, both approaches robustly exclude dominant PBH contributions at M_PBH ≳ 10² M_⊙, with f_PBH, max ≲ 0.001, reinforcing the strength of constraints derived from cosmic backgrounds in this mass range. However, our approach produces more detailed constraints and their dependence on different assumptions.

Our results also showed good agreement with existing multi-messenger bounds across several mass ranges. At M_PBH ∼ 1 M_⊙, the baseline and standard ADAF-only models in the isothermal gas profile case yield values consistent with the LIGO O2 subsolar merger limit (Abbott et al. 2019) and the OGLE 20 yr microlensing constraint (Mróz et al. 2024), respectively. At M_PBH ∼ 100 M_⊙, our constraints match the CMB accretion bounds from Poulin et al. (2017). Among all model variants, the baseline case provides the most restrictive limits across the full mass range, particularly for M_PBH ∼ 1 − 10 M_⊙, where it falls below all observational and theoretical bounds considered. This highlighted the importance of including sub-regime transitions in the accretion model and accounting for halo-by-halo physical variation when assessing PBH constraints from cosmic backgrounds.

In addition to the limits discussed above, our results can be directly compared with the PBH fractions identified in our previous study, CV2024, shown as black stars in Figure 7. That analysis focused on the impact of PBHs within individual halos at z ∼ 23, applying the same accretion and emission model used here, which includes all ADAF subregimes (eADAF, ADAF, LHAF) and the thin disc, but implemented in post-processing to halos extracted from one of the simulations in the CIELO suite (Tissera et al. 2025). Gas densities were modelled using a rescaled high-resolution profile from Liu et al. (2022), while characteristic velocities were extracted directly from the simulation particles. Unlike all other constraints shown in the figure, which correspond to maximum allowed values of f_PBH, the CV2024 markers represent the values at which significant feedback effects are already present in the gas. These include temperature increases by a factor of around 300 relative to the no-PBH case, with gas temperatures reaching up to 1.7 T_vir, and a suppression of the central neutral hydrogen abundance to below 1% of its unperturbed value. Although the two approaches differ in scale and methodology, they yield consistent results. For both the isothermal and Makino+1998 profiles, the baseline models presented in this work yield maximum allowed PBH fractions that lie just below the thresholds identified in CV2024. This agreement across independent frameworks reinforces the robustness of our findings and highlights the complementarity between local feedback signatures and large-scale limits in constraining the PBH abundance.

Finally, we verified that the emission from individual halos remains below current observational detection thresholds and is therefore not removed from the background through source masking. This result holds across all model variations. As a representative example, we considered PBHs with mass 33, M⊙, constituting all the DM, embedded in halos of 10⁶, M⊙ at z = 30. This scenario provides a conservative upper bound, as the total emission is dominated by halos with masses below 10⁶, M_⊙, as discussed in Section 2.3. Depending on the model variation, we find total luminosities ranging from approximately 10³⁵–10³⁹ erg/s in the 0.5–2 keV X-ray band, and from approximately 10³⁷–10³⁹ erg/s in the 0.044–2.07 eV infrared band. These values remain well below the detection thresholds of current deep field surveys, supporting the validity of comparing our predictions to the unresolved cosmic backgrounds (Cappelluti et al. 2017; Kaminsky et al. 2025).

4. Conclusion

Our study placed stringent constraints on the maximum fraction of DM that PBHs could constitute, and quantified their contributions to the CXB, LWB, and CRB across different mass ranges. A detailed summary is provided below for reference PBH masses:

• 1 M_⊙: If PBHs constituted 7 × 10⁻³ of the DM, they would fully explain the CXB background (99 per cent). A larger fraction would exceed the observed values. On the other hand, such an abundance of 1 M_⊙ PBHs only explained 33 per cent of the hard CXB. Their contribution to the LWB was significant but remained below the critical threshold (J_LW ∼ 1) needed to suppress star formation. Contributions to the CRB were negligible.

• 10 M_⊙: PBHs were constrained to a maximum DM fraction of 6 × 10⁻⁴. In this case, they could contribute up to 93 per cent of the soft CXB and 37 per cent of the hard CXB. Similarly, their LWB contribution was non-negligible but remained consistent with the J_LW ∼ 1 threshold, ensuring no significant disruption of molecular cooling. Their impact on the CRB remained insignificant.

• 33 M_⊙: PBHs could comprise up to 6 × 10⁻⁴ of DM. They would contribute in this case up to 80 per cent of the soft CXB and 33 per cent of the hard CXB. As with lower masses, their LWB contributions were limited by the J_LW ∼ 1 threshold, and their CRB contributions were negligible.

• 100 M_⊙: PBHs were limited to a maximum DM fraction of 7 × 10⁻⁴. These could explain up to 91 per cent of the observed soft CXB excess and 39 per cent of the hard CXB. Their LWB contributions were significant but adhered to the J_LW ∼ 1 limit, while their CRB contributions remained negligible.

Our study emphasised the importance of detailed modelling and the exploration of variations in key assumptions. We adopted an accretion model that included distinct subregimes–eADAF, standard ADAF, LHAF, and thin discs–which allowed us to capture the diversity of PBH accretion scenarios based on local halo and IGM conditions. Notably, including these subregimes relaxed the constraints on low-mass PBHs (1 M_⊙) by up to an order of magnitude compared to models without them, due to the enhanced accretion luminosities captured by this approach. We also tested the robustness of our results against variations in gas density profiles, velocity prescriptions, emission models, and the PBH impact on the halo mass function.

Building on this framework, we extended our analysis across the full mass range 10⁻⁸ − 10⁶ M_⊙, enabling a direct comparison with a wide array of observational and theoretical constraints. In particular, we found that our baseline model yielded constraints more restrictive than all existing bounds in the intermediate mass range M_PBH ∼ 1 − 10 M_⊙, for both density profiles considered. These were slightly stronger than those of Ziparo et al. (2022), though within the same order of magnitude. These results also demonstrated that assumptions regarding accretion subregimes, PBH spatial distributions, or velocity prescriptions could shift the allowed f_PBH values by factors of a few, with the most extreme case relaxing constraints by up to an order of magnitude. Altogether, these findings reinforced the value of a comprehensive and physically grounded modelling approach when evaluating the viability of PBHs as dark matter candidates.

While our results aligned with many previous studies, the combination of detailed accretion modelling, multi-frequency analysis, and systematic exploration of model variations provided a more nuanced understanding of PBHs as DM candidates. This approach refined the allowed parameter space for PBHs and emphasised the role of robust modelling in interpreting their astrophysical implications.

Our constraints, derived from their contributions to multiple cosmic radiation backgrounds, indicated that PBHs cannot constitute the totality of DM. However, their potential to uniquely impact early cosmic history further motivates their in-depth investigation, as their astrophysical signatures could provide key insights into both their nature and the broader evolution of the universe.

Acknowledgments

NP acknowledges support from Préstamo BID PICT Raices 2023 N° 0002. PBT acknowledges partial funding by Fondecyt-ANID 1240465/2024 and N°cleo Milenio ERIS NCN2021_017. BL gratefully acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foun- dation) under Germany’s Excellence Strategy EXC 2181/1 – 390900948 (the Heidelberg STRUCTURES Ex- cellence Cluster). This project has received funding from the European Union Horizon 2020 Research and Innovation Programme under Marie Skłodowska-Curie Actions (MSCA) grant agreement No. 101086388-LACEGAL. We acknowledge partial support by ANID BASAL project FB210003. We also thank Prof. Dr. Günther Hasinger for his valuable comments and suggestions, which have significantly improved this manuscript.

References

Appendix A: Derivation of radiation background equations

In this appendix, we present the derivation of the equations used to calculate the radiation background intensity arising from PBHs accreting within halos and the IGM, considering the cumulative contribution of sources distributed among cosmic structures.

The physical (proper) flux δF, defined as the energy received per unit time and per unit area by an observer at redshift z, is given by

$$\begin{array}{c}\hfill F=\frac{\delta E}{\delta A\delta t},\end{array}$$(A.1)

where δE represents the energy of a photon packet passing through an area δA over a time interval δt in the observer frame at redshift z. In the case of a diffuse source, the energy and time intervals in the source frame are modified due to cosmological redshift as

$$\begin{array}{c}\hfill \delta E^{\prime }=\frac{1+z\prime }{1+z}\delta E,\delta t^{\prime }=\frac{1+z}{1+z\prime }\delta t.\end{array}$$(A.2)

Thus, the observed flux can be rewritten as

$$\begin{array}{c}\hfill F={\left(\frac{1+z}{1+z\prime }\right)}^2\frac{1}{\delta A}\frac{\delta E^{\prime }}{\delta t^{\prime }}.\end{array}$$(A.3)

For a diffuse radiation field, the luminosity can be expressed in terms of the source frame physical emissivity ϵ(z′), which quantifies the energy emitted per unit volume and per unit time. Considering that the radiation originates from a diffuse medium occupying a physical volume of δV′=δa′δℓ′ at redshift z′, where δℓ′ is the line element and δa′ represents the cross-sectional area of the cylinder spanned by the light rays connecting the source to the observer, the luminosity takes the form

$$\begin{array}{c}\hfill \frac{\delta E^{\prime }}{\delta t^{\prime }}=ϵ(z\prime )\delta V^{\prime }=ϵ(z\prime )\delta a\prime \delta {\ell }^{\prime }.\end{array}$$(A.4)

Considering the cosmic expansion, the relationship between the observed cross-sectional area and the comoving area is given by

$$\begin{array}{c}\hfill \delta a\prime ={\left(\frac{1+z}{1+z\prime }\right)}^2\delta A.\end{array}$$(A.5)

Substituting these relations into Eq. (A.3), the flux for a diffuse source can be expressed as

$$\begin{array}{c}\hfill F={\left(\frac{1+z}{1+z\prime }\right)}^4ϵ(z\prime )\delta {\ell }^{\prime }.\end{array}$$(A.6)

The total flux, obtained by integrating the contributions from sources along the line of sight up to a maximum redshift z_max, is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill F& ={\displaystyle {\int }_z^{z_{max}}}{\left(\frac{1+z}{1+z\prime }\right)}^4ϵ(z\prime )\frac{d{\ell }^{\prime }}{dz\prime }dz\prime \hfill \\ \hfill & =-{\displaystyle {\int }_z^{z_{max}}}{\left(\frac{1+z}{1+z\prime }\right)}^4ϵ(z\prime )c\frac{d{t}_{\mathrm{U}}}{dz\prime }dz\prime \hfill \end{array}\end{array}$$(A.7)

where the physical line element dℓ′ is related to the cosmic age t_U through the relation dℓ′/dz′= − c dt_U/dz′.

To account for frequency dependence, the specific flux (energy per unit time, per unit area, per unit frequency) in the observer frame is defined as

$$\begin{array}{c}\hfill F_{\nu }\equiv \frac{\mathit{dF}}{d\nu },\end{array}$$(A.8)

where the observed frequency ν and the source-frame frequency ν′ are related by

$$\begin{array}{c}\hfill \nu =\frac{1+z}{1+z\prime }\nu \prime .\end{array}$$(A.9)

The specific emissivity of the diffuse source is then defined as

$$\begin{array}{c}\hfill {ϵ}_{\nu \prime }=\frac{dϵ}{d\nu \prime }.\end{array}$$(A.10)

Substituting these definitions, the total frequency-dependent flux from diffuse sources is expressed as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill F_{\nu }& =-{\displaystyle {\int }_z^{z_{\max }}}{\left(\frac{1+z}{1+z\prime }\right)}^4\frac{dϵ(z\prime )}{d\nu \prime }\frac{d\nu \prime }{d\nu }c\frac{d{t}_U}{dz\prime }dz\prime \hfill \\ \hfill & =-{\displaystyle {\int }_z^{z_{\max }}}{\left(\frac{1+z}{1+z\prime }\right)}^3{ϵ}_{\nu \prime }(z\prime )\frac{cd{t}_U}{dz\prime }dz\prime .\hfill \end{array}\end{array}$$(A.11)

In addition to cosmic age t_U, the flux equations can be reformulated in terms of the comoving distance r, defined from z = 0 to z′. Given the relation dr/dz′= − (1 + z′)c dt_U/dz′, the flux can be rewritten as

$$\begin{array}{c}\hfill F={\displaystyle {\int }_z^{z_{max}}}\frac{{(1+z)}^4}{{(1+z\prime )}^5}ϵ(z\prime )\frac{\mathit{dr}}{dz\prime }dz\prime ,\end{array}$$(A.12)

$$\begin{array}{c}\hfill F_{\nu }={\displaystyle {\int }_z^{z_{max}}}\frac{{(1+z)}^3}{{(1+z\prime )}^4}{ϵ}_{\nu \prime }(z\prime )\frac{\mathit{dr}}{dz\prime }dz\prime .\end{array}$$(A.13)

Finally, the radiation intensity J, defined as the energy per unit time, per unit area, and per unit solid angle, is related to the flux by the expression J = F/(4π). When considering frequency dependence, an additional dependence per unit frequency is introduced, leading to the definition of the specific intensity J_ν, which is related to the specific flux by J_ν = F_ν/(4π).

To derive Eq. (43), we begin with Eq. (A.12) and use the relation ϵ(z′) = ϵ_cm(z′)(1 + z′)³, where the comoving integrated emissivity ϵ_cm(z′) over a frequency range [ν₁, ν₂] is given by

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{cm},[{\nu }_1,{\nu }_2]}(z,z\prime )={\displaystyle {\int }_{{\nu }_1(z,z\prime )}^{{\nu }_2(z,z\prime )}}\frac{{ϵ}_{\nu \prime }(z\prime )}{{(1+z\prime )}^3}d\nu \prime ={\displaystyle {\int }_{{\nu }_1(z,z\prime )}^{{\nu }_2(z,z\prime )}}{\dot{\rho }}_{\nu \prime }(z\prime )d\nu \prime .\end{array}$$(A.14)

Substituting this definition into Eq. (A.12) yields

$$\begin{array}{c}\hfill J=\frac{1}{4\pi }{\displaystyle {\int }_z^{z_{max}}}\frac{{(1+z)}^4}{{(1+z\prime )}^2}{ϵ}_{\mathrm{cm}}(z\prime )\frac{\mathit{dr}}{dz\prime }dz\prime .\end{array}$$(A.15)

Evaluating the expression at z = 0 results in the cancellation of the factor (1 + z)⁴. By discretising the integral, we obtain Eq. (43).

To derive Eq. (44), we start from Eq. (A.7) and apply the relation ϵ(z′) = ϵ_cm(z′)(1 + z′)³. Substituting the cosmic time differential dt_U/dz′ with

$$\begin{array}{c}\hfill \frac{d{t}_{\mathrm{U}}}{dz\prime }=\frac{d{t}_{\mathrm{U}}}{da\prime }(-\frac{1}{{(1+z\prime )}^2}),\end{array}$$(A.16)

we arrive at the following expression for the background intensity

$$\begin{array}{c}\hfill J=\frac{1}{4\pi }{\displaystyle {\int }_z^{z_{max}}}\frac{{(1+z)}^4}{{(1+z\prime )}^3}{ϵ}_{\mathrm{cm}}(z\prime )c{\frac{d{t}_{\mathrm{U}}}{da\prime }|}_{z\prime }dz\prime .\end{array}$$(A.17)

Similarly, to derive Eq. (46), we start from Eq. (A.11) and incorporate the comoving specific emissivity ${\dot{\rho }}_{\nu \prime }(z\prime )$ along with Eq. (A.16). Substituting these expressions results in

$$\begin{array}{c}\hfill J_{\nu }=\frac{1}{4\pi }{\displaystyle {\int }_z^{z_{max}}}\frac{{(1+z)}^3}{{(1+z\prime )}^2}{\dot{\rho }}_{\nu \prime }(z\prime )c{\frac{d{t}_{\mathrm{U}}}{da\prime }|}_{z\prime }dz\prime .\end{array}$$(A.18)

All Tables

All Figures

	Fig. 1. Impact of PBHs on the power spectrum and HMF (assuming MPBH = 100 M⊙). Top panel: total power spectrum Ptot(k) as a function of wavenumber k for various fractions of DM in PBHs (fPBH = 0.0001, 0.001, 0.01, 0.1, and 1), compared to the standard ΛCDM case (fPBH = 0). Bottom panel: HMF at z = 10, showing the number density of halos as a function of halo mass M, also assuming MPBH = 100 M⊙.
In the text

Fig. 2. Comparison of hydrogen number density profiles nH(r) within DM halos of different masses (Mh = 104 M⊙, 106 M⊙, and 108 M⊙) at redshifts z = 6 (solid lines), z = 10 (dashed-dotted lines), and z = 40 (dotted lines). The three density profiles examined include a simple isothermal profile (magenta), the profile derived by Makino et al. (1998) (blue), and the rescaled profile from Liu et al. (2022) (green). The vertical grey lines denote the virial radii of the corresponding halos at each redshift.

In the text

Fig. 3. Contribution of accreting high-z BH sources to the present-day (z = 0) soft X-ray background (0.5–2 keV). The top row shows the cumulative evolution of the integrated background intensity, JX, as a function of redshift z, for PBH masses of MPBH = 1 M⊙, 10 M⊙, 33 M⊙, and 100 M⊙ in each column. Different values of fPBH are represented by distinct colours, as indicated in the legend, ranging from 10−4 to 1. The grey shaded region corresponds to the observed excess X-ray background reported in Cappelluti et al. (2017). Three halo density profiles are considered: simple isothermal (dashed lines), the profile from Makino et al. (1998) (dashed-dotted lines), and the profile from Liu et al. (2022) (dotted lines). The bottom row illustrates the relative contributions of halos ($J_{\mathrm{X}}^{\mathrm{h}}$) and the IGM ($J_{\mathrm{X}}^{\mathrm{IGM}}$) to the total X-ray background intensity ($J_{\mathrm{X}}^{\mathrm{tot}}$). These are represented as $J_{\mathrm{X}}^{\mathrm{h}}/J_{\mathrm{X}}^{\mathrm{tot}}$ (dashed-dotted lines) and $J_{\mathrm{X}}^{\mathrm{IGM}}/J_{\mathrm{X}}^{\mathrm{tot}}$ (solid lines), respectively. The relative contributions are computed using the same PBH masses and density profile configurations as in the top row.

In the text

Fig. 4. Background intensity of LW radiation (JLW) produced by accreting high-z PBH sources. The format of panels, colours, and line styles is identical to Fig. 3: each column corresponds to a different PBH mass, colours represent fPBH values, and line styles denote the density profiles considered (simple isothermal, Makino+1998, and Liu+2022). The grey dashed line represents the critical value JLW ∼ 1, where the destruction rate of H2 equals the formation rate, while the orange dashed-dotted line shows the contribution from stellar sources for reference. The bottom row illustrates the relative contributions of halos ($J_{\mathrm{LW}}^{\mathrm{h}}$) and the IGM ($J_{\mathrm{LW}}^{\mathrm{IGM}}$) to the total LW background intensity ($J_{\mathrm{LW}}^{\mathrm{tot}}$). These are represented as $J_{\mathrm{LW}}^{\mathrm{h}}/J_{\mathrm{LW}}^{\mathrm{tot}}$ (dashed-dotted lines) and $J_{\mathrm{LW}}^{\mathrm{IGM}}/J_{\mathrm{LW}}^{\mathrm{tot}}$ (solid lines), respectively. The relative contributions are computed using the same PBH masses and density profile configurations as in the top row. For fPBH > 0.0001, the halo contribution dominates at all redshifts, while for fPBH = 0.0001, the IGM contribution becomes significant at smaller redshifts, especially for larger PBH masses.

In the text

Fig. 5. Brightness temperature of the radio background (Tb) as a function of frequency for PBHs accreting in halos, evaluated at z = 0 (top row) and z = 18 (bottom row). Columns correspond to different PBH masses (MPBH = 1, 10, 33, 100 M⊙), while the line styles (solid and dashed-dotted) represent the simple isothermal and Makino+1998 density profiles, respectively. The colour scheme indicates fPBH: red for 1, green for 10−2, purple for 10−4, and gold for fmax. The maximum fraction of DM allowed for each PBH mass, fmax, is determined from CXB observational and LWB theoretical constraints: 7 × 10−3 for MPBH = 1 M⊙, 6 × 10−4 for MPBH = 10 M⊙ and 33 M⊙, and 7 × 10−4 for MPBH = 100 M⊙. In the top panels, the grey continuous line corresponds to the fit for the observed radio background excess at z = 0 from Condon et al. (2012).

In the text

Fig. 6. Cumulative soft X-ray and Lyman-Werner background intensities from PBHs accreting within halos under different model assumptions. Each column corresponds to a specific PBH mass (MPBH = 1 M⊙, 10 M⊙, 33 M⊙, and 100 M⊙) and the fPBH value derived from the baseline model where the observational (nsCXB) or theoretical (JLW ∼ 1) thresholds are exceeded. The top row shows the soft X-ray background intensity (J0.5 − 2 keV), while the bottom row displays the LW background intensity (JLW). The coloured lines represent different model variations: baseline model (blue), Ziparo+2022 velocity (green), standard ADAF without subregimes (yellow), CDM power spectrum (red), and isothermal PBH distribution (purple). Dashed and dashed-dotted lines correspond to results using the simple isothermal and Makino+1998 density profiles, respectively. The grey shaded region in the top row indicates the observed excess X-ray background from Cappelluti et al. (2017), and the dashed grey line in the bottom row marks the critical JLW ∼ 1, where molecular hydrogen dissociation balances its formation.

In the text

Fig. 7. Constraints on the maximum allowed PBH DM fraction for different monochromatic PBH masses, including results derived in this work under various modelling assumptions, as well as observational limits and previous theoretical results. We consider two gas density profiles within halos: a simple isothermal (left panel) and the Makino+1998 profile (right panel). For each case, we compute the allowed PBH abundance under different modelling scenarios: baseline (blue squares), standard ADAF without subregimes (yellow circles), isothermal PBH spatial distribution (purple hexagons), Ziparo+2022 velocity prescription (green downward triangles), and CDM power spectrum (red upward triangles). Maximum fractions are computed at one-decade intervals in mass across the range 10−8 M⊙ to 106 M⊙. Each coloured marker with a black edge corresponds to the maximum allowed PBH fraction obtained by applying the constraint methodology developed in this work, under a specific model variation. Solid lines connecting the markers and shaded regions, all shown in the same colour as their corresponding markers, are included to aid visual interpretation: although intermediate values were not explicitly computed, the region above each evaluated point can be reliably considered excluded under the corresponding assumptions. In both panels, our constraints are compared with observational limits, with the results from (solid black line; Ziparo et al. 2022), which represent independent theoretical bounds derived from modelling PBH emission contributing to the cosmic X-ray background using a different formalism, and with the results from CV2024 (black stars). While all other constraints correspond to upper bounds on fPBH, the CV2024 markers indicate fractions for which significant PBH-induced feedback effects are already present in the gas within halos at z ∼ 23. Observational constraints, all derived assuming a monochromatic PBH mass function, include microlensing limits from the Subaru M31 survey (Niikura et al. 2019a) updated by Kusenko et al. (2020), green dashed line), from the Expérience pour la Recherche d’Objets Sombres/MAssive Compact Halo Objects (EROS-2/MACHO) (Tisserand et al. 2007, green dot-dashed line), from the OGLE 5 yr survey (Niikura et al. 2019b, dotted green line), and from the OGLE 20 yr survey (Mróz et al. 2024, solid green line). A dotted black line shows the 95% confidence region assuming the six ultrashort-timescale OGLE events are due to planetary-mass PBHs (Niikura et al. 2019b). Additional constraints include the PBH fraction inferred under the assumption that all BH mergers detected during the third observing run (O3) of the LIGO-Virgo-KAGRA Collaboration originate from PBHs (Wong et al. 2021, blue point with error bars), subsolar-mass PBH binary merger limits from LIGO O2 (Abbott et al. 2019, solid blue line), cmB accretion limits from Poulin et al. (2017, orange dashed line), and Galactic Centre constraints from X-ray (magenta solid) and radio (purple dashed) observations (Manshanden et al. 2019).

In the text