Primordial Black Holes as Near Infrared Background sources

The near infrared background (NIRB) is the collective light from unresolved sources observed in the band 1-10 $\mu$m. The measured NIRB angular power spectrum on angular scales $\theta \gtrsim 1$ arcmin exceeds by roughly two order of magnitudes predictions from known galaxy populations. The nature of the sources producing these fluctuations is still unknown. Here we test primordial black holes (PBHs) as sources of the NIRB excess. Considering PBHs as a cold dark matter (DM) component, we model the emission of gas accreting onto PBHs in a cosmological framework. We account for both accretion in the intergalactic medium (IGM) and in DM haloes. We self consistently derive the IGM temperature evolution, considering ionization and heating due to X-ray emission from PBHs. Besides $\Lambda$CDM, we consider a model that accounts for the modification of the linear matter power spectrum due to the presence of PBHs; we also explore two PBH mass distributions, i.e. a $\delta$-function and a lognormal distribution. For each model, we compute the mean intensity and the angular power spectrum of the NIRB produced by PBHs with mass 1-$10^3~\mathrm{M}_{\odot}$. In the limiting case in which the entirety of DM is made of PBHs, the PBH emission contributes $<1$ per cent to the observed NIRB fluctuations. This value decreases to $<0.1$ per cent if current constraints on the abundance of PBHs are taken into account. We conclude that PBHs are ruled out as substantial contributors to the NIRB.


INTRODUCTION
The Near Infrared Background (NIRB) is the diffuse radiation of cosmological origin observed after subtracting the local foregrounds in the band 1-10 m (Kashlinsky et al. 2018).Since early studies by Partridge & Peebles (1967), the NIRB has been considered a valuable tool to investigate the emission from the first stars and galaxy populations, as ultraviolet (UV) and optical light from high-z sources is redshifted to the near-infrared band.
Actual measurements of the mean NIRB intensity (Tsumura et al. 2013;Matsumoto et al. 2015;Sano et al. 2015;Matsuura et al. 2017) give a lower bound  ≳ 10 nW m −2 sr −1 , in excess with respect to the contribution of known galaxy populations derived from galaxy number counts (Driver et al. 2016).However, direct measurements of the NIRB suffer from large uncertainties due to the subtraction of foregrounds (Leinert et al. 1998), namely interplanetary dust emission (zodiacal light), galactic stars light and galactic interstellar medium radiation (cirrus).
Being foregrounds smooth, a more robust technique is computing the power spectrum of NIRB fluctuations (Kashlinsky et al. 1996;Kashlinsky & Odenwald 2000), to which foregrounds contribution is limited.Moreover, from the power spectrum measurements, a lower limit to the   contribution from unknown sources can be derived (Kashlinsky et al. 2007).The latest measurements of the NIRB power spectrum (Cooray et al. 2012a;Kashlinsky et al. 2012) established an excess power on scales larger then ≳ 1 arcmin, irreconcilable with ★ E-mail: daniele.manzoni@sns.itemission from known galaxies up to  ∼ 5 (Helgason et al. 2012).
The origin of such a signal is still unknown.
Population III stars (PopIII) were one of the first hypothesis proposed about the sources of the NIRB excess (Santos et al. 2002;Salvaterra & Ferrara 2003).Although intriguing, this idea was soon after discarded because of the very high formation efficiency required (Madau & Silk 2005) and since it would overpredict the number of high-z dropout galaxies (Salvaterra & Ferrara 2006).Several works explored the possibility of high redshift galaxies ( ≳ 5) being the sources of the NIRB excess, but models failed to reproduce the required levels of fluctuations (Fernandez et al. 2010;Cooray et al. 2012b;Yue et al. 2013a;Helgason et al. 2016).
An alternative solution was the intrahalo light (IHL), i.e light from stars stripped from their parent galaxy (Cooray et al. 2012a;Cheng & Bock 2022).Despite its success in reproducing observations, such idea has to rely on poorly understood abundance of intrahalo stars (Ferrara 2012).Moreover, this model cannot account for the observed cross-correlation of the NIRB with the soft-X background (SXB) (Cappelluti et al. 2013(Cappelluti et al. , 2017)).Such a feature is difficult to explain even with galaxies spectra, but could be naturally justified by X-ray emission from accretion disks around black holes.Yue et al. (2013b) developed a model to explain both NIRB fluctuations and NIRB-SXB cross correlation with accreting direct collapse black holes (DCBHs), even though it is unclear whether the specific conditions of DCBHs formation are actually realized during cosmic evolution (Latif & Ferrara 2016).
Given the puzzling nature of the NIRB excess, a new scenario has been recently suggested, invoking Primordial Black Holes (PBHs) (Kashlinsky 2016;Cappelluti et al. 2022).PBHs are black holes formed deep into radiation dominated era from the collapse of overdensity peaks (Zel'dovich & Novikov 1967;Carr & Hawking 1974) and interest on them have been rejuvenated after the first detection of gravitational waves from black holes merger (Abbott et al. 2016;Bird et al. 2016;Blinnikov et al. 2016;Sasaki et al. 2016).The primordial origin of LIGO/VIRGO black holes is a viable solution to explain their observed mass spectrum and merger rates (Raidal et al. 2017;Ali-Haïmoud et al. 2017;Wong et al. 2021).Moreover, they could justify why most of the measured effective spins are close to zero (Abbott et al. 2019;De Luca et al. 2020) and could accomodate for black holes with masses in the pair-instability supernovae mass gap (45-120 M ⊙ ) (Abbott et al. 2020a;De Luca et al. 2021;O'Brien et al. 2021) and in the low mass gap (2.5-5 M ⊙ ) (Abbott et al. 2020c,b;Clesse & García-Bellido 2022).Finally, the recent evidence of a gravitational-wave background reported by the NANOGrav collaboration (Agazie et al. 2023) could directly probe PBH formation from high amplitude peaks of the primordial power spectrum (Clesse & Garcìa-Bellido 2017;Vaskonen & Veermäe 2021;Franciolini et al. 2023).
A key aspect about PBHs is that they were proposed as cold dark matter candidates (Chapline 1975).This hypothesis has been investigated in a plethora of studies, providing constraints on the fraction of DM comprised by PBHs (see Carr & Kühnel (2020) for a review).The presence of PBHs would entail a variety of astrophysical phenomena, such as gamma rays emission from evaporating PBHs (Laha 2019;Coogan et al. 2021), microlensing effects (Niikura et al. 2019;Blaineau et al. 2022) and disruption of wide binaries or ultra-faint dwarfs (Monroy-Rodríguez & Allen 2014;Brandt 2016).In addition, accreting PBHs would impact the CMB spectrum and anisotropies (Poulin et al. 2017;Serpico et al. 2020), the 21 cm power spectrum (Mena et al. 2019) and would produce radio and X-ray backgrounds (Cappelluti et al. 2022;Ziparo et al. 2022).
When deriving constraints on the abundance of PBHs, it is commonly assumed that PBHs have the same mass (i.e. a -function), although these constraints actually depend on the adopted PBH mass function (Kühnel & Freese 2017).In particular, PBH formation models in slow-roll inflation predict an approximately lognormal mass function (Dolgov & Silk 1993;Kannike et al. 2017), while latest simulations of PBH formation across the QCD epoch derived a mass function peaked around  PBH ∼ 1 M ⊙ , with a non-trivial shape departing from lognormal (Franciolini et al. 2022).
If PBHs constitute a fraction of dark matter, they would add a poissonian component to the matter power spectrum (Meszaros 1975;Afshordi et al. 2003;Ali-Haïmoud 2018), accelerating structure formation and consequently enhancing the abundance of haloes in which stars can form (Kashlinsky 2016).This effect on the star formation process is particularly relevant for what concerns the NIRB excess puzzle, since a higher star formation rate density at high-z can then provide the required levels of NIRB fluctuations (Cappelluti et al. 2022).Moreover, PBHs could directly contribute to the NIRB with the radiation emitted by accreting gas from their surroundings.
Hasinger (2020, hereafter H20) computed cosmic backgrounds from gas accretion onto PBHs and could recover only 0.3 per thousand of the NIRB with his model.However, H20 considered gas accretion only in the intergalactic medium (IGM), while PBHs could accrete matter also in dense virialized structures, i.e DM haloes.In particular, Ziparo et al. (2022, hereafter Z22) have shown that the contribution of PBHs accreting in DM haloes to X-ray and Radio backgrounds is > 60 per cent larger than those accreting in the IGM.
In this paper, following the model by Z22, we compute the NIRB produced by PBHs taking into account both PBH accretion in DM haloes and a self-consistent treatment of X-ray ionization and heating of the IGM.We further improve the Z22 model both considering the modification of the matter power spectrum induced by the presence of the PBHs, previously neglected, and generalizing the framework to extended mass functions.In Section 2 we summarize the basic model and present its extensions.In Section 3 we present the main results of this work.Finally, we state our conclusions in Section 4.

METHODS
To investigate the contribution of PBHs to the NIRB we rely on the formalism described in Z22.We first revisit their model in order to introduce the framework (Section 2.1, 2.2).We then compute the intensity and angular power spectrum of the NIRB in Section 2.3.In the last two Sections we extend the model to account for the modification of matter power spectrum induced by PBHs (Section 2.4 ) and extended mass functions (Section 2.5).

Cosmological distribution of PBHs
Assume that a DM fraction  PBH is made of PBHs of mass  PBH .DM distribution on cosmological scales can be described as a diffuse component with density equal to the mean DM density, and virialized regions where matter has collapsed into DM haloes.As PBHs are distributed as the DM, we decompose the number density of PBHs as where  IGM PBH ( ℎ PBH ) is the number density of PBHs in the intergalactic medium (haloes).The abundance of PBHs in haloes is related to the collapsed fraction of DM in haloes  coll , which can be computed as in the Press-Schechter formalism (Press & Schechter 1974).Here  crit () = 1.68/ () is the critical overdensity for collapse,  () is the growth factor and  2  is the mass variance of the linearly extrapolated matter overdensity field.Thus, the number density of PBHs in the IGM and in haloes are (3)

PBH distribution inside haloes
The distribution of PBHs inside haloes follows the DM density profile, here assumed to be NFW (Navarro et al. 1996): where  = / vir is the radial distance in virial radius units and  is the concentration parameter from Macciò et al. (2007).Following Z22, we model its redshift evolution as  ∝ (1 + ) −1 .The parameter   is a function of both the concentration parameter and the overdensity at the collapse redshift Δ  (Barkana & Loeb 2001): with Being PBHs distributed as DM, the number of PBHs within radius  an  + d is

PBHs accretion
To estimate the accretion rate of gas onto PBHs, we adopt the Bondi-Hoyle-Lyttleton formula (Bondi 1952;Edgar 2004): where   and   are the density and sound speed of the accreting gas, respectively,  BH is the relative velocity between the PBH and the gas, and  is the accretion parameter that accounts for non gravitational effects (i.e radiative feedback, gas pressure, outflows).Following Poulin et al. (2017), we adopt the value  = 0.01, which is a benchmark for an advection dominated accretion flow (ADAF, Yuan & Narayan 2014).Accretion conditions in the IGM and inside haloes differ substantially: in the following we describe the relevant physical quantities, i.e.   ,   and  BH , separately for the two cases.

Accretion in the IGM
Following Ricotti et al. (2008), we assume a uniform gas density in the IGM, equal to where  = 1.22 is the mean molecular weight for a gas of primordial composition and   is the proton mass.The sound speed of the gas is given by: where   is the Boltzmann constant and  IGM is the IGM temperature.The relative velocity between baryons and PBHs is gaussianly distributed on linear scales, hence its modulus follows a maxwellian distribution, with variance given by (Ali-Haïmoud & Kamionkowski 2017) To properly account for the distribution of relative velocities, it is useful to define an effective velocity  eff (Ricotti et al. 2008), whose analytical expression is (Mena et al. 2019) where  (, , ) is the confluent hypergeometric function of second kind.The accretion rate of PBHs in the IGM is finally obtained by substituting the relevant quantities computed above in equation ( 7).

Accretion within haloes
To model the internal structure of haloes, we assume that the gas is in thermal equilibrium at the virial temperature  vir .Moreover, we impose hydrostatic equilibrium between DM and gas.Given these assumptions, the density profile of gas is described by the following equation (Makino et al. 1998): where  esc is the escape velocity, given by: and  ,0 is a normalization constant set by imposing: 4 where  ℎ is the halo mass, Ω  and Ω DM are the total baryon and DM densities in units of the critical density.The sound speed in haloes can be computed via equation ( 9), substituting  IGM with  vir .As a consequence of hydrostatic equilibrium assumption, we set  BH = 0.

NIRB
To compute the specific luminosity of PBHs we follow Z22.Given the accretion rate , the bolometric luminosity of a single PBH is  =   2 , where  = 0.1 is the radiative efficiency.We assume that, as for astrophysical black holes, the spectrum of PBHs can be described by a double power-law with an exponential cut-off (H20): where the cut-off frequency is  cut = 200 keV and  = −0.7.Below the critical frequency   =   /, with   = 0.45( PBH /M ⊙ ) 0.4 m, synchrotron emission dominates and the power law index is  sync = 1.86 (H20).The above spectral shape is consistent with an ADAF accretion model with accretion rates  = /  EDD ≳ 10 −2 , which holds for those PBHs producing the bulk of the background radiation in our model.We fix the normalization of the spectrum by setting the bolometric correction in the 2-10 keV band to   = 0.1 (H20).Given the specific luminosity,   , the specific emissivity of a population of PBHs accreting in the IGM is The specific luminosity of an entire halo can be computed by: The specific emissivity of a population of PBHs accreting inside haloes is then given by integrating over the halo mass function (Murray et al. 2013): where  max =  ℎ ( vir = 10 4 ) is the minimum mass of haloes inside which stars can form and  min is the minimum mass of haloes required to form a baryon overdensity (Barkana & Loeb 2001): Inside haloes with  ℎ <  min , the gas density is close to the mean IGM one and therefore we consider their contribution in the IGM emissivity.
The background intensity in a given band [ 1 , 2 ] is related to the specific emissivity by (Fernandez et al. 2010;Yue et al. 2013a): where  ′ = (1 + ) and  () is the Hubble parameter as a function of redshift.The angular power spectrum of NIRB fluctuations from PBHs can be decomposed in a two-halo and a shot-noise term: The clustering component at frequency  and for the multiple moment  is given by (Cooray et al. 2004;Fernandez et al. 2010) where  () is the comoving distance and (, ) is the power spectrum of the underlying matter distribution.PBHs in the IGM correspond to DM in the linear regime and therefore  IGM (, ) =  lin (, ), where the right hand side is the linear matter power spectrum.Instead, haloes are biased tracers of the linear matter density field and their power spectrum can be written as  ℎ (, ) =  eff ()(, ), where the effective bias  eff is given by: where  ℎ (, ) is the halo bias, as derived in Tinker et al. (2010).The shot noise angular power spectrum is described by the following equation (Cooray et al. 2012a;Yue et al. 2013a): We note that in principle one should consider the one-halo term, given by (Cooray et al. 2012a): where ũ( = / (), ) is the Fourier transform of the NFW profile.
For the redshift and halo mass range of interest, we checked that ũ( = / (), ) ∼ 1 and therefore the one-halo term reduces to the shot noise term in equation ( 24).

Matter Power spectrum modified by PBHs
PBHs may constitute a fraction of DM, thus they would add a Poisson shot noise term to the linear matter power spectrum (Meszaros 1975;Afshordi et al. 2003): where  PBH,0 is the PBH number density at redshift  = 0.The total matter power spectrum can be then written as (Villanueva-Domingo & Ichiki 2023): where  () is the linear growth factor and  iso is the isocurvature transfer function.An approximate expression for  iso is given by (Peacock 1998): 10 4 10 5 10 6 10 7 10 8 10 9 M where  eq is the redshift of radiation-matter equality and  eq =  −1  ( eq )/(1 +  eq ).The contribution to the power spectrum from PBHs can be recast in the form (Villanueva-Domingo & Ichiki 2023): The PBH modification to the power spectrum affects the variance of the matter overdensity field and thus the halo mass function.Hereafter we will refer to a PBH-ΛCDM cosmology whenever adopting the power spectrum described by equations ( 28)-(30).In particular we consider a PBH-ΛCDM cosmology in our models PBH- and PBH-lognormal (see Sec. 3).In Fig. 1 we show the Press-Schechter halo mass function at  = 20 in the standard ΛCDM scenario and including the modification induced by PBHs, for different values of the parameter  PBH  PBH .At  = 20, when including the extra power on small scales due to PBHs, the halo mass function is a factor of 3 (40) higher for  ℎ = 10 5 M ⊙ ( ℎ = 10 7 M ⊙ ) with respect to the standard ΛCDM case, considering  PBH  PBH = 100 .In Fig. 2 we compare the bolometric emissivity, from both haloes and IGM, in the ΛCDM and PBH-ΛCDM cosmologies, as a function of redshift.As a consequence of the increased number of small haloes expected in the PBH-ΛCDM, the contribution to the total emissivity from accreting PBHs in DM haloes is enhanced by a factor of 2 (20) at redshift  = 30 (40).Moreover, also the collapsed DM fraction is higher and thus the relative contribution from PBHs accreting in the IGM is further lowered.The emissivity of PBHs accreting in haloes at redshift  = 20 (40) is roughly 10 (100) times the emissivity from PBHs in the IGM.We point out that, as a consequence of the aforementioned effects, in the PBH-ΛCDM cosmology halo emissivity dominates the IGM one at any redshift, unlike in the standard ΛCDM case.

Extended PBH mass function
The mass function of PBHs at the epoch of their formation is denoted by ( PBH ), and defined as: In the following, we generalize our formalism to extended mass functions.In particular, we consider the case of a lognormal mass function 1 : Here   is the critical mass that sets the position of the peak and  is the standard deviation of the distribution.Regarding PBHs in the IGM, their emissivity can be generalized to where the integral is performed over the PBH mass and  IGM PBH (, ) is taken from equation (3).Recalling that  IGM PBH ∝  −1 and   ∝  2 , we can write Therefore, when computing the emissivity of PBHs in the IGM, an extended mass function is equivalent to a -function centered at the mean mass M of the mass function 2 .For a lognormal mass function, the mean mass is Mlog =   exp  2 /2 .We will assume the benchmark value  = 1 throughout the rest of the paper and quote only the mean mass of the lognormal distribution.
Regarding PBHs accreting in haloes, we must specify how PBHs of different masses are distributed inside the halo.We note that the frictional acceleration exerted onto a body of mass  moving through a homogeneous distribution of particles of mass  ( ≪ ) with isotropic velocity distribution is ∝  (Binney & Tremaine 2008).Hence, PBHs with higher masses sink towards the centre of the halo before lighter ones.With this in mind and for simplicity, we then assume that more massive PBHs lie at smaller radii. 1 We choose a lognormal mass function to avoid fruitless complications.The main results of our work are unaffected by the exact shape of the mass function.
2 This is valid for a constant radiation efficiency.If  ∝   ∝  2 , then   ∝  2+2 and so the corresponding mean mass should be M = ∫ d  (  )  2+1 1/(2+1) .The mass  * () of PBHs at a given radius  can be derived by imposing that the mass  () enclosed in a sphere of radius  is equal to the integrated mass of all PBHs more massive than  * (): For a lognormal mass function, the right-hand side of the above equation can be computed analitically, giving where  ℎ is the halo mass.We show the resulting  * () for a lognormal mass function with Mlog = 30 M ⊙ in Fig. 3. Once specified  * (), the number of PBHs per unit length at radius  is given by d  d = 4    * () We can then substitute equation (37) into equation ( 17) and apply the same formalism described in Sec.2.3.In Fig. 4 we compare the bolometric luminosity of haloes in the case of a delta mass function with  PBH = 30 M ⊙ and of a lognormal mass function with Mlog = 30  ⊙ .Including the lognormal mass function boosts the halo luminosity by a factor of ∼ 10, because more massive PBHs accrete at smaller distances from the center, where the gas density is higher.

RESULTS
In this section we present the IGM temperature evolution, the mean NIRB intensity and the NIRB angular power spectrum obtained from three different models: (i) standard ΛCDM cosmology with a PBH delta mass function (ΛCDM-); (ii) PBH-ΛCDM cosmology with PBH delta mass function (PBH-) and (iii) PBH-ΛCDM cosmology with a PBH lognormal mass function (PBH-lognormal).We compare our predictions to observational data to test the hypothesis of accreting PBHs as sources of the NIRB.

IGM temperature and ionization evolution
X-ray emission from PBHs would heat and ionize the IGM well before galaxies start to reionize the Universe.To account for this ef- fect, we self-consistently derive the IGM temperature and ionization evolution following the formalism described in Sec. 3 of Z22.
In the pre-overlap phase of the cosmic reionization process, the Universe can be split in ionized and neutral regions.In the ionized regions, the redshift evolution of the free electron fraction   () is solved through equation (33a) in Z22, adopting the photoionization rate derived from the UV background in Puchwein et al. (2019).The free electron fraction traces the evolution of the volume filling factor of ionised regions, namely the fraction of volume occupied by ionized regions.In the same regions, we assume an IGM temperature  IGM,ion = 10 4 K.Such high temperature suppresses accretion onto PBHs due to high sound speeds (equation ( 9)).
In neutral regions, whose volume filling factor is 1−  (), the free electron fraction  , () evolves with redshift according to equation (35a) in Z22.Here, the photoionization rate calculation accounts for secondary ionizations due to X-rays emitted by PBHs: where ℎ min = 0.5 keV3 ,  th is the hydrogen ionization threshold,   is the hydrogen ionization cross section and  ion ∼ 0.3 is the fraction of the primary electron's energy going into secondary ionizations (Furlanetto & Stoever 2010).In neutral regions we also solve the redshift evolution of temperature  IGM,n () through equation (35b) in Z22, which takes into account the IGM heating due to the energy injected by X-rays (Mesinger et al. 2013), here assumed to be emitted by PBHs.The heating rate  PBH per baryon can then be computed as: where  heat ∼ 0.3 is the fraction of the primary electron's energy going into heat (Valdés et al. 2010).
Following the evolution of neutral regions  IGM,n is crucial because changes in the IGM temperature affects PBHs emissivity.On the one hand, if the IGM temperature increases, the effective velocity of PBHs accreting in the IGM increases as well (equation ( 11)): this lowers their luminosity and consequently their emissivity.On the other hand, to higher  IGM,n correspond higher  min (equation ( 19)): the integration interval in equation ( 18) is thus shortened, which reduces the emissivity of PBHs accreting in haloes.Therefore, if  IGM,n increases (decreases), the total emissivity of PBHs is lowered (enhanced).
We show the resulting temperature evolution in Fig. 5 for the three different models, adopting  PBH = Mlog = 30 M ⊙ and  PBH = 10 −4 , 10 −3 , 10 −2 , 10 −1 , 1.For  PBH = 1, heating from PBHs increases the IGM temperature in neutral regions at  ∼ 6 by a factor of ∼ 10, 20, 60 in model ΛCDM-, PBH-, PBH-lognormal respectively, with respect to a ΛCDM cosmology which does not include PBHs.While in model ΛCDM-  IGM,n starts increasing around  ∼ 30, in model PBH- it rises at higher redshifts ( ∼ 40), as the PBH emissivity is boosted by the higher number of small mass ( ℎ ≲ 10 6 -10 7 M ⊙ ) haloes.For  PBH ≲ 10 −2 , the effect of PBHs on the halo mass function is negligible and the evolution of  IGM,n in the two cases is almost identical.In model PBH-lognormal, the luminosity of haloes is further enhanced by the lognormal mass function (Fig. 4) and  IGM,n reaches ∼ 900K at  ∼ 6.
Before moving to the core results of this work, we briefly comment on the implications of IGM heating from PBHs.Firstly, our model does not affect the IGM temperature at  ≲ 5, where measurements from Lyman-alpha forest observations are obtained (Walther et al. 2019;Gaikwad et al. 2020).Below  ∼ 6, most of the Universe is ionized and IGM temperatures  IGM,ion ≳ 10 4 suppress the emission and thus the heating from PBHs in ionized regions.Moreover, the contribution from PBHs in neutral regions is also suppressed because their volume filling factor, i.e 1 −   in our model, approaches zero as cosmic reionization proceeds.Instead, at  ≳ 10, the Universe is mostly neutral and the volume filling factor of neutral regions is basically unity.Radiation from PBHs is then effective in heating the IGM above the adiabatic cooling temperature.Therefore, forthcoming 21 cm observations could provide stringent constraints on the heating from PBHs and hence on their abundance (Mena et al. 2019).

NIRB mean intensity and angular power spectrum
As already mentioned in the Introduction, direct measurements of the mean NIRB intensity are uncertain (Kashlinsky et al. 2018).However, constraints on the mean NIRB excess from unknown sources can be derived from the angular power spectrum measurement, as done by Kashlinsky et al. (2007), who found  2−5 m ≳ 1 nW m −2 sr −1 in the 2-5 m band.

NIRB mean intensity
We compute the mean NIRB intensity in the 2-5 m band, using equation (20).To test our predictions, in Fig. 6 we compare these results with data from Kashlinsky et al. (2007).We consider the limiting case  PBH = 1, which sets the upper limit for NIR flux produced by PBHs, for a mass range 1 M ⊙ ≤ M PBH = Mlog ≤ 10 3 M ⊙ .Our choice is driven, one the one hand, by the requirement of substantial accretion rates and hence luminosities and, on the other hand, by existing constraints on the abundance of higher mass PBHs.We show the results for the three different models in Fig. 6.We find that in the ΛCDM- (PBH-, PBH-lognormal) model, PBHs of cold dark matter (DM), we have computed the mean intensity and angular power spectrum of the NIRB arising from their accretion.Following the formalism by Ziparo et al. (2022), we account for PBH accretion both in the intergalactic medium (IGM) and in DM haloes, and we self-consistently derive the IGM temperature evolution, considering ionization and heating due to X-ray emission from PBHs.The Z22 model is based on the ΛCDM linear matter power spectrum, and considers a  function for the PBH mass distribution.
Besides this ΛCDM- model, we have considered the possibility that PBHs modify the matter power spectrum (PBH- model), and follow an extended lognormal mass function (PBH-lognormal model).In both PBH- and PBH-lognormal models we adopt a PBH-ΛCDM cosmology, accounting for the matter power spectrum modified by PBHs.
For each model, we have derived the intensity and angular power spectrum of the NIRB finding that PBHs contribute to the observed NIRB fluctuations to < 1 per cent, even in the most optimistic cases considered in this work.This conclusion is supported by these intermediate results: • The PBH modification to the power spectrum affects the variance of the matter overdensity field and thus the halo mass function, adding an extra power on small scales.In particular, at  = 20, in the PBH-ΛCDM model, the halo mass function is a factor of 3 (40) higher for  ℎ = 10 5 M ⊙ ( ℎ = 10 7 M ⊙ ) with respect to the standard ΛCDM case, considering  PBH  PBH = 100 .
• As a consequence of the increased number of small haloes expected in the PBH-ΛCDM, the contribution to the total emissivity from accreting PBHs in DM haloes is enhanced by a factor of 2 (20) at redshift  = 30 (40).Moreover, also the collapsed DM fraction is higher and thus the relative contribution from PBHs accreting in the IGM is further lowered.The emissivity of PBHs accreting in haloes at redshift  = 20 (40) is roughly 10 (100) times the emissivity from PBHs in the IGM.We point out that, as a consequence of the aforementioned effects, in the PBH-ΛCDM cosmology the halo emissivity dominates the IGM one at any redshift, unlike in the standard ΛCDM case.
• If the radiative efficiency is indipendent of the accretion rate, given an extended mass function, the emissivity of PBHs accreting in the IGM can be computed adopting a delta mass function with mass equal to the mean mass of the mass function Mlog .We compare the bolometric luminosity of haloes in the case of a delta mass function with  PBH = 30 M ⊙ and of a lognormal mass function with Mlog = 30  ⊙ .Including the lognormal mass function boosts the halo luminosity by a factor of ∼ 10, because more massive PBHs accrete at smaller distances from the center, where the gas density is higher.
• Considering 1 ≤  PBH [M ⊙ ] ≤ 10 3 and  PBH = 1, PBHs can produce at most ∼ 1 per cent of the flux required to explain NIRB fluctuations.The three models differ in their prediction by less than a factor of ≃ 2. Although in PBH-ΛCDM cosmology the total emissivity of PBHs at  ≳ 40 is ∼ 10 higher than in the standard ΛCDM scenario (for both delta and lognormal mass functions), the resulting NIRB is similar, because the gas heating from X-rays produced by PBHs damps their emissivity at lower redshifts.
• When accounting for current constraints on PBH abundance, the maximum relative contribution of PBHs to the NIRB is reduced to 0.1 per cent, for PBHs with  PBH ∼ 50 M ⊙ .
• None of our models is able to reproduce the NIRB angular power spectrum.At large angular scales ( ∼ 20 arcmin), fluctuations predicted by model ΛCDM- (PBH-, PBH-lognormal) are lower than the measured one by a factor of 1000 (400,200), in the most favorable case with  PBH = 10 3 M ⊙ .
Before concluding, we compare our findings with the results from Hasinger (2020, hereafter H20), whose model is adopted in Cappelluti et al. (2022).H20 predicted a NIR flux from PBHs of 10 −13 erg s −1 cm −2 deg −2 ∼ 3 × 10 −4 nW m −2 sr −1 .This corresponds to ∼ 0.3 per thousand of the NIRB flux required to explain NIRB fluctuations.Hence, we find a NIRB flux ∼ 10× higher then the one obtained in H20.We point out some substantial differences between the two models to fully grasp the discrepancy in the two results.Firstly, H20 adopt a non-linear relative velocity between gas and DM to capture the collapse of baryons into DM haloes.This approach does not account for the density profile of DM and gas inside haloes, which enhance the contribution of PBHs accreting in haloes, as gas densities are much higher than the mean baryon density.Secondly, H20 estimates as negligible the heating of accreting gas by X-rays produced by PBHs, which instead in our model provides a negative feedback on the PBH emissivity.Moreover, H20 adopts an extended mass function with a peak around 1 M ⊙ , but with broad tails reaching up to 10 9 M ⊙ .They conclude that the dominant contribution arises from PBHs with  PBH ∼ 10 4 M ⊙ , while we focused only on the range 1 ≤  PBH /M ⊙ ≤ 10 3 .A final difference concerns the accretion parameter, which they assume to be  = 0.05, i.e. 5× higher than the one adopted here.
To summarize, even if our modelling for the PBH contribution to the NIRB excess differs from the H20 one, we don't end up with a dramatic discrepancy.This is because the extra physical effects that we have included tend to balance each other.In fact, we should have expected a much higher NIRB flux due to the contribution of PBHs accreting in haloes and the boosted matter power spectrum due to the presence of the PBHs.However, these effects are balanced by the inclusion of the IGM heating from PBH X-ray emission that damps their emissivity at lower redshifts.

Figure 2 .
Figure2.Bolometric emissivity from IGM (dashed lines) and haloes (solid lines), in the standard ΛCDM scenario (black) and including the power spectrum modified by PBHs (green).In the case of PBH-ΛCDM cosmology, the halo signal dominates over the IGM one for all the redshift of interest.Here we adopt  PBH = 30 M ⊙ .

Figure 4 .
Figure 4. Bolometric luminosity of haloes as a function of the halo mass, in the case of a delta (black) and lognormal (blue) mass function, for  =20, 30, 40 (solid, dashed, dotted respectively).The PBH mass for the delta mass function and the mean mass of the lognormal distribution are both set to  PBH = Mlog = 30  ⊙ .With an extended mass function, the luminosity is boosted by a factor ∼ 10, because more massive PBH tend to sink towars the center where the gas density is higher.