Signatures of primordial black holes as seeds of supermassive black holes

It is broadly accepted that Supermassive Black Holes (SMBHs) are located in the centers of most massive galaxies, although there is still no convincing scenario for the origin of their massive seeds. It has been suggested that primordial black holes (PBHs) of masses $\gtrsim 10^{2} M_\odot$ may provide such seeds, which would grow to become SMBHs. We suggest an observational test to constrain this hypothesis: gas accretion around PBHs during the cosmic dark ages powers the emission of high energy photons which would modify the spin temperature as measured by 21cm Intensity Mapping (IM) observations. We model and compute their contribution to the standard sky-averaged signal and power spectrum of 21cm IM, accounting for its substructure and angular dependence for the first time. If PBHs exist, the sky-averaged 21cm IM signal in absorption would be higher, while we expect an increase in the power spectrum for $\ell~\gtrsim 10^2-10^3$. We also forecast PBH detectability and measurement errors in the abundance and Eddington ratios for different fiducial parameter configurations for various future experiments, ranging from SKA to a futuristic radio array on the dark side of the Moon. While the SKA could provide a detection, only a more ambitious experiment would provide accurate measurements.


I. INTRODUCTION
The idea that density fluctuations can provide the seeds for galaxy formation via gravitational instability and leave detectable traces in the Cosmic Microwave Background (CMB) [1,2] introduced the concept that the graininess in the Universe would be the seeds around which galaxies form [3]. Now we know that Supermassive Black Holes (SMBHs) inhabit the centers of most galaxies (see [4] for a review). Observations of quasars at z ∼ 6 − 7 indicate that, even at these early times, there were SMBHs with masses of several 10 9 M [5][6][7][8]. The existence of a population of Intermediate Mass Black Holes (IMBH) of masses around 10 2 − 10 6 M at z ∼ 20 − 15 would suffice [9] to seed them. The possible detection of a ∼ 10 5 M black hole in the Milky Way close to its center [10,11] may provide evidence for such a relic and support the argument that Intermediate Mass Black Holes are the seeds of SMBHs. Besides, IMBHs may inhabit the center of dwarf galaxies (e.g., [12] and references therein).
The optimal conditions in the relevant parameter space of the mass of the black hole and the gas density around it, that lead to fast growth of the black hole, were studied in [13,14]. This happens if the combined effects of the angular momentum and radiation pressure are ineffective in stopping the stream of gas flowing from large scales Marie Sk lodowska-Curie fellow towards the black hole. They find that this condition is fulfilled for M 10 4 − 10 5 M (where M stands for the mass of the seed) and large gas densities, for which the growth of massive seeds up to SMBHs is feasible. Even so, there is a limit on the maximum mass that the black hole can reach in an isolated halo, which depends on the total mass of the host halo and on the radiative efficiency of the accretion [15]. For smaller masses, the accretion is very inefficient, but fast enough growth can be achieved via mergers.
However, the origin and formation mechanism of the massive seeds are still uncertain (see [16,17] for a review). There are two main scenarios proposed to explain their origin: supercritical growth from stellar mass black holes formed from Population III stars and directly formed massive seeds at lower redshift. There are more exotic scenarios, such as IMBHs formed by a subdominant component of the dark matter being dissipative [18].
According to the first hypothesis, the seeds of SMBHs are remnants of Population III stars, formed with masses of tens of solar masses at z 20, which grow due to gas accretion and mergers [19][20][21]. However, in order to reach masses such as those observed at z ∼ 6 − 7 [22][23][24], the accretion needs to be supercritical over extended periods of time. Moreover, SMBH seeds growth is probably depressed due to the shallow gravitational potentials existing at those redshifts and the radiation pressure of the black hole emission. Indeed the recently discovered IMBHs in dwarfs are anorexic: apparently undermassive compared to the M BH − σ scaling relation [25,26] Besides, cosmic X-ray background observations impose con-straints on the growth of SMBHs, constraining the abundance of quasars with supercritical accretion [27] as well as of the abundance of miniquasars at high redshift [28]. Therefore, this scenario alone is very unlikely to account for the present abundance of SMBHs.
On the other hand, SMBH seeds might also be formed due to the collapse of gas clouds which do not fragment or form ordinary stars, but directly form a massive black hole (M ∼ 10 5 − 10 6 M ) at lower redshifts (z 15) [29][30][31][32][33][34][35][36]. This kind of seed is called a Direct Collapse Black Hole (DCBH). DCBHs may be realized if a metal-poor cloud is irradiated by non-ionizing ultraviolet light from nearby star-forming galaxies, which photodissociate molecular hydrogen and therefore prevent star formation. Hence, the gas can only cool via Lyman-α emission, which leads to a quasi-isothermal contraction without fragmentation until the gravitational collapse and the formation of an IMBH (see e.g. [37]). Conveniently, the DCBH radiation is very efficient in preventing the formation of H 2 . Therefore, a DCBH may trigger the formation of other DCBHs in a slowly-collapsing gas cloud more efficiently than galaxies [38].
Moreover, DCBHs are a good candidate for explaining the large-scale power spectrum of the Near Infrared Background and its cross correlation with the cosmic Xray background [39]. As DCBHs have a characteristic observational imprint [40], it can be possible to identify these seeds in deep multi-wavelength surveys [41]. Two promising candidates, whose infrared spectra require an exceptionally high star formation rate, were found at high redshift, with a predicted mass higher than 10 5 M [41]. These candidates are likely to be formed by direct collapse.
Nonetheless, the exact conditions and the probability of obtaining DCBHs are still uncertain; recent theoretical studies suggest that this mechanism might explain the abundance of the most luminous quasars at z ∼ 6 − 7, but not the general population of SMBHs [42][43][44].
In summary, neither of these two scenarios individually provide an entirely convincing explanation for the origin of the seeds of SMBHs. However, massive seeds could have been formed much earlier. This third possibility (see [45][46][47][48][49][50] and references therein), much less explored in the literature, considers Primordial Black Holes (PBHs) as the seeds which will grow to become SMBHs. If PBHs are formed with large enough masses, there is no need for supercritical accretion, as is the case for Population III stars.
The idea of the existence of PBHs [51] has recently regained popularity after they were suggested to be the progenitors (and to make up a sizeable fraction of the dark matter, see e.g. [52]) of the stellar mass black holes (∼ 30M ) detected by LIGO+VIRGO Collaboration [53].
Since then, a number of possible tests of the model have been performed with available data. They cover all of the theoretically allowed range, from the smallest masses constrained by black holes evaporation [54], to e.g. microlensing of stars [55,56], to larger masses, constrained by e.g. X-ray and radio emission [57], widebinaries disruption [58], and accretion effects [59][60][61]. More innovative tests that can be performed in the future have also been suggested, including using quantum gravity effects [62], the lensing of fast radio bursts [63], the cross-correlation of gravitational waves with galaxy maps [64][65][66], eccentricity of the binary orbits [67], the black hole mass function [68], the gravitational wave mass spectrum [69], merger rates [70] and the stochastic gravitational wave background [71][72][73][74].
The mass range required for PBHs to be the seeds of SMBHs is 10 2 M . In this range, the PBH abundance, f PBH = Ω PBH /Ω CDM , is strongly constrained by e.g., CMB observations [59,61], Ultra-Faint Dwarf Galaxies [75] and wide binaries [58]. However, most of these constraints have been derived in the context of a model in which PBHs comprise most of the dark matter, and they assume a delta function in their mass distribution (see [76] for updated constraints allowing for a wide mass distribution); if, on the other hand, PBHs of these masses are only required to be the seeds of SMBHs and not a substantial part of the dark matter, the high-mass tail of the PBH mass distribution can have a very small f PBH , satisfying all observational constraints.
Different scenarios for the SMBH seeds have different observational signatures. In this paper, we focus on their imprints on 21 cm Intensity Mapping (IM). The term 'IM' is sometimes dropped in the literature related to the emission from neutral hydrogen at large redshifts, in contrast with the studies of the emission lines from galaxies. However, we maintain it for the sake of clarity. 21 cm IM observations represent a promising future tool for cosmology (for a recent review, see [77]). In particular, observations of spin temperature maps in the dark ages provide a direct window into the matter density fluctuations free of complications such as galaxy bias and most astrophysical processes. It can be thought of as a series of CMB-like screens, and therefore, besides the auto-correlation signal, one can also consider the ISW effect [78] and lensing of 21 cm IM maps [79], including the possibility of performing tomographic analyses. It has recently been shown that 21 cm IM observations will give very powerful constraints on e.g., primordial non-gaussianity [80], inflationary models [81][82][83], scattering between dark matter and baryons [84], statistical isotropy [85] and annihilating and decaying dark matter [86,94].
21 cm IM will also be very powerful for setting observational constraints on PBHs. Using the power spectrum originating from Poisson fluctuations, the authors of [87] forecast constraints on f PBH based on future observations with SKA in the mass range M 10 −2 M . The abundance of PBHs of much lighter masses can be constrained by looking for the effects of Hawking evaporation on the Inter Galactic Medium (IGM) via 21 cm IM of the dark ages [88]. Minihalos have been also studied as interesting 21 cm emitters between reionization and the dark ages [89], although they are hard to differentiate from the standard diffuse signal emanating from the IGM [90].
The 21 cm IM sky-averaged signal can be also used to discriminate between the two main scenarios for the origin of SMBH seeds. Seeds formed from remnants of Popularion III stars dominate X-ray heating of the IGM and cause a rise in the 21 cm brightness temperature at z 20. An absence of such a signature might be due to the seeds being formed later, which would favor the DCBH scenario [91]. However, such a signature could originate not only via seeds formed from Population III star remnants but also by PBHs. Besides, at these redshifts, the 21 cm IM signal is affected by a large number of astrophysical uncertainties and its dependence on redshift changes considerably with different assumptions [92,93], making it difficult to identify the signal coming from the SMBH seeds.
Here we study the scenario in which PBHs are the seeds of SMBHs. In order to avoid the astrophysical uncertainties mentioned above, we concentrate on the dark ages (z 30). The detection of a signal corresponding to the predictions reported in this work would be an indication that massive miniquasars were already present in the dark ages and the most likely explanation would be that these black holes are primordial, hence the most straightforward candidate to be the seeds of SMBHs. In the standard scenario, during the so-called dark ages, the cosmic time previous to the formation of the first stars, there is no astrophysical feedback which contaminates the 21 cm IM signal and haloes are still not formed, so observations are free from galaxy bias and non-linearities in the clustering. Hence, the main uncertainties are only coming from the PBH sector. However, other exotic energy injections, such as that sourced in dark matter annihilation, might also heat up the IGM [86,94]. Nonetheless, we expect such signature to be distinguishable from the one of PBHs. A more quantitative evaluation of this issue will be presented elsewhere.
We assume that the dark ages end at z ∼ 30 (as it is standard convention), although in some scenarios star formation may start at earlier times and heat the IGM, hence changing both the sky-averaged and power spectrum of 21 cm IM (see e.g. [92,93]). In such cases, the uncertainties in the standard signal at z ∼ 30 would be larger and the identification of deviations as signatures of the presence of PBHs, more difficult.
We model the signature of massive PBHs, with abundances required to explain the current SMBH population in the 21 cm IM signal. We compute 2-point statistics of the fluctuations accounting explicitly for the temperature profiles around the PBHs in a comprehensive way, for the first time. We improve upon the work of [83,87,89,90] as we consider the scale-dependence of the PBH contribution to the spectrum (and not only a rescaling of the amplitude of the standard 21 cm IM signal or only the Poisson component).
After characterizing the PBH contribution to the stan-dard signal, we forecast the detectability with future experiments, ranging from the Square Kilometre Array (SKA, [95]) to a futuristic radio array on the dark side of the moon [96], which we refer to as the "Lunar Radio Array" (LRA). This paper is structured as follows. First, we review the standard 21 cm IM sky-averaged signal and power spectrum coming from the IGM in the dark ages, as well as the instrumental noise, in Sec. II. The effects of PBHs in the IGM and the spin temperature are characterized in Sec. III. Afterwards, the contribution to the 21 cm IM signal is modelled in Sec. IV and Sec. V for the sky-averaged signal and the power spectrum, respectively. Finally, forecasts for different future experiments are presented in Sec. VI. Discussions and conclusions can be found in Sec. VII. Throughout this paper, we assume the best fit values of the Planck 2015 TTTEEE+lowP power spectra [97] for the cosmological parameters.

II. STANDARD SIGNAL
We begin by reviewing the modelling of the standard 21 cm IM signal, i.e., without including the PBH contribution.

A. Sky-averaged signal
The optical depth of the IGM in the hyperfine transition is [98] where c is the speed of light, ν 0 = 1420.4 MHz is the rest-frame frequency of the hyperfine transition, A 10 = 2.85 × 10 −15 s −1 is the Einstein spontaneous emission rate coefficient for this transition, T s is the spin temperature of the gas, H(z) is the Hubble parameter, n H = 8.6 × 10 −6 Ω b h 2 (1 + z) 3 cm −3 is the hydrogen comoving number density [97], x H is the neutral fraction of hydrogen, k B is the Boltzmann constant and ∂ r v r is the comoving gradient of the peculiar velocity along the Line of Sight (LoS). We define T obs 21 as the observed differential brightness temperature between the 21 cm emission and the CMB: where δ b is the local baryon overdensity, h = H 0 /100 is the reduced Hubble constant and Ω m and Ω b are the mat-ter and baryon density parameters, respectively. Therefore, the sky-averaged 21 cm IM signal,T 21 , can be obtained from Eq. (2) by setting δ b = 0 and ∂ r v r = 0. We will mostly refer to the observed brightness temperature rather than to the local one, T loc 21 = T obs 21 (1 + z), throughout the paper, so we drop the superscript "obs" for simplicity.
Assuming that the background radiation includes only CMB photons, the spin temperature can be expressed as [99]: where T = 0.068K is the temperature correspondent to the 21 cm transitions, T k is the mean kinetic temperature of the IGM and y k and y α are the kinetic and Lyman-α coupling terms, respectively. We set T α ≈ T k , since it is a very good approximation when the medium is optically thick to Lyman-α photons [100], as in the case of study. The kinetic coupling term is due to the increase in the kinetic temperature by X-ray photon collisions with the gas: where C i are the de-excitation rates due to neutral hydrogen, electrons and protons, respectively . We use the fitting formulas of [101]: C e = n e γ e = n H (z)(1 − x H (z, r))γ e s −1 , where the number densities are in cm −3 and log(γ e /cm 3 /s) = −9.607 + 0.5 log T k exp(−(log T k ) 4.5 /1800) if T k ≤ 10 4 K, otherwise, γ e = γ e (T k = 10 4 ).
The coupling with the Lyman-α photons is described by the Wouthusyen-Field effect [99]. It depends on Lyman-α photons intensity,J 0 , given by: where ν α is the frequency of the Lyman-α transition, φ α is the fraction of the absorbed energy that goes into kinetic excitation of Lyman-α, N is the number of photons per unit area per unit time and σ is the absorption crosssection. We use the parametrization of [102], given by φ α = 0.48 1 − x 0.27 e 1.52 . Finally, the coupling term can be expressed as: where f 12 = 0.416 is the oscillator strength of the Lymanα oscillator.

B. Fluctuations
The optical depth and the spin temperature of a hydrogen cloud depend on its density and velocity divergence. Small anisotropies in these two quantities create fluctuations in T 21 . The 21 cm IM fluctuations power spectrum in the dark ages was computed in [103], and in [104] including the local velocity term. At the precision level we need in this work, given the uncertainties and assumptions in the modeling of the PBH contribution (see Sec. III, Sec. IV and Sec. V), it suffices to limit our computations to linear order. We follow the formalism developed in [105], which includes the effects due to supersonic relative velocities between baryons and dark matter [106]. This effect has been shown to help the formation of DCBHs at large redshifts [107], but it does not play a major role in the population of SMBHs at z ∼ 6 [108]. We refer the interested reader to [109,110] for a more detailed description of the 21 cm IM fluctuations, extending the formalism to higher order and including fluctuations in other quantities, such as the ionized fraction. Let . Then, at linear order, the fluctuations in the 21 cm IM signal can be expressed as: where α(z) = dT 21 /d δ b , including gas temperature fluctuations. The observed δT 21 in a directionn on the sky and at a certain frequency ν is given by where W ν (x) is the window function selecting the information at a certain frequency band centered in ν and x is the comoving distance along the LoS. This W ν (x) is a narrow function peaked at x(z) which depends on the experiment. Here we assume a Gaussian function of width ∆ν. In Fourier space, assuming that the baryons have caught up the dark matter and We can, therefore, define the transfer function of δT 21 as: where j is the spherical Bessel function with index , and we have defined J (kx) ≡ − ∂ 2 j (kx)/(∂kx) 2 , which can be written in terms of j , and j ±2 1 [104]. Given this, we can easily compute the 21 cm IM angular power spectrum at a certain frequency ν as: where P m (k) is the (isotropic) matter power spectrum. For computational efficiency, we will employ the flat-sky approximation [111] (for a pedagogical treatment, see e.g. [112,113]) for ≥ 10 3 .

C. Instrumental Noise
Although in the cosmic-variance limit the only source of noise is the variance arising by having a limited number of measurements of the power spectrum C , when considering an interferometer looking at the dark ages at a given frequency ν, there is an additional noise power spectrum [114][115][116][117]: where t o is the total time of observation, cover (ν) ≡ 2πD base /λ(ν) is the maximum multipole observable, D base being the largest baseline of the interferometer, f cover is the fraction of such baseline covered with antennas, and the amplitude T sys is the system temperature, which we assume to be the synchrotron temperature of the observed sky: found from extrapolating to lower frequencies the results of Ref. [118]. Therefore, the final uncertainty in the measurement of the C at the required multipole is: where f sky is the fraction of the sky observed by the experiment.

III. EFFECTS OF PBHS ON THE 21CM IM SIGNAL
The presence of PBHs affect the gas spin temperature: the PBH accretion triggers the emission of high-energy photons which heat and ionize the gas around the PBH. In this work, we present for the first time a computation of the 2-point statistics of the fluctuations accounting for the whole scale-dependence of the temperature profiles around the PBHs, focusing on linear perturbations in the dark ages.
An accreting PBH builds up a classical Bondi profile (i.e., r −3/2 ) around it. However, overdensities during the dark ages are still small and haloes are not formed yet. Therefore, as a first approach, we consider that there is no density profile in the gas around the PBH nor velocity inhomogeneities (δ b (r) = 0 and ∂ r v r = 0, respectively). Regarding the interaction between radiation and gas, we neglect radiative transfer effects (and limit ourselves to integrate over the frequency, as in Eq. (24)). Although these two effects might be relevant in some parameter configurations, they are competing: the former tends to reduce the volume affected by the PBH radiation, while radiative transfer increases the mean free path of high energy photons, hence increasing the distance to which X-rays can propagate and so the region heated by the PBH. While a more careful treatment will be needed, especially for comparison with observations and to assess their effective relative importance, here, for this initial exploration and signal-to-noise estimate, we assume that they compensate.
We assume that all processes are in equilibrium, given that their timescales are much smaller than the Hubble timescale. The steady-state approximation is very precise for masses M 3 × 10 4 M [119], but it breaks down for larger masses. Therefore, we limit our exploration to M ≤ 10 4 M . To explore a suitable mass range we consider three representative cases: M = 10 4 M , 10 3 M and 10 2 M . Given the slow growth of the PBHs at z 30, we assume that the PBH mass at different redshifts is the same when we perform a tomography analysis. Finally, we consider for simplicity that all PBHs have the same mass. This is an unrealistic scenario, but constraints for monochromatic mass distributions can be translated to any extended mass distribution using e.g., the methods proposed in [76,120].
We explain below the formalism we use to compute the temperature profiles around a PBH and show intermediate plots and results. Exact numerical calculations accounting for the time dependence can be found in [121].
A. Emission and neutral hydrogen fraction (xH (r)) IMBH emission is usually modelled by the combination of three components: a "multicolour disk black body spectrum" at low energies, a power-law spectrum from a surrounding "hot corona" at high energies and a small contribution from the reflected light from the corona by the gas around it. The contribution to the total emission due to the reflected radiation is small, but the light emitted by the disk produces a rather hard spectrum peaking at ∼ 1 KeV, as shown in e.g., [39,122] and references therein.
As the emission at low energies does not heat the gas around the PBH efficiently and sources at z > 22 contribute only little to the Near Infrared Background [39], we assume that gas accretion around the PBH powers only X-ray emission. Moreover we assume, as commonly done, that the emission is spherically symmetric. Therefore, a bubble with 21 cm IM signal different from the sky-averaged value is formed around the PBH. Finally, we can safely assume that PBHs of the masses we consider do not affect cosmic reionization [123].
Following [124], we assume that PBH accretion powers a miniquasar with a spherically symmetric power-law X-ray flux (limited to an energy range between 0.2 and 100 KeV). The difference between the heating of the gas by hard sources and those with a power-law spectrum may be significant (see e.g., [125,126]). However, we show in the Appendix A that the differences in the final angular power spectrum between a power-law spectrum and other more realistic choices (such as a piece-wise powerlaw [127] or including the emission from the disk as in [39]) are not significant with respect to the uncertainties in key parameters of the PBH population, i.e. their abundance, mass or Eddington ratio of the emission, as discussed below. Of course, in a refined application that goes beyond an initial feasibility analysis such as this paper, all these affects must be correctly modelled. Then, the spectrum of the photon emission, F (E), is given by: where A is a normalization factor chosen to have a luminosity L = λL Edd , where λ is the Eddington ratio and L Edd is the Eddington luminosity: Combining Eq. (18) and Eq. (19), it is easy to notice that λ and M are degenerate when computing the emission of the PBH, since A ∝ λM ≡ M. As explained below, relevant quantities, as x H or T 21 , only depend on the redshift and the intensity of the emission. Therefore, in order to illustrate how these quantities depend on both λ and M , we will show them in terms of M . The spectrum of Eq. (17) translates into number of photons per unit area per unit time at a comoving distance r from the source: where We use the fitting formula of [128] to compute the absorption cross section taking into account the contribution from helium and hydrogen atoms: with p = 2.65 if E < 0.25 keV and p = 3.30 if E > 0.25 keV. The emitted photons ionize the surrounding gas at a rate per hydrogen atom, Γ, as a function of the comoving distance r, given by: where x e (r)=(1 − x H (r)) is the ionized fraction, and the term E E0 φ(E, x e ) is introduced to take into account secondary ionizations. We apply the fitting formulas from [28] and [102] for E ≤ 0.5 KeV and E > 0.5 KeV, respectively.
Therefore, the neutral fraction is determined by the equilibrium between ionization and recombination rates: where α H = 2.6 × 10 −13 T −0.85 4 cm 3 /s is the recombination cross-section to the second excited atomic level, with T 4 = T k /10 4 K. For this computation, we assume T 4 = 1 (as in [124]).
The neutral fraction radial profile, x H (r), is shown in Fig. 1 for different redshifts and values of M. With increasing redshift, the hydrogen density increases; in a given volume at fixed photon flux, there are more atoms to ionize, hence the size of the ionized region decreases.
On the other hand, for increasing masses or Eddington ratios (i.e., larger M), as the PBH emission is more intense, the ionized region becomes larger.

B. Kinetic temperature
In addition to being ionized, the gas around the PBH is heated by the photons emitted by the miniquasar and cooled by the interaction with the CMB and the expansion of the Universe. The miniquasar heating affects the kinetic temperature, hence T k varies with the distance to the PBH. The heating rate per unit volume per unit time at a given comoving distance r from the source is: where f (x e (r)) is the fraction of the photon energy absorbed through collisional excitations. We use an extrapolation of the fitting formula of [102]: f = 0.9771(1−(1− x 0.2663 e ) 1.3163 ). As this fitting formula does not work well for a low-ionization medium (in reality f never goes to 0), we consider a floor f = 0.15 when x e ≤ 10 −4 [121].
Since the gas is exposed to Compton cooling by CMB photons, the heating rate per unit volume per unit time due to Compton processes is: where n e (z, r) = n H (z)x e (r) is the number density of electrons. On the other hand, the adiabatic expansion cooling per unit volume per unit time is ). Then, in equilibrium, H i = 0. Here, we do not consider Compton heating due to the emitted photons, because it is efficient only very close to the source [121]. Nonetheless, at those distances the hydrogen is totally ionized, so there is no signal in 21 cm IM and the results do not change. Moreover, those scales are far beyond the reach of 21 cm IM power spectrum resolution.
At large distances from the source, the gas is not affected by the PBH emission and its temperature is only determined by the adiabatic cooling due to the expansion of the Universe (there are no free electrons to scatter via Compton). Therefore, we need to set a contour condition by which T k (r → ∞) = T 0 k , the mean kinetic temperature of the IGM (without PBHs, which we take from the output of HyRec [129,130]). We include this condition in our computation of T k by a adding T 0 k x H (r) to the obtained T k . We will remove this contribution when computing T 21 of an isolated PBH.
We show gas temperature profiles as function of the comoving distance to the PBH in Fig. 2 for different redshifts and values of M. At large distances from the PBH, T k = T 0 k , hence the gas temperature is lower at lower redshifts. In the inner regions, the heating due to the emission of the PBH is coupled only to the neutral hydrogen, but, as the number of photons decays exponentially with the distance, this heating is more efficient close to the PBH. In these regions, PBH heating dominates over Compton and adiabatic cooling, so T k needs to be high to reach equilibrium. If Compton heating due to the emitted photons were considered, T k at distances tending to 0 would be much higher. However, as stated before, this would not change the signal in 21 cm IM because the hydrogen is totally ionized in those regions. At intermediate distances, PBH heating loses efficiency and T k drops even below T CMB until it reaches T 0 k . C. Spin temperature and differential brightness temperature Once we have computed the ionization fraction and gas temperature profiles (x H (r) and T k (r)), obtaining the spin and differential brightness temperature is straightforward using Eq. (3) and Eq. (2), respectively. T s may be driven whether by the collisional coupling or via the Wouthusyen-Field effect, whose weight is encoded in the coupling terms y k and y α in Eq. (3), respectively. We show radial profiles of y k and y α in Fig. 3 and Fig. 4, respectively, which make evident that T s is driven by collisional coupling in all the cases of study.
Spin temperature profiles can be seen in Fig. 5. T s behaves qualitatively similar to T k until T k ≈ T s < T CMB , where spin temperature coupling with CMB photons dominates and T s rises until T s ≈ T CMB , as can be seen in Fig. 3.
So far, we have applied the boundary condition that all quantities must match the standard values (i.e., without PBHs) when the distance to the PBH is large enough (e.g., T k (r → ∞) = T 0 k ). Nonetheless, we are interested on the isolated signal in 21 cm IM of a single PBH. Therefore, we subtract the contribution added due to the boundary condition in the same way that it was added before: where T 0 21 is the sky-averaged T 21 without PBHs. The T 21 (r) profile shown in Fig. 6 can be explained as follows. In the inner part, T 21 = 0 because all of the gas is ionized. The region with T 21 > 0 corresponds to the region where T k > T CMB and x H starts to grow; then, when T k drops because the PBH heating at those distances is less efficient, T 21 drops to negative values. Finally, T 21 rises again due to the collisional and Lymanα coupling of the photons to the source with the gas becomes totally inefficient and T s → T 0 s so T 21 → T 0 21 . Given that the PBH signal is isolated, at these distances, T 21 = 0.
In order to compute the fluctuations of T 21 , we need to compute also α profiles as a function of distance to the PBH, for which we follow the analytic expressions of [109]. Such profiles can be seen in Fig. 7.

IV. CONTRIBUTION TO THE SKY-AVERAGED SIGNAL OF 21 CM IM
Considerations about the minimum seed masses required [9,13], number of galaxies in the universe hosting SMBH [131], uncertainties on the accretion mechanisms, and CMB observations constraints on the maximum allowed dark matter fraction in PBH [59], lead us to consider a range of 10 −8 < Ω PBH < 10 −6 [123]. Key parameters of the model are largely unknown: SMBHs abundance and Eddington ratio (which is a proxy for the radiative efficiency) and mass. We consider here some representative values.
In addition to considering that all PBHs have the same mass, we also consider that all of them have the same Eddington ratio. This is an idealized case, since each kind of SMBH population (e.g. not active, type 1, type 2 and so on) has a different distribution of Eddington ratio (see e.g. [132][133][134]). It is customary to consider that SMBHs are active if λ 10 −4 , although this is an arbitrary limit, given that the Eddington ratio distribution is broad, and extends towards λ < 10 4 , as pointed by observations [135][136][137].
In any case, a characteristic value of the Eddington ratio is also largely unconstrained. Observational studies of X-ray selected SMBHs (which of course implies a selection bias favouring the most active luminous SMBHs) suggest large values of the Eddington ratio, i.e. λ ∼ 0.1. Nonetheless, one can consider that all SMBHs are active (not only those with λ 10 −4 ) and then, λ can take values 10 −4 [138]. Moreover, in [59], the evolution of PBH accretion under the most conservative assumptions was studied in a cosmological context assuming spherical accretion, finding much lower and mass dependent Eddington ratios. Besides, we assume for simplicity a duty cycle of unity, so the Eddington ratio would be smaller to match more realistic cases with lower duty cycles but higher luminosity.
Taking all this into account, we prefer to consider different parameter configurations to account for different possibilities spanning a wide range in the parameter space. We consider all the possible combinations of three masses (10 2 , 10 3 and 10 4 M ) and three abundances (Ω PBH = 10 −8 , 10 −7 and 10 −6 ). We also consider two possible scenarios with different choices of λ for each combination of M and Ω PBH : one with large Eddington ratio (λ ∼ 0.1 for astrophysical considerations) and another with small λ (see [59]). If a disk is formed and the accretion is not spherical, values of λ above this lower limit, but still below the astrophysical one, are expected [60]. Following, [59], as the change of λ with redshift for z 200 is small, we consider it constant and we take λ = 10 −4 for M = 10 4 M , λ = 10 −7 for M = 10 3 M and λ = 10 −10 for M = 10 2 M .
If we assume that there are PBHs present in the dark ages, their signal is superimposed to the standard one coming from the IGM and temperature fluctuations. We consider that the gas "bubble" around the PBH extends until the distance where |T 21 |< ∆T , which we set ∆T = 1 mK. This distance corresponds to the point in which T 21 (Fig. 6) becomes flat, and refer to it as r lim .
The differential flux per unit frequency received from the bubble can be expressed in terms of the differential brightness temperature as: where ∆Ω bubble = A/χ 2 (z), being A = πr 2 lim the comoving cross section of the bubble, and χ(z), the comoving distance to us. Furthermore, the line-integrated differential flux, δF , can be obtained multiplying the differential flux evaluated in the desired frequency, ν , by a redshift effective line width, ∆ν eff = (F(ν)dν ) /F(ν ). For an optically thin cloud, ∆ν eff can be approximated by: As both our gas and differential brightness temperature have radial profiles, we use an effective surface average defined as: The comoving number density of PBHs is: As was discussed in the previous section, λ and M are degenerate when considering the signal of an individual PBH. However, when considering the entire population, as the comoving number density of PBHs (Eq. (31)) only depends on Ω PBH and M , this degeneracy is broken. The average contribution of all the bubbles around the PBHs population to the differential flux per unit frequency is Finally, taking into account that ∆ν/∆z = ν 0 /(1+z) 2 and defining the beam-averaged effective differential brightness temperature, T 21 , using δF ν = 2 ν 2 rec k b T 21 ∆Ω beam /c 2 , we obtain [89]: We show the evolution of the sky-averaged signal with redshift in Fig. 8 for different cases with λ = 0.1. As can be seen, the contribution to the standard signal is positive (detected in emission) for z 130 (for which T 21 > 0 at any distance, Fig. 6), and negative (detected in absorption) for lower z. However, the contribution is only appreciable for z 50. For the same values of M, the PBH contribution is larger for larger n PBH , which is reasonable. On the other hand, for the same number density, the contribution is larger for larger intensity of the emission (i.e. larger M). Therefore, the contribution of PBHs to the standard sky-averaged signal in the cases where λ takes much smaller values will be negligible unless the number density is really high.

V. CONTRIBUTION TO THE ANGULAR POWER SPECTRUM OF 21 CM IM
In this section, we introduce how we compute the angular power spectrum of 21 cm IM, accounting for the first time for the emission of PBHs, the temperature profiles around them, and thus the full scale dependence of their contribution to the 21cm IM signal. As reference, note that the corresponding scale for the multipole at redshift z fulfills approximately kχ(z) = (using the Limber approximation); therefore, at z = 30, = 10 3 corresponds to k ∼ 0.09 Mpc −1 in a ΛCDM cosmology with the best fit parameters of Planck.
The modelling of the PBH signal in the fluctuations of T 21 is similar to that of the Sunyaev-Zel'dovich effect fluctuations from clusters of galaxies [139] or 21 cm IM from minihaloes before reionization [90]. In all these cases, there are extended sources tracing the peaks of the matter density field. In analogy, we use the halo model [140] to characterize the T 21 power spectrum during the dark ages in the presence of PBHs. A review of the formalism of the halo model can be found in [141]. Given that we only consider a monochromatic PBH population, all the integrals in mass that appear in the halo model formalism, which are of the type In the halo model, the power spectrum is the sum of two components: the correlation between points within the same halo or bubble is described by the 'one-halo' term, while the correlation between points in separate halos/bubbles is encoded in the 'two-halo' term. Hence P PBH (k) = P 1h PBH + P 2h PBH . In the same way, one can express the angular power spectrum in multipole coefficients as C PBH = C PBH(1h) + C PBH(2h) . We build on Eq. (11) to obtain the observed fluctuations of the 21 cm temperature fluctuations originated due only to the presence of PBHs in a directionn and in a frequency ν: where T 21 (r) and α(r) are the quantities obtained in Sec. III, r = √ x 2 + R 2 is the comoving distance to the center of the PBH and R = χ(z)/ is the comoving transverse distance to the center of the PBH. By using R = χ(z)/ , we assume a plane parallel approximation. This is justified because for low (where the plane parallel approximation breaks down), r r lim , hence δT 21,PBH (r) = 0. Once we have computed δT 21,PBH , we obtain the transfer function for the 21 cm IM fluctuations due to PBHs, T PBH , as in Eq. (12).
As the standard contribution in the linear regime without the PBHs comes from a continuum where there are no haloes, we consider that the one-halo term of the standard contribution vanishes. Therefore, we obtain the total angular power spectrum as the sum of the one-halo and two-halo terms, expressed as: where we assume that PBHs are completely correlated with the dark matter distribution and b is a scaleindependent bias. This is motivated by the following consideration. If PBHs are the seeds of the SMBHs, they are located at the centers of the potential wells so galaxies will form around them. We take the bias factor to be approximately the mean value of the galaxy bias. In explicit calculations we assume b = 1.25.
Given that the formation of a PBH is a rare event and PBHs spatial distribution is discrete, there is a Poissonian fluctuation in the number density of PBHs. Therefore, in addition to the standard matter power spectrum appearing in Eq. (36), there is a Poissonian power spectrum contribution. These fluctuations behave like isocurvature modes, as the formation of compact objects at small scales does not affect immediately the curvature at large scales [142]. The primordial power spectrum that describes them is: The isocurvature behaviour is enclosed in the transfer function of isocurvature modes, which is scaleindependent (T iso = 3 2 (1 + z eq ), where z eq is the redshift of matter-radiation equality, 1 + z eq ≈ 3400). Therefore, the power spectrum generated by the Poisson fluctuations is: where D(z) is the growth factor. The mass fraction, f PBH , appears because this contribution comes only from the fluctuation in number of PBHs and not all the matter. P Poisson should be added to the two-halo term multiplied only by T PBH . Nonetheless, given the ranges of f PBH we consider, the Poisson contribution is negligible at all scales. Only in studies exploring PBHs as a sizable fraction of the dark matter, where f PBH ∼ 1, it is found that the contribution of Eq. (38) dominates at small scales. In fact, Afshordi et al. (2003) [142] and Kashilinsky (2016) [143] propose to constrain the abundance of PBHs by looking for this scale independent contribution to the power spectrum in observations of the Lyα forest and the Cosmic Infrared Background anisotropies, respectively. Looking at Eq. (34), Eq. (35) and Eq. (36) it is easy to notice that the angular power spectrum will depend only on two quantities related to PBHs: n PBH and T PBH (C s also depend on other quantities not related with PBHs, such as the redshift). Therefore, although we do consider three parameters regarding PBHs, the relevant quantities are combinations of them: M = M λ and n PBH ∝ Ω PBH /M . The former is needed to compute the size of the bubble around the PBH (i.e., r lim ). Essentially, varying M shifts the features related with PBHs to different multipole ranges (via T PBH ). The latter is a rescaling of such contributions (C PBH(1h) ∝ n PBH and C PBH(2h) ∝ n 2 PBH ). Therefore, varying n PBH changes the   amplitude of the PBH features. These two effects are relevant to determine at which scale the PBH contribution starts to dominate. Thus, there is a degeneracy among the PBH parameters: where β is an arbitrary positive constant. All these effects can be seen in Fig. 9. In most of the cases, the PBH effects modify the standard power spectrum at ∼ 10 2 − 10 3 , with a large variation at ∼ 10 5 . As can be seen in Fig. 10, the PBH-induced deviation from the standard signal decreases with redshift because the size of the bubble decreases with redshift (see Fig. 6 and Fig. 7), so the multipole at which the deviation is appreciable at fixed n PBH and M increases.

VI. DETECTABILITY
We have characterized the imprints of massive PBHs in both the sky averaged signal and the power spectrum of 21 cm IM. The sky averaged signal requires dedicated single dipole experiments, such as EDGES [144], LEDA [145] or SARAS [146], to be measured. On the other hand, radio arrays as SKA aim to measure the fluctuations. As the PBH contribution on the sky-averaged signal is very small in most of the cases, we focus on the power spectrum and observations done with radio arrays.
We study the detectability of the signal and forecast constraints on massive PBH parameters assuming observations in the dark ages done with the SKA [95] and a futuristic Earth-based experiment, similar to the SKA but with much larger baseline and f cover , which we refer to as "SKA Adv ". Given that the atmosphere is opaque for frequencies 45 MHz, SKA will not be able to observe much further than z ≈ 30. Then, in order to observe well beyond the end of the dark ages (z 30) it will be necessary to observe from outside the Earth's atmosphere; a good candidate as a location for such observations is provided by the Moon [147,148]. This is why we also consider three different realizations of a futuristic radio array on the dark side of the Moon, that we call the "Lunar Radio Array" (LRA) [96]. The relevant specifications of the experiments considered can be found in Tab. I.

Spec
SKA SKA Adv LRA1 LRA2 LRA3  Because of its wide frequency coverage we consider that it will be possible to do tomography with LRA beyond z ∼ 30. We follow the arguments introduced in [80] to determine the redshift bins that can be considered independent when observing with ∆ν = 1 MHz between z = 30 and z = 200.
First of all, we compute the ∆ χ 2 considering SKA Adv as a function of Ω PBH and λ for two fiducial cases (M fid = 10 4 M , Ω fid PBH = 10 8 and λ fid = 0.1 (top) and λ fid = 10 −4 (bottom)). The 1σ, 2σ and 3σ contours are shown in Fig. 11 where the degeneracy among the parameters discussed in Sec. V can be appreciated.
We also forecast the errors on Ω PBH (σ Ω ) and λ (σ λ ) using Fisher matrices for all the fiducial cases we consider. The resulting forecasts are reported in Tab. II. The Fisher forecasts obtained should be considered as a rough estimate, especially for low fiducial values for λ: Fig. 11 shows that the constant ∆χ 2 contours are not well described by ellipses, which is what the Fisher approach assumes.
The PBH signal will be barely detected by SKA, since only for extreme cases in which n PBH is very large, the signal-to-noise ratio, S/N , for Ω PBH and λ is larger than unity.
As the amplitude of the power spectrum increases greatly at small scales, being able to resolve very small scales (i.e., large D base , which implies large cover ) will be key to detect the PBH signal and constrain the parameters. This is why the forecast uncertainties for SKA Adv are much smaller than for SKA (and similar considerations apply to LRA3 vs. LRA1).
On the other hand, the contribution of PBHs to the power spectrum decays with redshift (Fig. 10), hence the S/N between the case with PBHs and the standard one decreases fast with redshift, as shown in Fig. 12    information, cover has more impact in the final S/N . Generally, forecast errors for LRA1 are larger than for SKA Adv ; however, they are smaller for LRA2 than for SKA Adv , both having the same D base . This is true always except when both M and n PBH are large.
To summarize, although a detection of the PBH contribution in the dark ages might be achieved by SKA, in order to measure Ω PBH and λ accurately, a more ambitious experiment with a larger baseline is needed. Such measurements will be more precise if tomography is possible, for which experiments such as LRA are needed.

VII. DISCUSSION AND CONCLUSIONS
The origin and formation mechanism of SMBHs remains largely unknown. If the growth of the black holes happens only through (standard) accretion, in order to grow fast enough and reach M ∼ 10 9 M at z ∼ 7 [13] (and thus match the observed quasar abundance), massive seeds of ∼ 10 4 − 10 5 M need to be already present in regions with large gas densities at z ∼ 20. However, if mergers are also considered, the seeds can be lighter. Therefore, there are three candidates to be the seeds of SMBHs: remnants of Population III stars, DCBHs or intermediate mass PBHs.
In this work, we address the observational signatures that intermediate mass PBHs would have on 21 cm IM during the dark ages. We model this signal starting from the characterization of the radial profiles of T 21 around a single PBH to compute the contribution to the standard sky-averaged signal and to the angular power spectrum, using the halo model. This is the first time that the signature of PBHs accounting for its full scale dependence is modeled in the 21 cm IM power spectrum.
The values of the abundance of SMBHs (and therefore, of the seeds needed, Ω PBH ), the radiative efficiency (i.e., the Eddington ratio, λ) and the mass of the possible seeds, M , are largely unconstrained. Therefore, we consider several parameter configurations as fiducial cases. We forecast observational errors on λ and Ω PBH for each fiducial case assuming future observations made with SKA, a futuristic improved SKA-like experiment, and three different realizations of a futuristic radio array on the far side of the Moon (LRA).
We find that, although we consider three parameters (M , λ and Ω PBH ), the final power spectrum is only sensitive to two combinations of them: M λ = M (i.e., horizontal shifts of the one-and two-halo terms) and Ω PBH /M ∝ n PBH (i.e., changes in the amplitude of the PBH contribution to the power spectrum).
As a consequence, there is a degeneracy between the parameters, which can be expressed as C (M, λ, Ω PBH ) = C (M/β, λβ, Ω PBH /β) (with β being an arbitrary positive constant), as can be seen in Fig. 11. This perfect degeneracy is expected to be partially broken with more detailed modelling.
We find that the presence of PBHs increases the skyaveraged signal of 21 cm IM in absorption at z 50, but it is only appreciable when both M and n PBH are large (Fig. 8). With respect to the angular power spectrum, we find an enhancement of the signal for 10 2 − 10 3 (Fig. 9), which decays with redshift (because the size of the bubble around the PBH, i.e., the gas cloud producing a signal different from the standard sky-averaged value, is smaller for larger redshift, Fig. 6), as shown in Fig. 10. Although the enhancement is large, measuring λ and Ω PBH will be very difficult with SKA, as the effect is large only on small scales that can be reached only with a much longer baseline. On the other hand, as the signal-to-noise ratio decays fast with redshift, tomography does not add much information.
In this paper we have concentrated on the dark ages, given that they directly probe an epoch where the seeds should be present if they are primordial (and absent otherwise). Extending the analysis to lower redshift ranges would be interesting for experiments happening on a shorter timescale, although added complications due to astrophysics and degeneracy with other signals would be involved.
Our modelling makes several assumptions and simplifications, which we recap and discuss their resulting implications here. First of all, we consider that there is no overdensity surrounding the PBH (or that the profile around it does not affect drastically the signal) and also neglect radiative transfer effects. For this initial exploration of the subject, we assume that the effects of these two assumptions compensate, since accounting for density profiles would generate smaller bubbles but larger mean free paths of X-rays (consequence of the radiative transfer) would make the bubbles larger.
To compute the contribution of PBHs to the standard signal, we model a single PBH and afterwards we use the number density of PBHs, n PBH ∝ Ω PBH /M , to account for the full population. Hence, our formalism breaks down for large number densities. In such scenario, bubbles around different PBHs overlap, so PBHs can not be considered as isolated anymore. Besides, PBHs contribute significantly to cosmic reionization, advancing it if their number density is too high. In this case, more accurate modelling is needed. Actually, the cases with the largest n PBH considered here should be interpreted carefully due to the effects commented above. This caveat could also be more relevant if the actual bubbles are larger than considered here due to radiative transfer effects.
Moreover, we assume a simple power-law spectrum for the radiation emitted by the PBH accretion without modelling the full spectral energy distribution (although see Appendix A) and that all processes are in equilibrium (steady-state approximation, hence we are limited to M ≤ 10 4 M ). Although the effect of supersonic relative velocities between baryons and dark matter is included in our computation of the 21 cm IM fluctuations (see Sec. II B), it is not included in the modelling the heating of the IGM due to the PBH emission. Relative streaming velocities between gas and PBHs leave an imprint on T 21 radial profiles at the corresponding baryon acoustic oscillation scales, imprinting the corresponding features in the total angular power spectrum. The interested reader can find a study of the effects of relative velocities in the 21 cm IM power spectrum in the prereionization era, but after the first stars formed (hence at lower redshifts that those we are focused on) in e.g., [149]. We expect a similar qualitative behaviour for the case of PBHs at larger redshifts.
Finally, we have assumed an average value for the bias between the seeds and the dark matter distribution, while in reality its value might change with redshift and the mass of the seeds. However, given that the value of the bias is strictly related to the height of the peaks in the density field, it is also directly connected to the PBH initial mass and number. In principle, given that it affects the two-halo term contribution but not the one-halo term, variations of the bias would cause a slightly different signal, but we do not expect the final result of this paper to be substantially different.
Nonetheless, the impact of these assumptions and simplifications on the final power spectrum in the scenario under study are subdominant, given the magnitude of the uncertainties due to the PBH parameters. On top of this, we assume that a comprehensive characterization of the foregrounds which affect the detectability of the signal is possible, hence they do not affect the S/N or the forecast uncertainties in the PBH parameters reported in Sec. VI. In addition, we have considered that there is no other exotic energy injection during the dark ages and that star formation begins at z 30, although this might not be the case. In such cases, the identification of a signal as the product of the IGM heating due to the PBHs would be more difficult. A comprehensive study of these effects using simulations and radiative transfer codes to account for PBH distribution, clustering, relative velocities, gas accretion, mergers, and/or extended mass distributions of the PBHs as well as an estimation of how removing the foreground wedge or an early star formation affects the detectability is left for future work.
There are previous proposals to identify the seeds of the SMBHs from their observational signatures, e.g., with the 21 cm IM sky-averaged signals at 10 z 30 [91] to distinguish between black holes formed from remnants of Population III stars or DCBHs (although the signal at larger redshifts might also come from PBHs) or with spectral distortions, to ascertain if the seeds are primordial [46]. DCBHs are also one of the preferred candidates to explain the power spectrum of the Near Infrared Background and its cross correlation with the cosmic X-ray radiation, both at large scales [39]. The emerging spectrum from the DCBH environment is non-zero only in this window [123], which may be useful for identifying them with Chandra 2 or Athena 3 . Moreover, while the growth of remnants of Population III stars may remain undetectable for JWST 4 , the evolved stages of DCBHs might be identifiable [150]. However, these signatures might also have been produced by massive PBHs. We leave the study of this scenario for future work. With advanced gravitational wave detectors, such as LISA 5 , it will be possible to measure gravitational waves created in mergers of SMBHs at large redshifts, which will offer insights on the environments and history of such black holes and help to discriminate among the different candidates for being the seeds.
The advent of new experiments and corresponding observations will shed light on how SMBHs reached such huge masses and on the nature of the massive seeds needed to explain their existence. It is also possible that the three kinds of seeds discussed above coexist and give different signatures. We eagerly await observations that will open the window toward higher redshifts and will give us the opportunity to improve our understanding of some of the most extreme structures in the Universe. In this appendix we discuss the dependence of the final T 21 radial profile and power spectrum on the emitted spectrum assumed. We referred to the spectrum used in the main text (Eq. (17)) as Power-Law (PL) in opposition to a Power-Law with Low Energies (PL LE) in which the energy range is extended at the low energy limit (10.4 eV ≤ E ≤ 100 keV) and a Power-Law with High Energies (PL HE) in which the energy range is extended at the high energy limit (200 eV ≤ E ≤ 300 keV). All of them have the same exponent: −1.
We also consider a more elaborated spectrum, as the one introduced by Sazonov, Ostriker & Sunyaev (2004) [127]: 10.4 eV < E < 1 keV, E −1 , 1 keV < E < 100 keV, E −1.6 , 100 keV < E, (A1) with A(M) computed as in Eq. (18). We consider one case with a high energy cut of 100 keV (SOS LE) and another with the cut at 300 keV (SOS HE). Finally, we con-sider a more realistic spectrum which includes the contribution of the disk as a multicolor black body spectrum, added to a power-law with index −1 for energies larger than ∼ 3T max (where k B T max = (M/M ) −0.25 KeV) and with a high energy exponential cut off at 300 KeV, modelling the emission of the hot corona. We follow [39] and normalize each contribution to the total emission to have the same luminosity. The emission from the disk can be  We show the resulting T 21 (r) and C for all the emission models explained above in Fig. 13 and Fig. 14, respectively. Although the radial profiles of T 21 are different, the effect on the final power spectrum is small, compared with the uncertainties in the PBH parameters (i.e., Ω PBH , M and n PBH ). Moreover, as the dependence on the PBH parameters is the same for all the different emission spectra, significant changes on the forecasts reported on Tab. II or on the two dimensional confidence levels shown in Fig. 11 are not expected.