Constraining the primordial black hole abundance with 21cm cosmology

TL;DR

This paper explores how 21cm cosmology observations during reionization and cosmic dawn can tightly constrain the abundance of primordial black holes, which could significantly influence early universe properties and dark matter composition.

Contribution

It introduces a novel approach to constrain solar mass primordial black holes using 21cm observations, considering multiple effects on the intergalactic medium.

Findings

21cm observations can improve constraints on PBHs by over an order of magnitude.

PBHs cause distinct effects on the intergalactic medium, detectable via 21cm signals.

Upcoming experiments like SKA and HERA are key to these constraints.

Abstract

The discoveries of a number of binary black hole mergers by LIGO and VIRGO has reinvigorated the interest that primordial black holes (PBHs) of tens of solar masses could contribute non-negligibly to the dark matter energy density. Should even a small population of PBHs with masses $\gtrsim \mathcal{O}(M_\odot)$ exist, they could profoundly impact the properties of the intergalactic medium and provide insight into novel processes at work in the early Universe. We demonstrate here that observations of the 21cm transition in neutral hydrogen during the epochs of reionization and cosmic dawn will likely provide one of the most stringent tests of solar mass PBHs. In the context of 21cm cosmology, PBHs give rise to three distinct observable effects: ${\textit{(i)}}$ the modification to the primordial power spectrum (and thus also the halo mass function) induced by Poisson noise, ${\textit{(ii)}}$ a uniform heating and ionization of the intergalactic medium via X-rays produced during accretion, and ${\textit{(iii)}}$ a local modification to the temperature and density of the ambient medium surrounding isolated PBHs. Using a four-parameter astrophysical model, we show that experiments like SKA and HERA could potentially improve upon existing constraints derived using observations of the cosmic microwave background by more than one order of magnitude.

Tables

Table 1. Table 1: Fit coefficients characterizing the radiative efficiency of ADAF accretion, described by Eq. ( 8 ), in terms of the net accretion rate. From Ref.  Xie and Yuan ( 2012 ) , for δ = 0.1 𝛿 0.1 \delta=0.1 and for the viscous parameter α = 0.1 𝛼 0.1 \alpha=0.1 .

${\dot{M}}_{\mathrm{PBH}}/{\dot{M}}_{\mathrm{Edd}}$	${ϵ}_0$	$a$
$\lesssim 9.4\times {10}^{-5}$	0.12	0.59
$9.4\times {10}^{-5}-5\times {10}^{-3}$	0.026	0.27
$5\times {10}^{-3}-6.6\times {10}^{-3}$	0.50	4.53

Table 2. Table 2: Model parameters and priors adopted in the sensitivity analysis.

	Minimum value	Maximum value	Prior type
$f_{\mathrm{PBH}}$	${10}^{-8}$	1	Log
${\zeta }_{\mathrm{UV}}$	10	100	Log
${\zeta }_X$	$2\times {10}^{55}{M}_{\odot }^{-1}$	$2\times {10}^{57}{M}_{\odot }^{-1}$	Log
$T_{\min }$	${10}^4$ K	${10}^5$ K	Log
$N_{\alpha }$	$4\times {10}^2$	$4\times {10}^4$	Log

Equations

$$\frac{n_1}{n_0}=3e^{-T_0/T_S},$$

$$T_S=\frac{T_{\mathrm{R}}+y_k{T}_k+y_{\alpha }{T}_{\alpha }}{1+y_k+y_{\alpha }},$$

$$y_k=\frac{T_0}{A_{10}{T}_k}{n}_{\mathrm{H}}\left(x_{\mathrm{H}}{C}_{\mathrm{HH}}+(1-x_{\mathrm{H}})C_{\mathrm{eH}}+(1-x_{\mathrm{H}})C_{\mathrm{pH}}\right),$$

$$y_{\alpha }=\frac{16{\pi }^2{e}^2{f}_{12}}{27m_{e}c{A}_{10}}\frac{T_0}{T_k}{J}_{\alpha },$$

$$\delta T_b(\nu )=\frac{T_S-T_{\mathrm{R}}}{1+z}(1-e^{-{\tau }_{{\nu }_0}}),$$

$${\left(\frac{dE}{dVdt}\right)}_{\mathrm{inj}}=L_{\mathrm{acc}}{n}_{\mathrm{PBH}}=L_{\mathrm{acc}}\frac{f_{\mathrm{PBH}}{\rho }_{\mathrm{DM}}}{M_{\mathrm{PBH}}},$$

$$L_{\mathrm{acc}}=ϵ{\dot{M}}_{\mathrm{PBH}}.$$

$$ϵ={ϵ}_0{\left(\frac{100{\dot{M}}_{\mathrm{PBH}}}{{\dot{M}}_{\mathrm{Edd}}}\right)}^a,$$

$${\dot{M}}_{\mathrm{Bondi}-\mathrm{Hoyle}}=4\pi \lambda {\rho }_{\infty }\frac{(G{M}_{\mathrm{PBH}}{)}^2}{c_{s,\infty }^3},$$

$$\sqrt{⟨v_L^2⟩}\simeq \min \left[1,\frac{1+z}{1000}\right]\times 30km/s.$$

$$v_{\mathrm{eff}}\equiv {⟨\frac{1}{{\left(c_{s,\infty }^2+v_{\mathrm{rel}}^2\right)}^{3(1+a)/2}}⟩}^{-\frac{1}{3(1+a)}}=c_{s,\infty }{\left[{\left(\frac{3}{2}\right)}^{\frac{3}{2}}U\left(\frac{3}{2},1-\frac{3a}{2},\frac{3}{2}{\mathcal{M}}^{-2}\right){\mathcal{M}}^{-3}\right]}^{-\frac{1}{3(1+a)}}$$

$$v_{\mathrm{eff}}\simeq \{\begin{array}{ll}{c}_{s,\infty }{\mathcal{M}}^{1/(1+a)}{\left[3\sqrt{\frac{3}{2\pi }}B\left(\frac{3a}{2},\frac{3}{2}\right)\right]}^{-\frac{1}{3(1+a)}}& ,\mathcal{M}\gg 1\\ c_{s,\infty }& ,\mathcal{M}\ll 1,\end{array}$$

$${\dot{M}}_{\mathrm{PBH}}=4\pi \lambda {\rho }_{\infty }\frac{(G{M}_{\mathrm{PBH}}{)}^2}{v_{\mathrm{eff}}^3}.$$

$${\left(\frac{d{E}_c}{dVdt}\right)}_{\mathrm{dep}}=f_c(z){\left(\frac{dE}{dVdt}\right)}_{\mathrm{inj}}\equiv n_b{ϵ}_c^{\mathrm{PBH}},$$

$$f_c(z)=\frac{H(z)\int \frac{d\ln (1+z^{\prime })}{H(z^{\prime })}\int d\omega T_c(z,z^{\prime },\omega )L_{\mathrm{acc}}(z^{\prime },\omega )}{\int d\omega L_{\mathrm{acc}}(z,\omega )},$$

$$L(\omega )\propto \Theta (\omega -{\omega }_{\min }){\omega }^{\beta }\exp (-\omega /{\omega }_s),$$

$$\frac{d{x}_e}{dz}=\frac{dt}{dz}\left({\Lambda }_{\mathrm{ion}}-{\alpha }_{\mathrm{A}}C{x}_e^2{n}_b{\mathfrak{f}}_{\mathrm{H}}\right),$$

$${\Lambda }_{\mathrm{ion}}^{\mathrm{PBH}}={\mathfrak{f}}_{\mathrm{H}}\frac{{ϵ}_{\mathrm{HI}}^{\mathrm{PBH}}}{E_{\mathrm{HI}}}+{\mathfrak{f}}_{\mathrm{He}}\frac{{ϵ}_{\mathrm{HeI}}^{\mathrm{PBH}}}{E_{\mathrm{HeI}}},$$

$$\frac{d{T}_k}{dz}=\frac{2T_k}{3n_b}\frac{d{n}_b}{dz}-\frac{T_k}{1+x_e}\frac{d{x}_e}{dz}+\frac{2}{3k_B(1+x_e)}\frac{dt}{dz}\sum _j{ϵ}_{\mathrm{heat},j},$$

$$J_{\alpha ,\mathrm{PBH}}=\frac{c{n}_b}{4\pi }\frac{{ϵ}_{\mathrm{Ly}\alpha }^{\mathrm{PBH}}}{h{\nu }_{\alpha }}\frac{1}{H(z){\nu }_{\alpha }},$$

$$\delta T_b(\nu )\simeq 27x_{\text{H}}(1+{\delta }_b)\left(1-\frac{T_{\text{CMB}}}{T_S}\right)\left(\frac{1}{1+H^{-1}\partial v_r/\partial r}\right){\left(\frac{1+z}{10}\right)}^{1/2}{\left(\frac{0.15}{{\Omega }_{\mathrm{m}}{h}^2}\right)}^{1/2}\left(\frac{{\Omega }_{\mathrm{b}}{h}^2}{0.023}\right)\text{mK},$$

$$⟨{\tilde{\delta }}_{21}(\mathbf{k},z){\tilde{\delta }}_{21}^{\ast }({\mathbf{k}}^{\prime },z)⟩\equiv (2\pi {)}^3{\delta }^D(\mathbf{k}-{\mathbf{k}}^{\prime })P_{21}(k,z),$$

$$4\pi r^2\rho |v|$$

$$v\frac{dv}{dr}$$

$${\Lambda }_{\mathrm{ion}}(1-x_e)$$

$$r_{\mathrm{B}}=\frac{G{M}_{\mathrm{PBH}}}{v_{\mathrm{eff}}^2}.$$

$$F(E,r,z)=\frac{L_{\mathrm{acc}}(E,z)}{4\pi r^2}{e}^{-{\tau }_{\mathrm{ion}}(E,r,z)},$$

$${\tau }_{\mathrm{ion}}(E,r,z)={\int }_0^r{n}_{\mathrm{H}}(z)\sigma (E)x_{\mathrm{HI}}(r,z),$$

$${\Lambda }_{\mathrm{ion}}={\int }_{E_0}^{\infty }\frac{dE}{E}\sigma (E)F(E,r,z)\left[1+\frac{E-E_{\mathrm{th}}}{E_{\mathrm{th}}}{f}_{\mathrm{ion}}(E,x_e)\right],$$

$$\frac{d{T}_{b,\nu }(r)}{dr}=-{\kappa }_{\nu }(r)\left[T_{b,\nu }(r)-T_S(r)\right],$$

$${\kappa }_{\nu }(r)=\frac{h{\nu }_0}{4\pi }{\varphi }_{\nu }\left[B_{01}{n}_0(r)-B_{10}{n}_1(r)\right]\simeq \frac{3c^2{A}_{10}}{32\pi {\nu }_0^2}{x}_{\mathrm{H}}(r)n_{\mathrm{H}}{\varphi }_{\nu }\frac{T_0}{T_S(r)},$$

Full text

††preprint: IFIC/18-15

I Introduction

The canonical cosmological model assumes that dark matter is comprised of a cold gas of weakly interacting particles. Despite its simplicity, this minimal scenario provides an excellent fit to both cosmic microwave background (CMB) and large-scale structure measurements Adam et al. (2016); Ade et al. (2016); Aghanim et al. (2016); Alam et al. (2017). However, a precise understanding of the fundamental nature of dark matter is missing and remains at the forefront in the current list of unsolved problems in modern physics.

Although dark matter is usually interpreted in terms of a new elementary particle, other alternatives exist. Black holes (BHs) produced prior to big bang nucleosynthesis (BBN), i.e., primordial black holes (PBHs), represent such an alternative — remarkably, this solution is as old as particle dark matter Chapline (1975). In particular, this scenario has recently attracted much attention Bird et al. (2016); Clesse and García-Bellido (2017a, b); Sasaki et al. (2016); Raidal et al. (2017); Ali-Haïmoud et al. (2017); Kavanagh et al. (2018) in the context of the LIGO and VIRGO discoveries of several binary BH mergers Abbott et al. (2016a, b, 2017a, c, 2017b).

The idea that PBHs could be formed by strong accretion during the radiation-dominated epoch was first introduced five decades ago Zel’dovich and Novikov (1967). It was later suggested that initial fluctuations in the early Universe could give rise to a large number of gravitationally collapsed objects with masses above the Planck mass Hawking (1971), that however, would not grow substantially by accretion Carr and Hawking (1974). Very light PBHs would not survive until the present epoch, though; the prediction that any BH should emit particles with a black body spectrum Hawking (1974, 1975) implies that PBHs with masses below $10^{-18}M_{\odot}$ (with $M_{\odot}$ the Sun’s mass) would have evaporated on cosmological times scales Hawking (1974, 1975); Page (1976). While the evaporation of such light PBHs would have interesting observational consequences Khlopov (2010); Carr et al. (2010); Calmet et al. (2014), only heavier PBHs could constitute a non-negligible fraction of the dark matter content of the Universe Chapline (1975) (see, e.g., Refs Belotsky et al. (2014); Carr et al. (2016); Sasaki et al. (2018) for recent reviews). Alternative mechanisms for the formation of PBHs have also been proposed, such as the collapse of domain walls Sato et al. (1981); Maeda et al. (1982); Berezin et al. (1983) or cosmic strings Hogan (1984); Hawking (1989); Polnarev and Zembowicz (1991), first-order phase transitions Crawford and Schramm (1982); Hawking et al. (1982); Kodama et al. (1982); Hall and Hsu (1990); Moss (1994); Konoplich et al. (1999); Jedamzik and Niemeyer (1999); Khlopov et al. (2000), or fragmentation of scalar condensates Cotner and Kusenko (2017a, b); Cotner et al. (2018). Today, the most studied mechanism of PBH formation is gravitational collapse Carr (1975); Ivanov et al. (1994); Yokoyama (1997); Shibata and Sasaki (1999); Young et al. (2014) of the fluctuations predicted in inflationary scenarios Carr and Lidsey (1993); Carr et al. (1994); Ivanov et al. (1994); Yokoyama (1997); García-Bellido et al. (1996); Taruya (1999); Green and Malik (2001); Bassett and Tsujikawa (2001); Pi et al. (2018).

Following the recent observations of BH mergers by the LIGO and VIRGO experiments Abbott et al. (2016a, b, 2017a, c, 2017b), PBHs of several solar masses or heavier have been extensively examined and constrained by a number of different means (see, e.g., Refs. Mack et al. (2007); Ricotti et al. (2008); Josan et al. (2009); Carr et al. (2010); Capela et al. (2013); Clesse and García-Bellido (2017a); Carr et al. (2016); Green (2016); Bellomo et al. (2018); Kühnel and Freese (2017); Carr et al. (2017); Sasaki et al. (2018) for compilations of constraints, also at other mass scales): gravitational lensing Green (2017); García-Bellido and Clesse (2018); Zumalacarregui and Seljak (2018); García-Bellido et al. (2018), dynamical constraints Monroy-Rodríguez and Allen (2014); Brandt (2016); Green (2016); Li et al. (2017); Koushiappas and Loeb (2017); Kavanagh et al. (2018), radio and X-rays measurements Gaggero et al. (2017); Inoue and Kusenko (2017); Hektor et al. (2018a); Manshanden et al. (2019), and spectral and anisotropy distortions of the CMB Tada and Yokoyama (2015); Young and Byrnes (2015); Chen et al. (2016); Ali-Haïmoud and Kamionkowski (2017); Horowitz (2016); Aloni et al. (2017); Poulin et al. (2017); Bernal et al. (2017); Nakama et al. (2018); Deng et al. (2018). If energy is injected into the intergalactic medium (IGM) via the accretion of matter onto PBHs, the ionization and thermal evolution of the Universe could be affected and thus, both the spectral shape Carr and Lidsey (1993); Carr et al. (1994); Chluba et al. (2012); Kohri et al. (2014) and the anisotropies of the CMB Carr (1981); Ricotti et al. (2008, 2008); Chen et al. (2016); Ali-Haïmoud and Kamionkowski (2017); Aloni et al. (2017); Bernal et al. (2017); Horowitz (2016); Poulin et al. (2017). Namely, the work of Ref. Poulin et al. (2017) excludes PBH with masses $\gtrsim 2\,M_{\odot}$ as the sole component of dark matter, although uncertainties associated with the accretion make the determination of a robust upper limit difficult. Moreover, Poisson-noise fluctuations could become too large to be consistent with the power spectrum on scales of the Lyman-$\alpha$ forest, setting an upper limit on PBH masses near a few times $10^{4}\,M_{\odot}$ Afshordi et al. (2003), which has been recently improved to $M_{\rm PBH}\lesssim 60\,M_{\odot}$ Murgia et al. (2019).

The presence of PBHs could also change the process known as reionization, during which the first sources in the Universe emitted ultraviolet photons that ionized the neutral hydrogen. Collectively, observations from the CMB observations, Lyman-$\alpha$ emission in star-forming galaxies, and of the Gunn-Peterson optical depth from bright quasars indicate that the Universe should be fully ionized by $z\sim 6$. Nevertheless, only observations of the CMB are currently sensitive to the early part of reionization, and to the epoch during which the first stars formed (referred to as the cosmic dawn). The CMB, however, is only sensitive to a cumulative line-of-sight effect (observed via modifications to the integrated optical depth), and thus cannot provide detailed information about the state of the IGM during these epochs. The only proposed method for understanding the nature of the Universe on these time-scales is 21cm cosmology, which attempts to observe the hyperfine transition between the singlet and triplet levels of neutral hydrogen. To date, only one experiment working has claimed a putative signal; using a single dish radio telescope, the EDGES (Experiment to Detect the Reionization Step) collaboration Bowman et al. (2018) claimed detection of a large absorption dip in the globally averaged signal, appearing around $z\sim 17$, but whose amplitude cannot be explained within the standard $\Lambda$CDM scenario. Determining the validity of this signal, however, is notoriously challenging, as foregrounds at these frequencies are expected to be at least $\sim 10^{4}$ times larger than the signal of interest. Thus, the 21cm contribution can only be discriminated by systematically removing spectrally smooth foreground components from the observation. Next-generation radio interferometers such as HERA (Hydrogen Epoch of Reionization Array) Beardsley et al. (2015) and SKA (Square Kilometer Array) Mellema et al. (2013) will also attempt to observe these epochs. A major advantage of using an interferometer in these observations is that one also has direct information on the scales over which foregrounds emanate, allowing for a more robust removal of foregrounds and providing additional spatial information on the 21cm signal.

The field of 21cm cosmology should be considered particularly important for those in the field of particle physics, as the wealth of available information might allow physicists to probe any new process which, e.g., modifies structure formation, the energy in the IGM, star formation and evolution, and photon propagation. Studies have been performed in the context of various dark matter candidates, for example, illustrating unprecedented sensitivity to a wide array of dark matter candidates Valdés et al. (2013); Sitwell et al. (2014); Evoli et al. (2014); Carucci et al. (2015); Lopez-Honorez et al. (2016); Berlin et al. (2018); D’Amico et al. (2018); Liu and Slatyer (2018); Mitridate and Podo (2018); Pospelov et al. (2018); Aristizabal Sierra and Fong (2018); Boddy et al. (2018); Kovetz et al. (2018); Barkana et al. (2018); Barkana (2018); Fialkov et al. (2018); Escudero et al. (2018); Schneider (2018); Safarzadeh et al. (2018); Lidz and Hui (2018); Lopez-Honorez et al. (2019). Existing works looking into the effects of PBHs in 21cm cosmology have independently focused on understanding the effects arising from PBH accretion on the globally averaged signal Hektor et al. (2018b), Poisson-noise on the power spectrum Gong and Kitajima (2017, 2018), and local modifications to the medium around PBHs on the globally averaged signal and the angular power spectrum Tashiro and Sugiyama (2013); Bernal et al. (2018). In this work, we compare and contrast these three effects in the context of $\gtrsim\mathcal{O}(M_{\odot})$ PBHs, illustrating that the only process relevant for $M_{\rm PBH}\lesssim 10^{3}\,M_{\odot}$ is the heating and ionization of the IGM arising from X-rays produced during accretion. Using the 21cm power spectrum, which offers a far more robust and powerful probe the globally averaged signature, we illustrate the near-future sensitivity of both HERA and SKA.

The structure of the paper is as follows. In Section II we review the basics of 21cm cosmology, in particular focusing on the computation of the 21cm brightness temperature and power spectrum. In Section III we discuss the impact that PBHs have on the observable 21cm signal; in particular we study global heating and ionization as a consequence of accretion around PBHs (Section III.1), and how the process of accretion heats and ionizes the local medium surrounding the PBH (Section III.2), and we also address how Poisson-noise induced in the power spectrum by the discrete nature of PBHs modifies the signal arising in minihalos (Section III.3). We comment on the sensitivity analysis in Section IV, and then present the results and conclude in Section V.

II 21cm cosmology

Below, we present a brief overview of the fundamentals of 21cm cosmology. A more extensive overview can be found, e.g., in Refs. Madau et al. (1997); Furlanetto et al. (2006); Pritchard and Loeb (2012); Furlanetto (2015).

The intensity of the redshifted 21cm line provides direct information about the fraction of neutral hydrogen residing in excited and ground states. This is conventionally characterized by the so-called spin temperature $T_{S}$, and is formally defined as

[TABLE]

where $n_{0}$ and $n_{1}$ are the number density of neutral hydrogen atoms in the ground (singlet) and excited (triplet) states, $k_{B}T_{0}=h\nu_{0}$ is the energy of the 21cm photon, and the factor of 3 comes from the degeneracy of the triplet excited state. The spin temperature of the 21cm line is controlled by three processes: (i) absorption of and stimulated emission induced by a background radiation field, (ii) collisions of neutral hydrogen with hydrogen atoms, free protons, or free electrons, and (iii) indirect excitations/de-excitations induced by scattering with ambient Lyman-$\alpha$ photons. Consequently, the spin temperature can be expressed directly in terms of temperatures characterizing the efficiency of each of these processes,

[TABLE]

where $T_{\rm R}$, $T_{k}$, and $T_{\alpha}$ are the temperature of the background radiation, kinetic temperature of the gas, and the temperature associated with the intensity of the Lyman-$\alpha$ emission (for all environments of interest, it is sufficient to take $T_{\alpha}\simeq T_{k}$), and $y_{k}$ and $y_{\alpha}$ are their respective couplings.

The collisional coupling, $y_{k}$, is given by

[TABLE]

where $A_{10}=2.85\times 10^{-15}{\rm s}^{-1}$ is the Einstein coefficient of spontaneous emission rate, $n_{\rm H}$ is the hydrogen comoving number density, $x_{\rm H}$ is the neutral fraction of hydrogen, and $C_{\rm HH}$, $C_{\rm eH}$ and $C_{\rm pH}$ are the de-excitation rates due to collisions between a hydrogen atom and another hydrogen atom, an electron, or a proton, respectively. These quantities are computed using the fitting formulas of Ref. Kuhlen et al. (2006).

The Lyman-$\alpha$ coupling, $y_{\alpha}$, can be written in terms of the Lyman-$\alpha$ flux $J_{\alpha}$,

[TABLE]

where $e$ and $m_{e}$ are the charge and mass of the electron, and $f_{12}=0.416$ is the oscillator strength of the Lyman-$\alpha$ transition.

Rather than working directly with observed intensity, it is conventional for radio astronomers to use the effective brightness temperature, $T_{b}(\nu)$, which is proportional to the specific intensity (or spectral radiance) of a blackbody in the Rayleigh-Jeans limit, $\nu\ll T_{b}(\nu)$. The change in intensity relative to the background radiation field is thus given by the differential brightness temperature:

[TABLE]

where $\tau_{\nu_{0}}$ is the optical depth of the 21cm line. Throughout this work we will implicitly assume the background radiation field to simply be the CMB, $T_{\rm R}=T_{\rm CMB}$. We note, however, that in the light of the abnormally deep absorption trough observed by EDGES Bowman et al. (2018), several scenarios predicting an unexpected excess of radio background have been explored Feng and Holder (2018); Ewall-Wice et al. (2018); Fraser et al. (2018); Mirocha and Furlanetto (2019); Pospelov et al. (2018); Dowell and Taylor (2018); Jana et al. (2019); Fialkov and Barkana (2019); Ewall-Wice et al. (2019) (see, however, Ref. Sharma (2018)).

III The effects of primordial black holes

The existence of a population of $M_{\rm PBH}\gtrsim M_{\odot}$ PBHs has the potential to modify the 21cm signal in a number of different ways. In order to understand exactly how and where these effects would arise, let us consider for the moment a single isolated PBH. Matter surrounding the PBH is accreted at some rate $\dot{M}_{\rm PBH}$. The in-falling matter heats and ionizes in a localized region surrounding the PBH. Some of the energy gained via accretion is radiated by the newly formed plasma, primarily in the form X-rays produced by thermal synchrotron and bremsstrahlung emission, and modified by comptonization. Additionally, energy resulting from the in-fall of matter can be viscously dissipated, directly heating the accreting matter rather than radiating X-rays. This type of energy dissipation is at the heart of so-called advection-dominated accretion flow (ADAF) (see, e.g., Ref. Yuan and Narayan (2014) for a review). It has recently been argued Poulin et al. (2017) that accretion onto $\sim\mathcal{O}(M_{\odot})$ PBHs at high redshifts is likely to proceed via a thick inflated disk, best described within the ADAF framework — we will thus follow Ref. Poulin et al. (2017) in adopting this formalism. Typically, ADAF accreting models radiate extremely inefficiently (i.e., only a small fraction of the accreted energy is radiated and escapes into the IGM), while the remainder is dissipated locally in the form of heat through viscous interactions of the plasma. The relative efficiency with which X-rays are emitted in the IGM can be characterized by a parameter $\epsilon$ ($\epsilon\rightarrow 1$ being the limit that the accretion medium radiates perfectly, and $\epsilon\rightarrow 0$ being the limit in which most of the accreted energy goes into local heating via advection), which depends intimately on the geometry and strength of accretion. X-rays that escape into the medium have a large mean free path and would deposit their energy approximately uniformly on all scales, affecting both the temperature and ionization fraction of the IGM. On the other hand, the remainder of the energy is deposited locally (and is eventually swallowed by the PBH), modifying the temperature and ionization profile of the medium directly surrounding the PBH. Thus, accretion onto a single BH modifies the state of the IGM both globally (Section III.1) and locally (Section III.2), effects which for the sake of simplicity we treat separately in what follows.

In addition to the effects arising from accretion, the presence of a large population of PBHs can modify structure formation. The discrete nature of massive PBHs implies that shot noise fluctuations in the PBH number density are not necessarily negligible. These fluctuations manifest as isocurvature perturbations in the linear matter power spectrum, leading to an enhancement in structure at small scales Afshordi et al. (2003). For viable populations of PBHs (i.e., those evading current constraints), the effect is only expected to be pronounced for halos below the star-forming limit — these halos are typically referred to as minihalos. Within the context of $\Lambda$CDM, minihalos contribute minimally to the 21cm signal Iliev et al. (2002) (the signal is expected to be at the level of $\mathcal{O}(1)\,$mK in $\Lambda$CDM). However, it has been argued that the enhancement in the number density of these halos by a population of PBHs can make their contribution observable Gong and Kitajima (2017, 2018). We address this claim in Section III.3, and show that the contribution of minihalos is typically suppressed when the global effects of accretion are self-consistently included in the calculation.

III.1 Global heating and ionization

Here we address the global impact of accretion by PBHs on the IGM, i.e., the effect that arises from radiation that escapes the local vicinity of PBHs and deposits its energy approximately uniformly in the IGM. In what follows we consider the simplifying assumption of a monochromatic mass function for PBHs, broadly used in the literature. Note, however, that this is not necessarily a realistic assumption, since many mechanisms of PBH formation in the early Universe would generically lead to extended mass distributions, which could bias the obtained constraints Carr et al. (2016); Green (2016); Carr et al. (2017); Bellomo et al. (2018).

The rate of energy injected into the medium per unit volume is given by Ali-Haïmoud and Kamionkowski (2017); Poulin et al. (2017)

[TABLE]

where $L_{\rm acc}$ is the luminosity of radiation emitted from matter accreting onto a PBH of mass $M_{\rm PBH}$, $n_{\rm PBH}$ is the number density of PBHs of mass $M_{\rm PBH}$, $\rho_{\rm DM}$ is the density of dark matter, and $f_{\rm PBH}$ is the fraction of dark matter that is comprised of PBHs of mass $M_{\rm PBH}$. The accretion luminosity can be expressed in terms of the product of the accretion rate $\dot{M}_{\rm PBH}$ and the radiative efficiency $\epsilon$,

[TABLE]

In general, both the accretion rate and the radiative efficiency are complicated unknown functions of both the PBH properties and the evolution of the medium. However, various approximations have been developed that can be used to estimate each of these functions.

We begin by addressing the radiative efficiency factor $\epsilon$ in Eq. (7), which describes the fraction of energy radiated away into the IGM. In principle, it depends on the accretion rate, accretion geometry, and the properties of the medium. Here, we assume that the radiative cooling is inefficient and the dynamics of the accreting gas are controlled by advective currents (ADAF) Ichimaru (1977); Rees et al. (1982); Narayan and Yi (1994); Blandford and Begelman (1999); Narayan (2002); Pellegrini (2005); Xie and Yuan (2012); Yuan and Narayan (2014). It is worth noting that if the radiative cooling is efficient and forms a thin accretion disk (rather than the thick disk or tori that may appear in the ADAF scenario), the radiative efficiency could be orders of magnitude larger Shakura and Sunyaev (1973). For the efficiency function, we adopt the fitting formula obtained in Ref. Xie and Yuan (2012),

[TABLE]

where $\dot{M}_{\rm Edd}=10\,L_{\rm Edd}=10\times 4\pi\,GM_{\rm PBH}\,m_{p}/\sigma_{\rm T}\simeq 1.28\times 10^{39}(M_{\rm PBH}/M_{\odot})$ ergs/s is the Eddington rate (and $L_{\rm Edd}$ is the Eddington luminosity), and the functional form and coefficients are indicated in Tab. 1, assuming the viscous parameter to be $\alpha=0.1$ and the parameter which dictates the fraction of dissipation that heats the electrons directly to be $\delta=0.1$. Larger (smaller) values of $\delta$ are expected to increase (decrease) the radiative efficiency by a factor of a few in the most extreme scenarios.

We now shift our attention to the accretion rate. Accretion onto a point mass moving through some homogeneous medium was first studied by Hoyle and Lyttleton in the first half of the $20^{\rm th}$ century Hoyle and Lyttleton (1939, 1940a, 1940b, 1941). Nevertheless, this early treatment only considered the gravitational effects (neglecting, for example, the thermodynamic pressure), and was thus only useful as a first-step approximation within a particular limiting regime (far away from the point source and with a reasonably large relative velocity). Bondi and Hoyle were later able to correct for the previously neglected effect of pressure and calculated an exact solution to the problem of spherical accretion onto a stationary point mass in a homogenous medium Bondi and Hoyle (1944),

[TABLE]

where $\rho_{\infty}$ and $c_{s,\infty}$ are the density and speed of sound far away from the point source, $G$ is the gravitational constant, and $\lambda$ is an ${\cal O}(1)$ parameter quantifying non-gravitational effects Bondi and Hoyle (1944), which also allows to treat non-spherical accretion. Given that simulations suggest that radiative outflows reduce the amount of matter that ultimately accretes in the inner region, potentially reducing the observed luminosity by a factor of $\mathcal{O}(100)$ Xie and Yuan (2012), we account for this by taking $\lambda$ to be $0.01$.

In order to address the problem of accretion onto a BH moving with respect to the ambient medium with relative velocity $v_{\rm rel}$, Bondi suggested the substitution of $c_{s,\infty}$ by an effective velocity $\sqrt{c_{s,\infty}^{2}+v_{\rm rel}^{2}}$ Bondi (1952). The relative velocity between PBHs and the ambient medium can be described by the relative velocity between dark matter and baryons. Unfortunately, deriving the relative dark matter-baryon velocity at scales relevant for accretion is non-trivial Tseliakhovich and Hirata (2010); Dvorkin et al. (2014). On linear scales, the velocity field is Gaussian and thus, the modulus of the three-velocity field follows a Maxwellian distribution, with the square root of the variance approximated as Ali-Haïmoud and Kamionkowski (2017)

[TABLE]

The accretion luminosity is the result of averaging over the distribution of the relative velocities of PBHs with respect to the ambient gas. The accretion rate is inversely proportional to $\left(c_{s,\infty}^{2}+v_{\rm rel}^{2}\right)^{3/2}$ (after the substitution above) Ricotti et al. (2008); Ali-Haïmoud and Kamionkowski (2017) and, for the usual (simplified) case considered in the literature, the individual PBH accretion luminosity is proportional to the square of the accretion rate (i.e., $L_{\rm acc}\propto\dot{M}_{\rm PBH}^{2}$). It is important to emphasize that this applies only to scenarios in which $\epsilon\propto\dot{M}_{\rm PBH}$. Generalizing the accretion rate to scenarios in which the radiative efficiency has an arbitrary power law dependence (i.e., $\epsilon\propto\dot{M}_{\rm PBH}^{a}$) is obtained by expressing the effective velocity as

[TABLE]

with $\mathcal{M}\equiv\frac{\sqrt{\left<v_{L}^{2}\right>}}{c_{s,\infty}}$ and $U(x;y;z)$ the confluent hypergeometric function of second kind, or Tricomi’s function. Up to first order, the two limiting cases are:

[TABLE]

with $B(x,y)$ the beta function. Notice that this expression reproduces the limiting case in which $a=1$ found in Refs. Ricotti et al. (2008); Ali-Haïmoud and Kamionkowski (2017).

Therefore, using the effective velocity rather than the speed of sound, the accretion rate can be written as

[TABLE]

Even if more accurate expressions have been obtained El Mellah and Casse (2015), this modified accretion rate provides a reasonable description of the problem. It is worth noting that the Bondi accretion rate is strictly valid only for spherical accretion. It was recently argued, however, that accretion onto PBHs with $M_{\rm PBH}\gtrsim M_{\odot}$ likely results in the formation of an accretion disk Poulin et al. (2017). Note that this statement ran contrary to previous assumptions that PBH accretion was spherical Ricotti et al. (2008); Horowitz (2016); Ali-Haïmoud and Kamionkowski (2017); Aloni et al. (2017). Given that there exists a relatively large amount of uncertainty in the dynamics of disk accretion, some of which serves to increase and some of which serves to decrease the relative accretion rate Krumholz et al. (2005, 2006); El Mellah and Casse (2015), Eq. (13) can likely be thought of as a reasonable ball-park estimate. Therefore, despite all these caveats regarding the potential legitimacy of Eqs. (11) and (13), though, we will follow previous works Ali-Haïmoud and Kamionkowski (2017); Poulin et al. (2017) in adopting the expression of $v_{\rm eff}$ for $a=1$.

Finally, it is important to note that the quantity of interest is not the rate of energy injected into the medium, given in Eq. (6), but rather the rate of energy deposited into the medium. These two quantities are related by functions known as energy deposition factors,

[TABLE]

where $\epsilon_{c}^{\rm PBH}$ represents the energy deposition rate per baryon, and the subscript $c$ denotes the channel in which energy is deposited — i.e., ionization of hydrogen (HI) or helium (HeI), heating of the medium (heat), or excitations (${\rm Ly}\alpha$). These energy deposition factors are computed by using the following expression:

[TABLE]

where the $T_{c}(z,z^{\prime},\omega)$ are the transfer functions taken from Ref. Slatyer (2016). As is done in Ref. Poulin et al. (2017), the spectrum of the luminosity is taken from ADAF models of Ref. Yuan and Narayan (2014); namely, we adopt a simple parameterization given by

[TABLE]

where $\omega_{\rm min}\equiv(10\,M_{\odot}/M_{\rm PBH})^{1/2}$ eV, $\omega_{s}\sim\mathcal{O}(m_{e})$ (taken explicitly to be $\omega_{s}=200$ keV), and $\beta=-1$ (with reasonable values of $\beta\in[-1.3,-0.7]$). In principle, the exact form of the spectrum changes with the accretion rate. However, the results found here are largely independent of these details.

Eqs. (7) and (16) provide the ingredients necessary to compute the energy injected into the IGM from a population of PBHs as a function of redshift, and Eq. (15) describes how and where the energy is deposited. With this in hand, we now turn to incorporating this energy injection into the equations governing the evolution of the ionization fraction, the gas temperature, and the Lyman-$\alpha$ flux.

The quantity $x_{e}({\bf x},z)$, governing the evolution of the local ionized fraction of the neutral IGM, is given by

[TABLE]

where $n_{b}=\bar{n}_{b,0}(1+z)^{3}(1+\bar{\delta_{b}}({\bf x},z))$ is the baryon number density, $\Lambda_{\rm ion}$ is the ionization rate, $\alpha_{\rm A}$ is the case-A recombination coefficient, $\mathfrak{f}_{\rm H}=n_{\rm H}/n_{b}$ is the hydrogen number fraction and $C\equiv\langle n_{e}^{2}\rangle/\langle n_{e}\rangle^{2}$ is the clumping factor (set to one as default), with $n_{e}$ the electron number density. Ionization from the PBH accretion mechanism would lead to an additional contribution of the form

[TABLE]

where $\mathfrak{f}_{\rm H}=n_{\rm{H}}/n_{\rm{b}}$ and $\mathfrak{f}_{\rm He}=n_{\rm{He}}/n_{\rm{b}}$ are the hydrogen and helium number fractions, and $E_{\rm HI,HeI}$ are the ionization energies for hydrogen and helium. Here, we neglect the effect of secondary ionizations111Since the differential brightness temperature is proportional to the neutral hydrogen fraction, rather than the free electron fraction, the effect of ionizations from PBH accretion is expected to be significantly subdominant to that of heating. Thus, neglecting secondary ionizations is a safe assumption.. The evolution of the kinetic temperature of the gas $T_{k}$ is computed via

[TABLE]

where the last term accounts for the heating/cooling processes, with $\epsilon_{{\rm heat},j}$ the heating rate per baryon for the process $j$ (Compton cooling, X-ray heating and PBH heating).

Finally, there is also a contribution to the Lyman-$\alpha$ flux resulting from collisional excitations due to energy deposition by PBH accretion in the IGM, which can be written as

[TABLE]

where $\nu_{\alpha}$ is the emission frequency of a Lyman-$\alpha$ photon.

The aforementioned modifications to the free electron fraction, temperature of the gas, and Lyman-$\alpha$ flux are incorporated by modifying the publicly available codes cosmorec/Recfast++ Chluba and Thomas (2011); Ali-Haïmoud and Hirata (2010); Chluba et al. (2010), relevant at high redshifts when astrophysical mechanisms can be neglected, and 21cmFAST Mesinger et al. (2011), relevant at low redshifts when astrophysical mechansims, such as ionization and heating from stars and X-ray binaries, are relevant. Figures 1 and 2 show the redshift evolution of the free electron fraction and the kinetic gas temperature for a population of PBHs with mass $M_{\rm PBH}=100\,M_{\odot}$ and different relative abundances $f_{\rm PBH}$.222This range of values is simply intended to illustrate the dependence of these observables on $f_{\rm PBH}$. Note that the cases with the largest abundances are in tension with CMB constraints Poulin et al. (2017). Notice that the effect of PBHs accretion on the free electron fraction in Fig. 1 is clearly visible: the presence of the extra heating and ionization terms from PBHs accretion changes the redshift evolution of $x_{\rm e}$, increasing this quantity from the early recombination era, below $z\sim 1000$, until the late reionization era. The kinetic gas temperature would also be increased by the presence of the energy injection in the IGM (see Fig. 2). Similar to the case in which there is energy injection from dark matter annihilations Evoli et al. (2014); Lopez-Honorez et al. (2016), PBH accretion leads to an earlier and more uniform heating of the IGM, which is larger for an increasing fraction of dark matter in the form of PBHs, until stellar sources turn on and start to ionize the medium (around $z\sim 15$ in Fig. 1).333Note that although the spatial and redshift PBH distribution follows that of matter, it is different from the distribution of X-ray sources, i.e., star-forming halos beyond a threshold for atomic cooling. These results illustrate that even small abundances of PBHs could have dramatic effects on the properties of the IGM.

Before continuing, we would like to emphasize that the treatment of accretion adopted in this work is rather conservative. For the redshifts relevant for 21cm cosmology, the conditions necessary for the formation of accretion disks around PBHs seem likely. Within the context of disk accreting models, ADAF accretion is among the lowest in the radiative efficiency of X-rays. Adopting a larger radiative luminosity or accretion rate would correspondingly enhance the observable signatures associated with global heating and ionization of the IGM.

III.1.1 Global 21cm line signal

Since the optical depth at the frequencies of interest is small, one can expand the exponential term in Eq. (5) and express the global differential brightness temperature as

[TABLE]

where $\delta_{b}$ is the baryonic overdensity, $H$ is the Hubble parameter, $\partial v_{r}/\partial r$ is the peculiar velocity gradient along the line of sight, which introduces redshift space distortions, and $\Omega_{\rm m}$ and $\Omega_{\rm b}$ are the matter and baryon abundances of the Universe. A number of current and future experiments are devoted to detect the 21cm global signal averaged over all directions in the sky. Examples are EDGES and the future LEDA (Large Aperture Experiment to Detect the Dark Ages) Greenhill and Bernardi (2012) and DARE (Moon space observatory Dark Ages Radio Experiment) Burns et al. (2012). The observation of an absorption profile located at a redshift of $z\sim 17$ has been recently claimed by the EDGES experiment, with an amplitude which is twice the maximum predicted within the context of the $\Lambda$CDM model Bowman et al. (2018). This has motivated a larger number of studies in the literature Muñoz and Loeb (2018); McGaugh (2018); Barkana (2018); Barkana et al. (2018); Fraser et al. (2018); Kang (2018); Yang (2018); Pospelov et al. (2018); Costa et al. (2018); Slatyer and Wu (2018); Falkowski and Petraki (2018); Muñoz et al. (2018); Fialkov et al. (2018); Berlin et al. (2018); D’Amico et al. (2018); Safarzadeh et al. (2018); Hill and Baxter (2018); Clark et al. (2018); Cheung et al. (2019); Hektor et al. (2018b); Liu and Slatyer (2018); Hirano and Bromm (2018); Mitridate and Podo (2018); Mahdawi and Farrar (2018); Feng and Holder (2018); Ewall-Wice et al. (2018); Mirocha and Furlanetto (2019); Pospelov et al. (2018); Witte et al. (2018); Dowell and Taylor (2018); Lopez-Honorez et al. (2019); Jana et al. (2019); Fialkov and Barkana (2019); Ewall-Wice et al. (2019).

To compute the redshift evolution of the 21cm global signal we use of the publicly available code 21cmFAST Mesinger et al. (2011). Nevertheless, the code only evolves the IGM for redshifts $z\lesssim 35$. Thus, initial conditions for the mean ionization fraction $x_{e}$ and for the gas temperature $T_{k}$ are required inputs. We produce these inputs using a modified version of cosmorec/Recfast++ Chluba and Thomas (2011); Ali-Haïmoud and Hirata (2010); Chluba et al. (2010).444Notice that these initial values are different from the ones in the default tables, as the presence of PBHs would change significantly the recombination process and its outputs. The cumulative effect of a population of $100\,M_{\odot}$ PBHs on the global differential 21cm brightness temperature is shown in Fig. 3 for various PBH abundances. As discussed above, accretion onto PBHs significant increases both the free electron fraction and the temperature of the IGM (as shown in Figs. 1 and 2). Consequently, the effect of PBHs on the 21cm global differential brightness temperature can be summarized as a suppression of the absorption trough before X-ray heating dominates the signal. Notice the difference with the standard scenario, $f_{\rm PBH}=0$, depicted by the solid black lines in Fig. 3.

It is also possible to measure the power spectrum of the differential brightness temperature, i.e., to extract the fluctuations in the all-sky averaged differential brightness temperature. This method could in principle provide a cleaner foreground removal. Indeed, this is the major goal of experiments such as GMRT (Giant Metrewave Radio Telescope), LOFAR (LOw Frequency ARray), MWA (Murchison Widefield Array) and PAPER (Precision Array for Probing the Epoch of Reionization), which has improved the upper limits at $z=8.4$ Ali et al. (2015). As previously stated, future high-redshift 21cm experiments include SKA and HERA.555The 21cm line will also be observed at lower redshifts, predominately in the post-reionization era ($z\lesssim 3$), via the so-called intensity mapping technique, which will attempt to shed light on the spatial distribution emanating from dense clumps of neutral hydrogen Wyithe et al. (2008); Chang et al. (2008); Loeb and Wyithe (2008); Villaescusa-Navarro et al. (2014). This is the target of the GBT-HIM project Chang et al. (2016), CHIME (Canadian Hydrogen Intensity Mapping Experiment) Bandura et al. (2014), the Tianlai project Chen (2015) and SKA-mid frequency (see, e.g., Ref. Bull et al. (2015)). The dimensionless 21cm power spectrum is given by

[TABLE]

where $\delta^{D}$ is the Dirac delta function, the brackets indicate average quantities, and $\widetilde{\delta}_{21}(\mathbf{k},z)$ is the Fourier transform of ${\delta}_{21}(\mathbf{x},z)={\delta T}_{b}(\mathbf{x},z)/\overline{\delta T_{b}}(z)-1$. In what follows, we shall work with the 21cm differential brightness temperature power spectrum, $\Delta^{2}_{21}(k,z)=(k^{3}/2\pi^{2})P_{21}(k,z)$, and show results for $\overline{\delta T_{b}}^{2}(z)\Delta^{2}_{21}(k,z)$ in Section IV.1.

III.2 Local heating and ionization

As mentioned in the previous section, within the ADAF framework only a small fraction of the energy gained via accretion is radiated into the medium. In this section we attempt to model the contribution to the 21cm line from the local environment of the PBH, estimating the radial profiles.

III.2.1 Radial profiles

The difficulty in estimating the true contribution to the 21cm signal from the local medium surrounding PBHs arises from the fact that solving for the radial density, $\rho$, velocity, $v$, temperature, $T$, and ionization profiles, $x_{e}$, of the accreting medium is non-trivial, and can only be solved analytically in very specific scenarios and certain regions of the parameter space (for an analytic solution at small radii within the ADAF framework, see, e.g., Ref. Narayan and Yi (1995)). As will be shown below, the dominant contribution to the 21cm signal comes from radii near the Bondi radius (the Bondi radius, to be defined, approximately demarcates where the solutions asymptote to the IGM value). This motivates studying the simple solution to the problem of static spherically symmetric accretion around a point, in which we neglect heating from the radiated emission (likely a valid approximation for the ADAF framework since the efficiency is small) and viscosity (an assumption which certainly brakes down near the BH, but which plays a less significant role as one approaches the IGM).

In this static spherically symmetric case, the fluid equations for determining the radial profiles are given by

[TABLE]

where $\alpha_{\rm B}$ is the case-B recombination coefficient, $P=\rho\,(1+x_{e})T/m_{p}$ is the pressure, and $m_{p}$ is the mass of the proton (we neglect the helium contribution in this computation for the sake of simplicity). The first of these equations arises from integrating the continuity equation, the second one is the Euler equation (accounting for momentum conservation), and the final one determines the balance between ionization and recombination. As in Section III.1, we adopt the Bondi-Hoyle accretion rate given in Eq. (13). The relevant scales are the Bondi velocity, defined as in Eq. (11), and the Bondi radius

[TABLE]

Note that these definitions differ from those originally defined by Bondi, a consequence of the fact that we have modified the accretion rate via the substitution $c_{s}\rightarrow v_{\rm eff}$, as discussed above. For reference, a value of $v_{\rm eff}\sim c_{s}$ and $T_{k}\sim T_{\rm ad}\propto(1+z)^{2}$ lead to a Bondi radius of $r_{\rm B}\sim 3\times 10^{-5}\,(M_{\rm PBH}/M_{\odot})(30/(1+z))^{2}$ kpc. For distances greater than the Bondi radius, the effects of accretion are negligible, asymptoting to their background values.

As mentioned in Section III.1, the accretion of matter onto the central BH leads to a highly energetic radiation flux which is capable of ionizing the surrounding medium. In order to solve Eq. (25), it is necessary to identify the flux of ionizing photons as a function of the radial distance to the BH. Here, we assume that most of the radiation is produced at very small radii (where the surrounding medium can reach temperatures in excess of $10^{9}$ K), and thus, one can write the flux in terms of the comoving distance coordinate $r$ from the source by using a point source approximation:

[TABLE]

where $L_{\rm acc}(E,z)$ is the accretion luminosity given by combining Eqs. (7) and (16). The ionization optical depth $\tau_{\rm ion}$ is given by

[TABLE]

where the ionization cross section $\sigma(E)$ is taken to be the fitting function obtained in Ref. Zdziarski and Svensson (1989). For most of the radii of interest, the optical depth is found to be negligible. The local ionization rate is given by

[TABLE]

where $E_{\rm th}$ is the ground level energy of the hydrogen. We parameterize the energy deposition function $f_{\rm ion}$666Note that this is conceptually the same as parameter defined in Eq. (15). However, in the case of local ionization, it is appropriate to apply the on-the-spot approximation, and thus the functional forms and dependencies of these parameters are different. using different fitting formulas for $E\leq 0.5$ keV Dijkstra et al. (2004) and $E>0.5$ keV Shull and van Steenberg (1985).

Since, at the scales of interest, we can safely neglect heating/cooling effects, we get the usual adiabatic relation between the temperature and the density $T\propto\rho^{2/3}$. The set of Eqs. (23–25) can be easily solved numerically given the boundary conditions $\rho_{\infty}$ and $T_{\infty}$, which we take as the IGM values, that depend on both $M_{\rm PBH}$ and $f_{\rm PBH}$, and are calculated using the formalism presented in the previous section (Section III.1).

Is is important to bear in mind that after obtaining the temperature and density profiles, one must resolve the radiative transfer equation radially in order to obtain the differential brightness temperature. For these small enough length scales, the time evolution and expansion of the Universe can be neglected, leading to the following kinetic equation in terms of the brightness temperature at a frequency $\nu$, $T_{b,\nu}$:

[TABLE]

with $\kappa_{\nu}$ being the absorption coefficient

[TABLE]

where $B_{01},B_{10}$, and $A_{10}$ are the Einstein coefficients between the triplet and singlet states, and $\phi_{\nu}$ is the line profile, given by the Doppler broadening

[TABLE]

with $\Delta\nu_{\rm D}=\nu_{0}\sqrt{2\,T_{k}/m_{p}}$.

The brightness temperature along a line of sight at a radial distance $r$ from the PBH, Eq. (30), can be formally solved with the boundary condition at $r=0$, leading to

[TABLE]

with the optical depth $\tau_{\nu}(r^{\prime},r)$ defined as

[TABLE]

At distances near the PBH the medium would be fully ionized, leading to a vanishing $\kappa_{\nu}$. Therefore, it is appropriate to adopt the boundary condition $T_{b,\nu}(0)=T_{\rm CMB}$. Defining the differential brightness temperature as $\delta T_{b,\nu}(r)\equiv(T_{b,\nu}(r)-T_{\rm CMB})/(1+z)$, we get

[TABLE]

The spin temperature can be calculated using Eq. (2), with the collisional coupling coefficient, $y_{k}$, defined in Eq. (3) and the Lyman-$\alpha$ coefficient, $y_{\alpha}$, in Eq. (4). The Lyman-$\alpha$ flux in $y_{\alpha}$ contains two contributions, $J_{\alpha}(r)=J_{\alpha}^{\rm loc}(r)+J_{\alpha}^{\rm back}$. The first of these comes from the local emission in the local medium, given by

[TABLE]

We parameterize the energy deposition fraction $f_{\alpha}(x_{e})$ using the fitting formula obtained in Ref. Shull and van Steenberg (1985), given by $f_{\alpha}=0.48\,(1-x_{e}^{0.27})^{1.52}$. The second component to the Lyman-$\alpha$ flux comes from the background contribution, which is given the mean average over the IGM, produced by a combination of stellar sources and distant PBHs. This global flux $J_{\alpha}^{\rm back}$ for a given astrophysical model is computed by implementing the formalism of Section III.1 into 21cmFAST, as discussed in the previous section. Thus, the Lyman-$\alpha$ coupling coefficient can then be determined summing both components.

In Fig. 4, we show radial profiles of the density (top-left panel), gas temperature (purple line in the top-right panel), spin temperature (blue and red lines in the top-right panel), and the collisional ($y_{k}$, solid lines in the bottom panel) and Lyman-$\alpha$ ($y_{\alpha}$, dashed lines in the bottom panel) couplings at various redshifts for $M_{\rm PBH}=1\,M_{\odot}$. As designed, both the density profiles and the kinetic temperature profiles asymptote to the IGM values, which have been taken from the customized version of 21cmFAST (see Section IV.1 for more details). The subsequent 21cm temperature profile, computed using the values in Fig. 4, is depicted in Fig. 5. Note that the contribution must be truncated at large radii to avoid double counting the IGM contribution — we define this truncation to occur when the temperature profile is within $1\%$ of its asymptotic value, although the results are insensitive to a reasonable choice of this value. Note that the suppression of $y_{\alpha}$ (bottom panel in Fig. 4) between $z=20$ and $z=10$ can be explained by its dependence with the inverse of the temperature of the gas, $y_{\alpha}\propto T_{k}^{-1}$, which has been already heated significantly by this epoch.

III.2.2 Averaged 21cm line signal

The prescription provided above allows one to compute the 21cm profile for an individual PBH. Nevertheless, the observable quantity of interest is the contribution to the global 21cm line from all the PBHs in the Universe. Here, we re-derive the formula that allows one to compute the globally averaged brightness temperature $\overline{T_{b}}$ from the radial 21cm profile Iliev et al. (2002). The brightness temperature is related to the flux per unit frequency and solid angle $\frac{dF}{d\nu d\Omega}$ arising from a PBH of mass $M_{\rm PBH}$ by

[TABLE]

The averaged value of this quantity on a frequency interval $d\nu$ and a solid angle $d\Omega$ can be written as

[TABLE]

where $dV$ is the differential comoving volume, $n_{\rm PBH}$ is the number density of PBHs, and $\langle F\rangle$ is the area-averaged flux for a single PBH:

[TABLE]

with $A=\pi r_{\rm max}^{2}$, and where $r_{\rm max}$ is the maximum radius at which the PBH influences the IGM (defined as above — specifically, this is the radius at which $T_{k}$ is within $1\%$ of the asymptotic value). Using the fact that $|d\nu/dz|=\nu_{0}/(1+z)^{2}$, and writing the differential comoving volume as $dV=D^{2}dDd\Omega$, where $D$ is the comoving distance, we find

[TABLE]

where we have used $dD/dz=c/H(z)$. The integrated flux over frequency and solid angle is then given by

[TABLE]

where $\Delta\Omega=A/D_{A}^{2}$ is the solid angle subtended by the halo around the PBH, with $D_{A}=D/(1+z)$ the angular diameter distance, and $\Delta\nu_{\rm eff}=((1+z)\,\phi_{\nu}(\nu_{0}))^{-1}=\sqrt{\pi}\,\Delta\nu_{\rm D}/(1+z)$ Iliev et al. (2002). One can then revert to the averaged brightness temperature using

[TABLE]

From Eqs. (40) and (42), we find the globally averaged brightness temperature to be given by

[TABLE]

For the area-averaged differential brightness temperature $\delta\overline{T}_{b}$, Eq. (43) is applied but replacing $T_{b,\nu_{0}}$ by $\delta T_{b,\nu_{0}}$ computed from Eq. (35). We depict in Fig. 6 the area-averaged 21cm signal for one single PBH as a function of redshift (left panel), and the globally averaged signal (right panel), for various values of $M_{\rm PBH}$ and $f_{\rm PBH}$. The dashed and solid lines in the left panel indicate when the signal would be seen in absorption or emission. Interestingly, and perhaps counter-intuitively, the global contribution is often negative (i.e., it is seen in absorption). At first glance this appears strange, as the process of accretion heats and ionizes the medium — thus, one might naively expect to see this contribution in emission. However, the asymptotic signal is observed in absorption during this epoch and, while the process of accretion does heat the local medium, the density rises faster than the temperature (recall that $\rho(r)\propto T(r)^{3/2}$). Thus, there exists a range of radii near the Bondi radius for which the kinetic temperature has not yet been heated above the background temperature, but for which the density has noticeably increased — the net effect is an amplification of a slightly suppressed absorption dip.

The large drop in the signal for $z\lesssim 10$ in Fig. 6 appears because stellar sources have begun ionizing the medium. We emphasize that for the masses and dark matter fractions adopted here, the local contribution is never greater than $0.1$ mK. Given that the usual IGM signals for these epochs range between $1$ and $100$ mK (see section IV.1), the local contribution from PBHs is negligible for the models of interest.

III.2.3 Comment on the local contribution

Before continuing, a discussion on the treatment and impact of the local contribution is in order. We begin by addressing a number of concerns regarding the formalism adopted in the section above.

It should be emphasized that the spherically symmetric treatment of accretion adopted here neglects a number of potentially important features. First, our formalism only accounts for the adiabatic evolution of the medium, leading to the relation $T\propto\rho^{2/3}$. Thus, we have neglected the local heating of the gas by the highly energetic radiation generated during accretion. This heating term is proportional to the luminosity, and therefore, to the radiative efficiency and the PBH mass. However, the effect of such a heating would be the enhancement of the local temperature, which would suppress the absorption signal that arises near the Bondi radius, subsequently reducing the global contribution to the 21cm signal. Thus, from this point of view, neglecting the heating term is a conservative assumption. We have also chosen to neglect viscosity in our treatment, despite the fact that viscosity is known to play an important role in ADAF accretion Xie and Yuan (2012); Narayan and Yi (1994). We believe that this choice is justified as the derived signal seems to be dominated by the behavior at larger radius, where the effects of viscosity are reduced. Nevertheless, a more rigorous numerical treatment would be required to verify this statement.

Another concern about the adopted formalism is the assumption of the static solution of an isolated PBH in a uniform background. In reality, this is almost certainly not the case. The matter accreted onto the PBH must come from the IGM. We have implicitly assumed, however, that the asymptotic condition of the IGM remains unaltered by the presence of the PBH. Conceptually, one would expect that an underdensity at $r\gtrsim r_{\rm B}$ is necessary in order to counteract the overdensity observed for $r\lesssim r_{\rm B}$. Neglecting momentarily the effects of heating and ionization, the contribution of the underdensity should identically cancel that of the overdensity. The fact that the signal is predominately seen in absorption indicates that it is dominated by the state of the IGM near the Bondi radius, where not properly accounting for these changes in density could be important. Thus, it is entirely possible that the observed absorption feature is actually artificial, and arises from the fact that we have not correctly accounted for the influence of the PBH on the IGM. Nevertheless, this implies that our formalism is likely to significantly overestimate the true contribution of local accretion, again suggesting the estimate is conservative.

Before continuing, it is important to address why the results found here differ by orders of magnitude from previous PBH calculations Tashiro and Sugiyama (2013); Bernal et al. (2018). We think there are multiple reasons for this discrepancy,777In addition to the differences in formalism, discussed below, the de-excitation rates (taken from Ref. Kuhlen et al. (2006)) appear to be incorrectly transcribed in Ref. Bernal et al. (2018). Most notably, the logarithms in the e-H de-excitation rate should be base 10, rather than natural logs. We have verified, using the formalism of Ref. Bernal et al. (2018), that this mistake does indeed enter the calculations. Implementing the correct de-excitation rate into the same formalism completely removes the absorption dip responsible for the majority of the signal (see Fig. 6 of Ref. Bernal et al. (2018)), and thus significantly alters the contribution of minihalos to the global 21cm signal. The check we have performed agrees with the result of Ref. Tashiro and Sugiyama (2013), which appears to have a formalism nearly identical to that of Ref. Bernal et al. (2018), but without the artificial absorption dip. which we outline below.

First, note that the assumed PBH luminosity in Ref. Bernal et al. (2018) is near the Eddington limit (specifically they adopt a luminosity of $L=0.1\,L_{\rm Edd}$). Such a large luminosity is not consistent with ADAF accretion, but rather solutions in which one forms a thin disk that radiates extremely efficiently ($\epsilon\sim 1$). In the event that only the radiative efficiency is modified, local heating would be reduced, and global heating enhanced. However, there also exists a significant difference in the normalization of the luminosity, as the dimensionless accretion rate adopted in our framework is roughly $L/L_{\rm Edd}\sim 10^{-10}$. The difference in adopted luminosity changes the relative importance of heating, implying that within our formalism one cannot neglect the derivative terms in the accretion equations, whereas for large luminosities these terms may be subdominant to local heating. Next, it is slightly concerning that Ref. Bernal et al. (2018) seeks a static solution to the local profile on scales comparable to and higher than the Hubble length (note that $r_{\rm H}=c/H(z)\simeq 4.6\times 10^{4}\left[30/(1+z)\right]^{3/2}$ kpc). Hence, the expansion of the Universe can no longer be neglected, and thus the static framework is applied out of the range of validity. This does not appear to be a problem with the ADAF framework, as the smaller luminosities imply smaller scales of interest.

Finally, it is important to note that one should not directly apply Eq. (21) to obtain the brightness temperature profile around the PBH. The Eq. (21) is only valid when considering large scales; namely it is valid when one can neglect the radial dependence of both the spin temperature and the optical depth, and when the optical depth is exclusively given by the optical depth of the IGM. Here, neither of these statements are true, and thus one must solve the radiative transfer equation as a function of the radial coordinate using Eq. (35) (as was done in previous works, see, e.g., Refs. Iliev et al. (2002); Gong and Kitajima (2017)). This profoundly impacts the strength of the signal. In Fig. 7 we illustrate, using the formalism of the previous section, that incorrectly applying Eq. (21) leads to an enhancement of the signal by roughly $\sim 2$ orders of magnitude.

Therefore, we conclude that it is unlikely that the local contribution from PBHs could significantly add to the 21cm signal. Nevertheless, it is likely the case that complex high-resolution simulations will need to be performed for this contribution to be truly understood.

III.3 Shot noise

Since PBHs are point sources, they contribute to the power spectrum as Poisson white noise (i.e., they contribute in a scale-independent manner) Afshordi et al. (2003). The associated perturbations are isocurvature modes and thus, their contribution only affects scales smaller than those corresponding to the matter-radiation equality era, $k>k_{\rm{\rm eq}}\equiv a_{\rm eq}H_{\rm eq}$, being negligible at larger scales. Approximating the isocurvature transfer function as $T_{\rm iso}=3/2\,(1+z_{\rm eq})$ for $k>k_{\rm{\rm eq}}$ and [math] otherwise, where $z_{\rm{eq}}$ is the redshift of matter-radiation equality Peacock (1999), this term can be written as Afshordi et al. (2003)

[TABLE]

where $\Omega_{\rm DM}$ is the fraction of dark matter energy density over the critical density $\rho_{c}$ and $D(z)$ is the linear growth factor of density perturbations. The inclusion of this contribution modifies the halo mass function on small scales. For the masses studied here (namely $M_{\odot}\lesssim M_{\rm PBH}\lesssim 10^{3}M_{\odot}$), however, the effects are only important for halos which are not large enough to cool and collapse to form stars (such halos are known as ‘minihalos’). Since star-forming halos are expected to produce the dominant 21cm signal in the epochs of interest, one is naively driven to conclude that the Poisson noise in the power spectrum is not relevant and can be neglected (even though, we have included it in the calculation of the global signal). It was recently argued Gong and Kitajima (2017, 2018), however, that the enhanced number density of minihalos arising from the shot noise induced by $\sim\mathcal{O}(M_{\odot})$ PBHs could sufficiently elevate the contribution of minihalos to a discernible level. We briefly outline the formalism for estimating the minihalo contribution below, and show that the computation of Ref. Gong and Kitajima (2017) failed to self-consistently account for the heating of the IGM induced by PBH accretion. Correctly accounting for this effect tends to suppress the minihalo contribution to the point where it can be neglected.

III.3.1 Brightness temperature inside minihalos

A dark matter halo at redshift $z$ can be characterized by one parameter, the overdensity $\Delta_{c}(z)$ relative to the critical density $\rho_{c}(z)$, so that its mass $M_{\rm h}$ and virial radius $R_{\rm vir}$ are related by

[TABLE]

Here, we use the virial overdensity based on spherical collapse Bryan and Norman (1998)

[TABLE]

where

[TABLE]

In the calculation below, we assume the neutral hydrogen fraction, $x_{\rm H}$, to be given by the IGM value, and the kinetic temperature in the halo by $T_{k}=T_{\rm vir}$, where the virial temperature of a truncated isothermal halo is Shapiro and Iliev (1999)

[TABLE]

being $\mu=1.22$ the mean molecular weight of neutral IGM. The average number density of neutral hydrogen of the halo is

[TABLE]

where $\rho_{b}$ the baryon energy density. For the sake of consistency with the constant temperature assumption, we consider that neutral hydrogen in halos follows a truncated (singular) isothermal distribution, normalized to provide the mean value given in Eq. (49),

[TABLE]

We have also verified that reasonable modifications to this assumption (e.g., taking an NFW profile Navarro et al. (1997) or a truncated non-singular isothermal profile Shapiro and Iliev (1999)) give rise to only small changes on the cumulative minihalo signal.

Since the flux of Lyman-$\alpha$ photons predominantly arises from star formation, the spin temperature inside minihalos, whose virial temperatures are below threshold for atomic hydrogen cooling, is dictated exclusively by the collisional coupling, $y_{k}$, given by Eq. (3). The brightness temperature for photons with frequency $\nu$ through a minihalo of given mass and redshift is obtained by solving the radiative transfer equation. Here, however, one has to be careful so as to not double count factors contained in the IGM contribution. To clarify, consider decomposing the IGM and minihalo contribution to the spin temperature such that

[TABLE]

where $\tau_{\rm IGM}$ and $\tau_{\rm mh}$ are the optical depth of the IGM and the minihalo, respectively, and we define

[TABLE]

where $r$ is the radial distance from the center of a minihalo with maximum radius $R_{\rm vir}$. The general solution to the radiative transfer equation (neglecting scattering) is given by

[TABLE]

Using Eqs. (51–54), one finds

[TABLE]

where $s$ is the line-of-sight. In order to avoid double counting, one must subtract off the contribution of the IGM, given by

[TABLE]

implying the solution for the minihalo brightness temperature is given by

[TABLE]

where we have added some explicit dependences. The optical depth of the minihalo, $\tau_{\rm mh}(\nu,s)$, is given by

[TABLE]

and $s$ and $\ell$ are radial comoving distances, $s$ being the component along the line-of-sight and $\ell$ connecting $s$ to the center of the minihalo, $\ell^{2}=s^{2}+b^{2}R_{\rm vir}^{2}$ Chongchitnan and Silk (2012), with $b$ the impact parameter. One might be inclined to obtain the differential brightness temperature by subtracting $T_{\rm CMB}$ from Eq. (58), however note that it is given by

[TABLE]

and thus the proper contribution to the total differential brightness temperature is simply given by Eq. (58).

Obtaining the global contribution from all minihalos can be done in a similar manner to what was shown in Section III.2. Specifically, one can use Eq. (43), but with the generalization that $n_{\rm PBH}\rightarrow\int\,dM_{\rm h}\,dn/dM_{\rm h}$, i.e.,

[TABLE]

where, in analogy to the area-averaged signal in Section III.2, the brackets represent the average over the minihalo cross section (c.f., Eq. (39)), $A=\pi R_{\rm vir}^{2}$. The integration of Eq. (61) is bounded from above by the mass threshold for star formation, which we take as the one corresponding to a virial temperature, Eq. (48), of $T_{\rm min}=10^{4}$ K. The lower limit of the integral is chosen as the maximum between the PBH mass and the Jeans mass $M_{\rm J}$ (although for all cases of interest here, $M_{\rm PBH}<M_{\rm J}$), which is given by the mass enclosed in a region with diameter of the Jeans length, $\lambda_{\rm J}=\sqrt{(5\pi\,T_{k})/(3\,G\,\bar{\rho}\,\mu\,m_{p})}$,

[TABLE]

where $\bar{\rho}$ is the mean matter density prior to collapse. Here, we calculate $T_{k}$ using the formalism presented in Section III.1.

The number density of halos per unit mass is given by the halo mass function,

[TABLE]

where $f(\sigma)$ represents the fraction of mass that has collapsed to form halos per unit interval in $\ln\sigma^{-1}$, with $\sigma$ the root-mean-square of the matter density fluctuation smoothed with a real-space top hat filter over the virial radius,

[TABLE]

where $D(z)$ is the growth factor of linear perturbations, $P(k)$ is the linear matter power spectrum (here, modified by including the shot noise term from the PBH population, Eq. (44)), and $W(k;M_{\rm h})$ is the Fourier transform of the real-space top hat filter. Here we use the Sheth-Tormen halo mass function Sheth and Tormen (1999); Sheth et al. (2001),

[TABLE]

with fit coefficients $A=0.3222$, $a=0.707$, and $p=0.3$, and the critical collapse density $\delta_{c}=1.686$. The halo mass function at $z=20$ for several PBH cases is depicted in Fig. 8. Note that, since the additional shot noise contribution from PBH to the matter power spectrum depends on the combination $f_{\rm PBH}M_{\rm PBH}$ (see Eq. (44)), so does the modification of the halo mass function with respect to the standard $\Lambda$CDM case.

Using Eqs. (58), (61), and (63), one can calculate the globally averaged 21cm signal from the cumulative population of minihalos. In Fig. 9, we depict the contribution of minihalos assuming a monochromatic population of PBHs with mass of either $10~{}M_{\odot}$ or $10^{3}~{}M_{\odot}$ and a dark matter fraction of $f_{\rm PBH}=10^{-2}$ or $f_{\rm PBH}=10^{-4}$. Notice that the redshift dependence of the minihalo contribution can be rather non-trivial, as the $T_{k}$ dependence of the Jeans mass is modified by the accretion model, and the mass derivative of the halo mass function can be either enhanced or suppressed relative to $\Lambda$CDM at this mass scale. We can see that the enhancement on the formation of structures due to the Poisson noise would produce a larger signal in the case $M_{\rm PBH}=10~{}M_{\odot}$ and $f_{\rm PBH}=10^{-2}$. However, contrary to previous claims Gong and Kitajima (2017, 2018), the other cases shown in Fig. 9 are not enhanced but suppressed. For small PBH masses and low values of $f_{\rm PBH}$, the minihalo contribution saturates to the standard $\Lambda$CDM prediction; that is to say, the number density of minihalos is negligibly modified, and the result is equivalent to that of $\Lambda$CDM. Meanwhile, the signal is suppressed for models with larger values of $M_{\rm PBH}\times f_{\rm PBH}$, where the enhancement of the halo mass function is significant. The reason for the discrepancy between the result presented here and that found in Refs. Gong and Kitajima (2017, 2018) lies on the treatment of the Jeans mass. In those previous works, only adiabatic cooling was considered for the evolution of the IGM temperature, while here the heating by standard X-ray sources and by accretion onto the population of PBHs (see Section III.1) has also been included. This implies a larger Jeans mass and thus, a reduction of the relevant mass range of integration in Eq. (61). For instance, for the models depicted in Fig. 2, X-ray radiation would heat up the temperature around an order of magnitude with respect to the adiabatic cooling case at $z\sim 10$, increasing therefore the Jeans mass by a factor of $10^{3/2}\sim 30$.

IV Methodology

IV.1 Numerical simulations

As already mentioned above, in order to compute the redshift evolution of both the 21cm global signal and power spectrum we make use of the publicly available package 21cmFASTv1.2 Mesinger et al. (2011). This code produces realizations of the evolved density, ionization, peculiar velocity and spin temperature fields from semi-analytic calculations. It depends on a number of parameters describing the different processes taking place at the reionization period. We use here a minimal set of four parameters: the ultraviolet (UV) ionization efficiency, $\xi_{\rm UV}$, the number of X-ray photons per solar mass, $\xi_{\rm X}$, the minimum virial temperature of halos hosting galaxies, $T_{\rm min}$, and the number of photons per stellar baryon between Lyman-$\alpha$ and the Lyman limit, $N_{\alpha}$. 888The default value of $N_{\alpha}$ in 21cmFAST is obtained by assuming Pop II stars Barkana and Loeb (2005) and normalizing their emissivity to $\sim 4400$ ionizing photons per stellar baryon. The phenomenological parameter $\xi_{\rm UV}$ is assumed to be proportional to: (i) the fraction of ionizing photons escaping their host galaxies, (ii) the number of ionizing photons per stellar baryons inside stars, and (iii) the fraction of baryons that form stars. The efficiency for ionization, heating, and Lyman-$\alpha$ production by X-ray sources depends on the total X-ray emission rate, which is proportional to the star formation rate and to the number of X-ray photons per solar mass in stars, $\xi_{\rm X}$.999A value of $\xi_{\rm X}=10^{56}\,M_{\odot}^{-1}$ implies $N_{\rm X}\simeq 0.1$ X-ray photons per stellar baryon. The temperature $T_{\rm min}$ is the minimum of the virial temperature of a halo, below which gas cannot cool efficiently, and thus star formation is suppressed.101010Note that the definition of the virial temperature here is taken to be as defined in Ref. Barkana and Loeb (2001), which differs from Eq. (48) by a factor of $2.44$. This difference arises from the fact that in Eq. (48) we consider $\mu=1.22$ and that the virial temperature depends upon the adopted halo profile, although these changes do not significantly impact any of the results of this section. The default value in the 21cmFAST code is $T_{\rm min}=10^{4}$ K, which has been identified in the literature with the atomic cooling threshold Evrard (1990); Blanchard et al. (1992); Tegmark et al. (1997); Ciardi et al. (2000); Haiman et al. (2000), and using Ref. Barkana and Loeb (2001), it corresponds to a minimum halo mass of $M_{\rm min}\simeq 8\times 10^{7}\,M_{\odot}$ at a redshift $z=10$. In the following, we present various examples of the 21cm signal using a fiducial $\Lambda$CDM model defined by $(\xi_{\rm UV},\,\xi_{\rm X},\,T_{\rm min},\,N_{\alpha})=(50,\,2\times 10^{56}\,M_{\odot},\,5\times 10^{4}\,K,\,4000)$.

In Fig. 3 we showed the global differential 21cm brightness temperature for a population of $100\,M_{\odot}$ PBHs, for a range of abundances, as a function of redshift. In Fig. 10 we depict the 21cm power spectrum as a function of redshift at scales $k=0.15$ Mpc*-1* (left panel) and $k=0.4$ Mpc*-1* (right panel), which are expected to be reasonably free from foregrounds Pober et al. (2013a). Notice that the 21cm power spectrum in the standard scenario (solid black curves) exhibits three characteristic peaks, associated to the the epochs of reionization, heating from X-ray sources and Lyman-$\alpha$ pumping Pritchard and Furlanetto (2007); Mesinger et al. (2011); Baek et al. (2010), from lower to higher redshifts. The presence of PBHs is translated into a suppression of the X-ray heating peak in the power spectrum, which is obviously more pronounced as one increases $f_{\rm PBH}$. Additionally, the Lyman-$\alpha$ pumping peak could even disappear for $f_{\rm PBH}\gtrsim 10^{-2}$. In Fig. 11 we show the 21cm power spectrum as a function of the scale for two fixed redshifts: $z=12.2$ (left panel) and $z=17.5$ (right panel). The red areas indicate the scales where the signal is expected to be contaminated by foregrounds. The effect of PBHs arises primarily from an increase in the flux of X-rays, which can induce both scale-dependent and redshift-dependent features in the power spectrum.

In Figs. 10 and 11 we also include the forecasted errors associated to future HERA and SKA measurements of the 21cm power spectrum. These errors have been estimated by means of the publicly available code 21cmSense111111https://github.com/jpober/21cmSense Pober et al. (2014, 2013b) (see also Ref. Parsons et al. (2012)). The total noise is given by

[TABLE]

with two contributions, one from thermal noise ($N$) plus a second one, a sample variance error ($S$), $\overline{\Delta}^{2}_{21}\equiv\overline{\delta T_{b}}^{2}\Delta^{2}_{21}$. The thermal noise depends on the solid angle, the integrated observation time and the temperature of the system. We focus here on the future HERA Beardsley et al. (2015) and SKA-low frequency Mellema et al. (2013) experiments. We consider both the intermediate and final HERA configurations, with 127 and 331 antennas, which we refer to henceforth as Hera 127 and Hera 331. For SKA-low frequency, we follow the design presented in the SKA System Baseline Design Document Mellema et al. (2013). Finally, we assume an exposure of 1080 hours and a bandwidth of 8 MHz, as these are the default values for these parameters in 21cmSense.

One of the fundamental difficulties associated with full statistical 21cm analyses arises from the fact that consistent and accurate theoretical predictions for the 21cm signal require time consuming semi-analytical (or even fully hydrodynamical) simulations. Restricting our attention momentarily to a single PBH mass, obtaining a coherent interpolation over a five dimensional parameter space (defined by $f_{\rm PBH},\,\zeta_{\rm UV},\,\zeta_{\rm X},\,T_{\rm min}$, and $N_{\alpha}$) requires enormous amounts of computing time. Recently, various groups have developed techniques which rely on simplified calculations (at the expense of numerical accuracy), principle component decompositions, and machine learning algorithms, to reduce the computational expense involved in obtaining a comprehensive 21cm parameter scan (see, e.g., Refs. Kern et al. (2017); Shimabukuro and Semelin (2017); Gillet et al. (2019); La Plante and Ntampaka (2018); Jennings et al. (2019); Doussot et al. (2019)).

In this work, we circumvent this numerical issue by exploiting a class of feed-forward neural networks known as multilayer perceptrons (MLPs). Specifically, for a fixed PBH mass, we compute the full 21cm power spectrum for $\sim 600$ choices of astrophysical parameters (constrained to the range defined in Tab. 2), and construct MLPs with two hidden layers, each containing $\sim 50$ hidden nodes, to emulate the output of 21cmFAST for arbitrary choices of parameters. The MLP is trained on $\sim 70\%$ of the computed power spectra while the remaining $30\%$ are simultaneously used to ensure the neural network is not over-learning. We find that this procedure reproduces the computed datasets (both the trained dataset and test dataset), as well as various power spectra computed with randomly generated points in parameter space after the training of the neural network. We choose to assess the relative accuracy of the neural net by defining an accuracy statistic, which later serves as the quantity directly entering our likelihood, as

[TABLE]

where $\bar{\Delta}^{2}_{21}(k,z)$ is the power spectrum computed by 21cmFAST, and ${\Delta}^{2}_{21,NN}(k,z)$ is the power spectrum as computed by the neural net. Following Ref. Greig and Mesinger (2015), we have chosen to add in quadrature an additional modeling error (here we conservatively adopt a $40\%$ error, although we note that this factor is somewhat adhoc) that is intended to capture the approximate level of disagreement between 21cmFAST calculations and full hydrodynamic simulations. Using Eq. (67), one finds that the mean error of the neural net is $\lesssim 1\%$.

IV.2 Statistics

We adopt a multivariate Gaussian likelihood with a diagonal covariance matrix, the elements of which are set using Eq. (67), and further assume that measurements are obtained in six log-spaced bins in $k-$space from $k=0.15~{}{\rm Mpc}^{-1}$ to $k=1~{}{\rm Mpc}^{-1}$ and 9 log-spaced measurements in redshift from $z\sim 8.3$ to $z\sim 19.5$. While it is likely that experiments such as SKA will be capable of achieving far better resolution in both $k-$ and $z-$space, adopting a finer grid without properly accounting for correlated errors risks the possibility of significantly overestimating the sensitivity of these experiments.

We adopt the four-parameter astrophysical model discussed in Section IV.1 (implying a total of 5 model parameters for each fixed value of $M_{\rm PBH}$), where the priors on each parameter are as given in Tab. 2. For each fixed value of $M_{\rm PBH}$, we perform a Markov chain Monte Carlo (MCMC) assuming the true measurement is given by a fiducial astrophysical model with $f_{\rm PBH}=0$.

V Results and conclusions

The $2\sigma$ upper limits as derived from the MCMC analysis described in Section IV.2 are shown in Fig. 12 for HERA (red line), assuming 331 antennas, and SKA (green line). Also shown in Fig. 12 are current constraints from microlensing surveys Allsman et al. (2001); Tisserand et al. (2007); Wyrzykowski et al. (2011) (blue line), and the CMB Ali-Haïmoud and Kamionkowski (2017); Horowitz (2016); Poulin et al. (2017) (purple and black lines), where the purple contour denotes the bound derived assuming an identical accretion model to the one adopted here. As discussed in the introduction, there exist additional constraints in this region, however these are not shown for clarity. The MCMC analysis was only performed for PBH masses in the range $(1-10^{3})~{}M_{\odot}$, however this is a consequence of the computational difficulty associated with 21cm analyses. In principle, these experiments will be sensitive to both smaller and larger masses, and the relative importance of the shot noise, minihalos, and the modification to the IGM, in general, varies with mass.

The results shown in Fig. 12 suggest that near-future 21cm experiments will be able to increase the sensitivity to $\mathcal{O}(M_{\odot})$ PBHs relative to that of the CMB by more than one order of magnitude. Moreover, 21cm surveys provide a highly complementary probe that is sensitive to the dynamics of accretion at much lower redshifts. Thus, should a positive detection be made with both probes, one may be able to shed light on the dynamics and evolution of accretion in these systems. For completeness, we also show in Fig. 13 the result of the MCMC performed for $M_{\rm PBH}=10^{3}\,M_{\odot}$ using the SKA telescope array. Fig. 13 illustrates that many of the astrophysical parameters adopted in this analysis could be well constrained by these experiments.

In this work we have focused on understanding the sensitivity that near-future radio interferometers may have to a population of $\gtrsim\mathcal{O}(M_{\odot})$ PBHs using the 21cm power spectrum. In particular, we have jointly analyzed three effects: namely, heating and ionization arising from accretion on both local and global scales, and the signal arising from minihalos (whose number density is enhanced with respect to that of $\Lambda$CDM, a consequence of having non-negligible Poissonian noise in the distribution of PBHs). Contrary to previous claims in the literature Tashiro and Sugiyama (2013); Bernal et al. (2018), we find that the local effect of accretion is negligible for all redshifts, masses, and PBH fractions discussed in this work. Moreover, for the PBH masses studied here, the addition of shot noise to the power spectrum modifies the halo mass function below the mass threshold for efficient star formation. Consequently, this term would only be relevant if minihalos gave sizable contributions, or if one considers extremely heavy PBHs (however it is worth noting that there unavoidably exists uncertainty in the high redshift halo mass function within $\Lambda$CDM which could produce degenerate effects to minihalos). While previous work has claimed that the enhanced number density of minihalos, produced by viable populations of PBHs, could give rise to an observable signal Gong and Kitajima (2017, 2018), we show here that these calculations had not self-consistently accounted for the global heating of the IGM, which significantly raises the Jeans mass and tends to suppress the minihalo contribution. Thus, the dominant effect on the 21cm power spectrum arises exclusively from the modifications to the heating and ionization of the IGM. We find experiments like HERA and SKA could significantly improve upon existing limits from the CMB, should observations be inconsistent with a $f_{\rm PBH}=0$ cosmology (i.e., no PBHs).

Finally, we emphasize that 21cm observations will have access to much more information than just the two-point correlation function. Should foregrounds be removed to a sufficiently high degree, these experiments may be able to exploit higher dimensional n-point correlation functions, or perhaps one could exploit the power of convolutional neural nets to shed light of the highly non-Gaussian behavior of neutral hydrogen during these epochs.

Acknowledgments

OM, PVD and SW are supported by the Spanish grant FPA2017-85985-P of the MINECO. SPR is supported by a Ramón y Cajal contract, by the Spanish MINECO under grant FPA2017-84543-P, and partially, by the Portuguese FCT through the CFTP-FCT Unit 777 (UID/FIS/00777/2019). All the authors also acknowledge support from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreements No. 690575 and 674896. This work was also supported by the Spanish MINECO grant SEV-2014-0398.