<figure><img src="image_1.png" /></figure>

HAL
open science

# Cosmic microwave background bounds on primordial black holes including dark matter halo accretion

Pasquale D. Serpico, Vivian Poulin, Derek Inman, Kazunori Kohri

► To cite this version:

Pasquale D. Serpico, Vivian Poulin, Derek Inman, Kazunori Kohri. Cosmic microwave background bounds on primordial black holes including dark matter halo accretion. Phys.Rev.Res., 2020, 2 (2), pp.023204. 10.1103/PhysRevResearch.2.023204. hal-02518035

HAL Id: hal-02518035

https://hal.science/hal-02518035

Submitted on 13 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Cosmic microwave background bounds on primordial black holes
including dark matter halo accretion

Pasquale D. Serpico●,¹,* Vivian Poulin,² Derek Inman●,³ and Kazunori Kohri⁴,⁵,⁶,⁷

¹Univ. Grenoble Alpes, USMB, CNRS, LAPTh, F-74940 Annecy, France

²LUPM, CNRS & Université de Montpellier, F-34095 Montpellier, France

³Center for Cosmology and Particle Physics, Department of Physics, New York University, 726 Broadway, New York, New York 10003, USA

⁴Theory Center, IPNS, KEK, Tsukuba 305-0801, Ibaraki, Japan

⁵The Graduate University of Advanced Studies (Sokendai), Tsukuba 305-0801, Ibaraki, Japan

⁶Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Department of Physics, Oxford, OX1 3PU, United Kingdom

⁷Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

<figure><img src="image_3.png" /></figure>

(Received 2 March 2020; accepted 29 April 2020; published 26 May 2020)

Even if massive ($10 M_{\odot} \lesssim M \lesssim 10^4 M_{\odot}$) primordial black holes (PBHs) can only account for a small fraction of the dark matter (DM) in the universe, they may still be responsible for a sizable fraction of the coalescence events measured by LIGO/Virgo, and/or act as progenitors of the supermassive black holes (SMBHs) observed already at high redshift ($z \gtrsim 6$). In the presence of a dominant, non-PBH DM component, the bounds set by CMB via an altered ionization history are modified. We revisit the cosmological accretion of a DM halo around PBHs via toy models and dedicated numerical simulations, deriving updated CMB bounds which also take into account the last Planck data release. We prove that these constraints dominate over other constraints available in the literature at masses $M \gtrsim 20–50 M_{\odot}$ (depending on uncertainty in accretion physics), reaching the level $f_{\text{PBH}} < 3 \times 10^{-9}$ around $M \sim 10^4 M_{\odot}$. These tight bounds are nonetheless consistent with the hypothesis of a primordial origin of the SMBH massive seeds.

DOI: 10.1103/PhysRevResearch.2.023204

# I. INTRODUCTION

Despite several decades of direct, indirect, or collider searches, we are still ignorant of the nature of mark matter (DM) of the universe. It is even unclear if it is made of a single species or if, just like its baryonic matter counterpart, it is constituted by different components. Based on current constraints, primordial black holes (PBH) formed in the early universe remain a viable DM candidate only in the mass window $10^{-16} M_{\odot} \lesssim M \lesssim 10^{-11} M_{\odot}$ (see, e.g., Refs. [1–3]). However, PBH of much higher masses, even if not constituting the bulk of DM, can still have other interesting cosmological and astrophysical consequences. For instance, PBH with masses of tens of $M_{\odot}$ could be responsible for some if not most of the heavy BH mergers discovered by LIGO/Virgo, even if they only contribute a fraction $f_{\text{PBH}} \sim \text{few} \times 10^{-3}$ to the DM [4–6].

One wonderful probe for such putative pristine objects is the cosmic microwave background (CMB). PBH would start accreting gas soon after matter-radiation equality ($z \simeq 3500$); the gas heats up to the point that ionizing radiation is

emitted, which alters the opacity of the gas in the long period between recombination and reionization. CMB temperature and polarization fluctuations are extremely sensitive to this phenomenon and can be used to constrain it. This argument has been used by several authors in the past, from the pioneering study [7] to the amended and more recent calculations in Refs. [8–12].

In this paper, we extend our previous calculation [11] to account for the accretion of the dominant, non-PBH DM particles onto PBH, which enhances the baryonic accretion and eventually the PBH luminosity constrained by the CMB, simply as a result of the increased gravitational potential felt by the baryons. Besides providing a more realistic assessment of the CMB bounds on stellar-mass PBH, we are also motivated by the possibility that intermediate mass PBH ($10^2 M_{\odot} \lesssim M \lesssim 10^4 M_{\odot}$) may provide the seeds leading to the supermassive black holes (SMBHs) observed at redshift $z \gtrsim 6$ (with the current record-holder of $10^{8.9} M_{\odot}$ at $z = 7.54$ [13]) whose origin constitutes a long-standing mystery, see, e.g., Ref. [14]. To the best of our knowledge, nobody has assessed the viability of this hypothesis in the light of current CMB anisotropy bounds.

This paper is structured as follows. In Sec. II, we start by reviewing the formalism and hypotheses entering the bounds on PBHs set by the CMB. Since the main systematic uncertainty in these bounds consists of the treatment of accretion physics, we derive the bounds in a couple of physically motivated benchmarks, which should bracket the uncertainties. In Sec. III, we discuss the role played by DM halos

<footnote>* serpico@lapth.cnrs.fr</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
accreting around PBH, treating them both via semianalytical toy models and dedicated numerical simulations. Our CMB bounds are reported in Sec. IV, providing an upgrade to the bounds presented in Ref. [11] in the light of the last Planck data release and the updated treatment of the energy release via the new EXOCLASS package [15]. Also, we present an extension of the bounds up to masses of $M \simeq \text{few} \times 10^4 M_\odot$, beyond which some working hypotheses break down. Section V briefly reviews the puzzle concerning the origin of SMBH and discusses the implications of our CMB limits. We can anticipate that, even under Eddington accretion conditions for the PBH surrounded by DM halos, CMB angular power spectra are not capable of testing the hypothesis that the SMBH detected already at high redshift are (seeded by) PBH. Hence this remains a viable conjecture with interesting consequences for the cosmology of the dark ages on which we briefly comment upon in Sec. VI, where we also report our conclusions.

## II. THE LUMINOSITY OF ACCRETING PBH

In the following, we assume that the ionization of gas in the dark ages due to accretion onto PBHs and probed by the CMB can be considered homogeneous. A justification is provided in Appendix A. The key input needed to compute the CMB bound is the total energy injection rate per unit volume:

$$ \frac{d^2E}{dVdt} = L_{\text{acc}}n_{\text{PBH}} = L_{\text{acc}}f_{\text{PBH}}\frac{\rho_{\text{DM}}}{M}, \quad (1) $$

where $L_{\text{acc}}(M, t)$ is the bolometric (in general, time dependent) accretion luminosity onto a PBH of mass $M$ and is the main unknown. A milder uncertainty (within a factor $\lesssim 2$) comes from the spectral distribution of the emitted radiation, which in turns determines the energy deposited in the medium—what actually matters—and for which we make use of the transfer functions<fnref n="1" /> from Refs. [17,18] as implemented in the EXOCLASS package [15]. As in Ref. [11], we assume that the energy-differential spectrum of $L_{\text{acc}}$, $L_\omega$, is dominated by Bremsstrahlung emission (see Refs. [8,19,20]), with a mildly decreasing frequency dependence over several decades, and a cutoff given by the temperature of the medium near the Schwarzschild radius $T_s$. Namely, we adopt

$$ L_{\omega} \propto \omega^{-a} \exp(-\omega/T_s), \quad (2) $$

where $T_s \sim \mathcal{O}(m_e)$ (we use 200 keV in the following, like in Ref. [11]) and $|a| \lesssim 0.5$, again like in Ref. [11] ($a=0$ was used in Ref. [8]). Before addressing the question of the effective PBH mass evolution, let us review typical benchmarks for $L_{\text{acc}}$. Also, note that if the time dependence in $M(t)$ is not negligible, then $f_{\text{PBH}}$ may become time dependent as well. As we will argue below, this is not the case for the redshift range of interest. If comparing the CMB bounds with other, low-$z$

bounds, one should take into account that our $f_{\text{PBH}}$ denotes the initial DM mass fraction in the form of PBH.

### A. The Eddington limit

The Eddington luminosity $L_E$ is the luminosity at which accretion is balanced by radiation pressure in a spherical system and can be simply computed as

$$ L_E \equiv \frac{4\pi \mu G M m_p c}{\sigma_T} \approx 1.26 \times 10^{38} \frac{M}{M_\odot} \text{erg/s}, \quad (3) $$

with $\sigma_T$ the Thomson cross section, $m_p$ the proton mass, and $\mu$ a chemical composition dependent parameter, equal to unity for pure hydrogen. Under the rather extreme hypothesis that $L_{\text{acc}} = L_E$, we see from Eq. (1) that the energy injected per unit volume per unit time scales like a constant times the matter dilution factor of the universe, so that the bounds should be independent of $M$. Note that $L_E$ is a quantity linear in the accreting object mass and allows one to introduce a fundamental accretion timescale (independent of the mass of the object) known as Salpeter time,

$$ \tau_E \equiv \frac{Mc^2}{L_E} = \frac{\sigma_T c}{4\pi \mu G m_p} \approx 0.4 \text{ Gyr.} \quad (4) $$

As we will illustrate more quantitatively in Sec. VA, under standard assumptions the PBH mass accretion timescale is about one tenth of Eq. (4) for a BH shining at the Eddington luminosity. On the other hand, the age of the universe at the most relevant redshifts for CMB bounds, $300 \lesssim z \lesssim 600$, is 1–3 Myr, more than a factor 100 smaller than Eq. (4). For less extreme accretion rates, the hierarchy between the accretion timescale and the age of the universe is a fortiori even bigger. As a result, the PBH mass can be considered constant and equal to its initial value at the epochs relevant for deriving CMB bounds.

### B. More realistic accretion scenarios

The function $L_{\text{acc}}$ in Eq. (1) is usually parameterized in terms of the two following quantities:

(i) $\dot{M}$, the matter accreted per unit time onto the PBH.

(ii) $\epsilon$, the overall efficiency of conversion of accreted matter into radiation, in terms of which one writes:

$$ L_{\text{acc}} \equiv \epsilon \dot{M} c^2. \quad (5) $$

Concerning $\dot{M}$, an analytical theory exists in two limiting cases, both applying to a homogeneous gas of mass density $\rho_\infty$: The hypothesis of a stationary, spherical symmetric accretion of a body at rest [21] (Bondi), and the purely ballistic limit (i.e., accounting only for gravity, no hydrodynamical nor thermodynamical effects included) of a point mass moving at a constant speed $v_{\text{rel}}$ in the gas [22–24] (Hoyle-Lyttleton). These limiting cases justify the following “Bondi-Hoyle-Lyttleton” (BHL) parametrization

$$ \dot{M}_{\text{BHL}} = 4\pi\lambda\rho_{\infty} \frac{(GM)^2}{v_{\text{eff}}^3}, \quad (6) $$

where

$$ v_{\text{eff}}^2 \equiv v_{\text{rel}}^2 + c_s^2, \quad (7) $$

<footnote n="1">¹More recent tools have been developed for the computation of the energy deposited in the cosmological plasma during the cosmic dark ages [16], but it has been shown that the effect of energy injection onto the CMB bounds is accurately described by the tabulated transfer functions.</footnote>
$c_s$ is the speed of sound in the homogeneous matter of density $\rho_\infty$, and $\lambda$ is a dimensionless coefficient dependent upon environmental parameters, of $\mathcal{O}(1)$ for $v_{\text{rel}} \gg c_s$ (Hoyle-Lyttleton), and a calculable function assuming values of $\mathcal{O}(0.1-1)$ (see, e.g., Ref. [8]) in the limit $c_s \gg v_{\text{rel}}$ (Bondi). For more realistic situations, Eq. (6) is often used, but with $\lambda$ now intended as an adjustable parameter or function fitted, e.g., to simulation results. The cross section for spherical accretion onto a pointlike particle is usually described in terms of the Bondi radius, which is the distance from the center at which the escape velocity equals the sound speed. Hence, it makes sense to define a “generalized” Bondi radius,

$$r_B \equiv \frac{GM}{v_{\text{eff}}^2}, \qquad (8)$$

in terms of which Eq. (6) writes

$$\dot{M}_{\text{BHL}} = 4\pi\lambda\rho_{\infty} v_{\text{eff}} r_{\text{B}}^2, \qquad (9)$$

and which, just like Eq. (8), reduces to Bondi’s results for $v_{\text{eff}} \to c_s$, while smoothly interpolating to the Hoyle-Lyttleton regime for larger velocities.

Concerning $\epsilon$, it can be computed in spherical symmetry under some assumptions for the radiative processes, the most up-to-date treatment being provided in Ref. [8]. In that case, it assumes rather small values, of the order of $10^{-5}$. In the case of disk accretion, a typical benchmark value considered in the literature is 0.1. In the following, we adopt for $\lambda$ and $\epsilon$ the same prescriptions used in Ref. [11] and already implemented in the EXOCLASS package.

One of the main unknowns in the cosmological problem at hand is the actual relative velocity of PBH and baryons. The most obvious velocity scale in the problem is the sound speed in the baryon fluid, $c_s$. At large spatial scales, a larger velocity $v_L$ is predicted in linear (but nonperturbative) theory [25], but it is questionable if that applies down to the small scales relevant for accretion, where the PBH potential dominates and the DM fluid approximation breaks down [11]. Despite some first studies in that respect suggesting negligible effects in the epoch of interest [26], gas simulations are definitely required to account for the dissipative nature of the baryonic gas accreting on the proto-halos. We remain agnostic on the question and consider two cases:

(i) If $v_L$ is a reasonable proxy, the PBH-baryon motion at the relevant redshifts is supersonic ($M \sim 2-5$) and, as argued for instance in Refs. [27,28], the accretion should then be disklike.<fnref n="2" /> In this case, our fiducial parametrization closely follows [11], with

$$v_{\text{eff}} \approx \sqrt{c_s \sqrt{\langle v_L^2 \rangle}}, \qquad (10)$$

where, at $z \lesssim 1000$, we adopt [8,25]

$$c_s \approx 6 \frac{\text{km}}{\text{s}} \sqrt{\frac{1+z}{1000}}, \quad \sqrt{\langle v_L^2 \rangle} \approx 30 \frac{\text{km}}{\text{s}} \left( \frac{1+z}{1000} \right). \qquad (11)$$

(ii) If, on the other hand, the relevant small-scale motions are subsonic (i.e., $M \lesssim 1$), for the low values of $f_{\text{PBH}}$ considered in the following the accretion should be better approximated by a spherical one. In this case, a conservative, spherical accretion scenario is adopted, with Bondi accretion and $v_{\text{eff}} \simeq c_s$, with $c_s(z)$ obeying Eq. (11). We also assume the limiting case of purely collisional ionization, known to yield the most conservative bounds [8].

### C. Regimes of validity

The treatment just described has two limitations:

(i) The hypothesis of stationarity, i.e., the system settles down in the Bondi steady state fast compared to the cosmological expansion:

$$\frac{r_B}{v_{\text{eff}}} H(z) < 1. \qquad (12)$$

As already argued in Ref. [29], this leads to the requirement $M \lesssim \text{few} \times 10^4 M_\odot$, which is the upper limit for which we present our results.

(ii) How to deal with super-Eddington accretion, i.e., $L_{\text{acc}} > L_E$, which in our formalism can be attained for sufficiently heavy PBH. It is still debated what happens under these accretion conditions, which are sensitive to multidimensional effects and realistic disk radiation spectra (for a recent study, see Ref. [30]). It is clear however that several phenomena come into play: For instance, radiative and kinetic feedback can break stationary conditions, with outflows and episodic periods of very high luminosity alternating with a long period of low accretion and luminosity. Or quasi-steady-state, super-Eddington mass accretion can take place, with a corresponding drop in efficiency in order to satisfy $L_{\text{acc}} \lesssim L_E$. While we will comment again on this regime in Sec. V, we address the reader to reviews such as Ref. [31] for details and a more complete picture. In the following, we will adopt the prescription to cap luminosity at $L_E$ whenever the formalism yields nominally $L_{\text{acc}} > L_E$.

Interestingly enough, super-Eddington accretion is attained in our formalism for $M \gtrsim 10^4 M_\odot$, so that both conditions above yield similar limitations, albeit by coincidence. The homogeneous approximation discussed in Appendix A is also valid in the same range of interest. It is worth clarifying that CMB anisotropy bounds are expected to exist also at higher masses, but they become rather uncertain and definitely the formalism above is insufficient to tackle them. Fortunately, for such high masses, other bounds become relevant, as discussed in Sec. V B.

In Appendix B, we also check that the dynamical friction that a rather massive PBH experiences moving supersonically in the cosmological baryonic gas is negligible for the masses and redshifts of interest for this work.

### III. INCLUDING COSMOLOGICAL DM HALOS

It has been argued in the past [32–34] that, due to the PBH gravity, a DM halo would form around massive PBH, boosting their accretion. Note that the Eddington luminosity benchmark only applies to baryons, subject to radiation pressure. As far as baryonic accretion is the only one considered, the Salpeter timescale suggests that the PBH mass remains

<footnote n="2">²Note that this argument is independent of the angular momentum due to binaries or small-scale motion discussed in Ref. [11], hence applies also to small $f_{\text{PBH}}$.</footnote>
essentially constant down to the redshifts of interest for CMB bounds. Hence, we can safely consider $f_{\text{PBH}}$ constant, while DM halos affect the phenomenology via the altered accretion rate.

Although the original Bondi problem was considering accretion onto a point particle, a natural generalization of the notion of Bondi radius for an extended distribution of mass can be written as [35]:

$$ \frac{G_N M_{\text{PBH}}}{r_{\text{B,eff}}} - \Phi_h(M_{\text{PBH}}, r_{\text{B,eff}}, t) = v_{\text{eff}}^2(t), \quad (13) $$

where $r_{\text{B,eff}}$, the effective Bondi radius, is the unknown, $M_{\text{PBH}}$ is the initial PBH mass and $\Phi_h$ the (time-dependent) gravitational potential associated to the DM halo. Our treatment of the problem consists of adopting Eq. (9) but with $r_B$ replaced by $r_{\text{B,eff}}$, the solution of Eq. (13) with the gravitational potential of the halo estimated analytically or numerically. In this section, we revisit our estimates accounting for the DM capture phenomenon. First, we consider a toy model, which we believe determines an upper limit to the effect. Then, we improve over this estimate with the help of numerical simulations.

### A. Toy model

The most optimistic scenario for PBH growth is that DM is exactly cold and with no dispersion, and the PBH is the only center of attraction in the whole universe. This is a spherically symmetric problem. In order to calculate the time evolution of a radius $r$ of a mass shell around a PBH which encloses different species, we solve the following differential equation,

$$ \frac{d^2 r}{dt^2} = -\frac{4G_N}{\pi} 3r \left[ \rho_{\text{PBH}} + \sum_i (\rho_i + 3p_i) \right], \quad (14) $$

where $\rho_i$ and $p_i$ are the energy density and pressure of a component “i,” respectively, and we defined the energy density of the PBH as $\rho_{\text{PBH}} = 3M_{\text{PBH}}/(4\pi r^3)$. Here $i$ runs on all the components of radiation and matter (and dark energy if it were effective). The physical radius $r$ is represented by $r = a(t)x$ where $a(t)$ is the scale factor normalized to be $a(t_0) = 1$ at the present Universe ($t = t_0$), and $x$ is the co-moving coordinate. At each time $t$, the bound (or halo) mass is equivalent to the DM density up to the radius $r_s$ defined by

$$ \frac{dr_s}{dt}(t) = 0, \quad (15) $$

although a similar value would be found if one derives $r_s$ from the condition that the density at $r_s$ is twice the cosmological background one, as in Ref. [33]. Under the above-mentioned approximations, we find good agreement with the results reported in Ref. [33], namely:

(i) A time evolution given by

$$ M_{\text{halo}} \approx \left( \frac{3000}{1+z} \right) M_{\text{PBH}}. \quad (16) $$

(ii) A density profile proportional to $\propto r^{-3}$, as illustrated in Fig. 1, down to the distances (not resolved in Fig. 1) where a free-fall profile $r^{-3/2}$ takes over.

Equation (16) should be understood as an upper limit to the mass growth of the PBH via DM accretion proceeding self-similarly once a DM halo of mass larger than the PBH is

<figure><img src="image_4.png" /><figcaption>FIG. 1. Density profile of halo around a PBH and integrated halo mass as a function of comoving radius. We plot the cases at $z = 30$ for $M_{\text{PBH}} = 1, 10, 10^2, 10^3, 10^4, 10^5, 10^6 M_{\odot}$.</figcaption></figure>

accumulated. Its breakdown is only expected at very late times (e.g., when dark energy kicks in) or when the hypothesis of isolated PBH breaks down (which further requires high $f_{\text{PBH}}$). On the other hand, the radial profile crucially depends on the free-fall boundary condition at the center. Not accounting for the DM angular momentum is however a very crude approximation. It was speculated in Ref. [29], Sec. 4, that other scaling solutions like the ones described in the seminal paper by Bertschinger [32] may provide a better description of the results. To verify this conjecture, we turn to the results of N-body simulations.

### B. Numerical simulations

We perform cosmological N-body simulations using a version of the CUBEP$^3$M code [36] modified to include PBHs as a separate particle species co-evolved with a generic collisionless DM candidate [37]. We opt for homogeneous DM initial conditions at $a = 10^{-6}$ and select cosmological parameters consistent with Planck: $\Omega_c = 0.26$, $\Omega_b = 0.05$, and $z_{eq} = 3374$, with $\Omega_{\text{PBH}} = f_{\text{PBH}}\Omega_c$. Baryons are not evolved and instead are assumed to be homogeneous inducing errors of order $\Omega_b/\Omega_m \simeq 15\%$. With this setup, simulations are invariant to the numerical masses of the particles, instead being sensitive to their ratio:

$$ \frac{M_{\text{PBH}}}{M_{\text{DM}}} = \frac{f_{\text{PBH}}}{1 - f_{\text{PBH}}} \frac{N_{\text{DM}}}{N_{\text{PBH}}}. \quad (17) $$

After running a simulation, we can set $M_{\text{PBH}}$ to a physical value which then fixes the volume, $L^3$, of the simulation:

$$ M_{\text{PBH}} = \bar{\rho}_{\text{cr}} \Omega_c f_{\text{PBH}} \frac{L^3}{N_{\text{PBH}}}, \quad (18) $$

where $\bar{\rho}_{\text{cr}}$ is the comoving critical density.
<figure><img src="image_5.png" /><figcaption>FIG. 2. DM halo profiles (with mass scaling in x axis) around a PBH at $z = 99$ and $z = 999$, for different values of the PBH abundance. A power-law profile $r^{-2.25}$ is also shown for comparison.</figcaption></figure>

We would like our simulations to have the best possible length resolution (i.e. the smallest box size $L$). From Eq. (18) we see that, at fixed $M_{\text{PBH}}$ and cosmological parameters, this is achieved by minimizing the number of PBHs in the simulation, i.e. $N_{\text{PBH}} = 1$. By doing this we no longer accurately follow Poisson fluctuations in the PBH density field; however, from the Epstein mass function [38,39] describing the Poisson distribution we can deduce that PBHs rarely interact when $f_{\text{PBH}} \ll (1+z) \times 10^{-4}$. We performed an explicit test of this by running a simulation with $f_{\text{PBH}} = 10^{-5}$ and $N_{\text{PBH}} = 100$ and found the resulting profiles comparable, but noisier, to the single PBH case. To accurately model the isolated halo growth we require that the DM halo be composed of many DM particles when it becomes comparable to $M_{\text{PBH}}$:

$$f_{\text{PBH}} \gg \frac{1}{1 + N_{\text{DM}}} \approx N_{\text{DM}}^{-1} \approx 4 \times 10^{-9}, \quad (19)$$

where the last equality utilises our maximum number of particles: $2 \times 512^3$. We therefore consider $f_{\text{PBH}} = 10^{-5}, 10^{-6}$, and $10^{-7}$ with $N_{\text{DM}} = 2 \times 384^3, 2 \times 512^3$, and $2 \times 512^3$. The simulations are run from $a = 10^{-6}$ to $a = 10^{-2}$. We have tested how much halting the energy injection at this redshift affects our constraints and found it to be around a percent.

Given our assumption that the PBHs are isolated, we expect the DM halo to be independent of $f_{\text{PBH}}$ apart for numerical resolution effects. We show the density profile as a function of $f_{\text{PBH}}$ at $a = 10^{-3}$ and $10^{-2}$ in Fig. 2. As expected, we find that the profiles are independent of $f_{\text{PBH}}$. We also find that the best resolution is obtained for the $f_{\text{PBH}} = 10^{-5}$ simulation, despite the fact that it has slightly fewer particles. We therefore use the results from this simulation. We also see the effects of redshift on the profile. Over almost two decades in radius, the profile at early times matches the $r^{-2.25}$ power law predicted by Ref. [32]; at late times, it is only slightly steeper, moving closer to $r^{-2.5}$. This power law profile is consistent with those found in (independent) numerical

<figure><img src="image_6.png" /><figcaption>FIG. 3. Gravitational potentials of the PBH (dotted inclined lines, with the straight line representing the analytical result), the halo (curved dashed line) and the sum (solid line) compared to the typical velocity scales of the problem (horizontal lines) at $z = 99$.</figcaption></figure>

simulations by Adamek et al. [40], who also find a smooth transition to standard NFW-like profile at large radii.

Our chief goal is determining $r_{\text{B,eff}}$. To do this, we first interpolate the particles to a grid using the Cloud-in-Cell method and then solve Poisson’s equation for the gravitational potential:

$$\phi_i = (4\pi G)^{-1} \nabla^{-2} \rho_i, \quad (20)$$

where $i$ can indicate PBH and CDM separately and $\nabla^{-2}\rho$ is evaluated in Fourier space. We then find $r_{\text{B,eff}}$ by plugging the potential above in Eq. (13). We show the obtained potentials at $z = 99$ in Fig. 3, with $r_{\text{B,eff}}$ defined as the intersection of the total potential (solid black) with the horizontal gray line. We remark some numerical artifacts: For the PBH potential we know the exact result, $-\phi_{\text{PBH}}(r)/(GM_{\text{PBH}}) = 1/r$. However, we see that the result differs from this analytical result on both small and large scales: On small scales this is due to the interpolation error, whereas on large scales it is due to periodic boundary conditions. Typically, the PBH contribution to the halo is negligible whenever the DM distribution is important. The periodicity artifact leads to a slight underestimates of $r_{\text{B,eff}}$. At high redshifts the PBH is much more relevant, $r_{\text{B,eff}}$ is smaller and the numerical solution may lead to a slight overestimate of the solution. It is interesting to note that the numerical results lead to an estimated halo mass which about 60% of the simple result of Eq. (16), with a similar scaling with redshift, although with a different mass profile. These results motivate the semianalytical model described in the following section, which we later use to cover PBH masses whose Bondi radii are not resolved by our simulations.

Note that it is safe to neglect the “ordinary” DM halos feedback onto the halos growing around PBH, since the former ones only grow at much later times (typically $z \lesssim 30$ in a $\Lambda$CDM cosmology) than those of concern for us. A fortiori, the feedback of the baryons can also be neglected. The bulk of baryons has significant velocity at the epochs of interest, and they are still kinetically coupled to the CMB. Most of them are essentially unbound to halos, and their ratio to the DM in the
growing halos around PBH is much smaller than the baryon to dark matter cosmological density ratio of ~15%. Hence, objections on the realism of power law DM density profiles around BH surviving in the current universe [41] do not apply to the pristine configurations considered here.

### C. Semianalytical model

In the specific case of a pointlike potential due to the PBH plus the power-law matter distribution around it, with density $\rho(r) \propto r^{-\alpha}$ up to a distance $r_h$ and total mass $M_h$, Eq. (13) rewrites

$$v_{\text{eff}}^2(z) = \frac{G_N M_{\text{PBH}}}{r_{\text{B,eff}}} + \frac{G_N M_h}{r_{\text{B,eff}}} \left\{ \Theta(r_{\text{B,eff}} - r_h) + \frac{\Theta(r_h - r_{\text{B,eff}})}{1-p} \left[ \left( \frac{r_b}{r_h} \right)^p - p \left( \frac{r_{\text{B,eff}}}{r_h} \right) \right] \right\}, \quad (21)$$

where $p = 3 - \alpha$, and $M_h$ and $r_h$ depend from $\{M_{\text{PBH}}, z\}$.

We adopt Eq. (16) for the halo mass within the turnaround radius, where the turnaround radius is (see, e.g., Sec. 4 in Ref. [29])

$$r_{\text{t.a.}} \approx 58 \, \text{pc} (1+z)^{-1} \left( \frac{M_h (M_{\text{PBH}}, z)}{M_{\odot}} \right)^{1/3}. \quad (22)$$

We identify $r_h = r_{\text{t.a.}}$, in order to have a self-consistent normalization of the mass.

Equation (21) admits either the solution

$$r_{\text{B,eff}} = \frac{G_N (M_{\text{PBH}} + M_h)}{v_{\text{eff}}^2} \approx \frac{G_N M_h}{v_{\text{eff}}^2} \equiv r_{\text{B,h}}, \quad (23)$$

which holds if $r_h < r_{\text{B,h}}$; otherwise, if $r_h > r_{\text{B,h}}$, neglecting the PBH mass one has

$$r_{\text{B,eff}} \approx r_h \left[ (1-p) \frac{r_h}{r_{\text{B,h}}} + p \right]^{\frac{1}{p-1}} \le r_h. \quad (24)$$

Note that Eq. (24) tends to $r_{\text{B,h}}$ when $p \to 0$, as expected: When the DM halo profile is very steep and/or the halo is very compact, as far as accreting baryons are concerned they simply see a BH whose effective mass is the sum of the PBH and the DM halo mass. If the halo is fluffy or large, only a fraction of the mass of the halo contributes to the accretion. In any case, the condition $r_{\text{B,eff}} \ge r_{\text{B,PBH}}$ must hold. This constraint must be verified and eventually imposed by hand as a lower limit if using the approximated Eq. (24) or the RHS of Eq. (23). We have found that the CMB constraints obtained using this model are in agreement within 50% with the ones obtained from results of the numerical simulations in the mass range covered by the simulations.<fnref n="3" /> We thus use this model with $p = 0.75$ to compute the impact of PBH accretion onto the CMB.

<figure><img src="image_7.png" /><figcaption>FIG. 4. Impact of accreting PBH with and without DM halos on the lensed CMB TT [panel (a)], TE [panel (b)], and EE [panel (c)] power spectra. The reference is $\Lambda$CDM with $\{\omega_b = 0.02218, \omega_{\text{cdm}} = 0.1205, 100 * \theta_s = 1.04069, \tau_{\text{reio}} = 0.055, \ln(10^{10}A_s) = 3.056, n_s = 0.9619\}$. We consider a disk accretion scenario and set the PBH fraction to the constraints at 95% C.L. derived in this work.</figcaption></figure>

## IV. CMB CONSTRAINTS

### A. Impact of accretion onto DM halos

We modify the branch EXOCLASS [15] of the public code CLASS to include the effect of DM halo. We compare the effect of accretion with and without halos. In practice, our simulations only have the necessary resolution to solve Eq. (13) for redshifts $100 < 1+z < 1000$, which encompasses the redshift range at which e.m. energy deposition has the biggest impact on the CMB, $300 < 1+z < 600$ [17]. We extrapolate with constant values of the lower (higher) boundary at lower (higher) redshifts. We checked that this has subpercent impact by turning off injection at $z < 100$, while the effect of energy ejection is naturally turned off at higher $z$ since the plasma is still mostly ionized. We show the effect of the e.m. energy injection from accretion of matter around PBH on the CMB power spectra in Figs. 4 and 5 for $M_{\text{PBH}}/M_{\odot} = 100, 1000$. Since we find that including the halos increases the impact of the PBH on the power spectra by up to ~2 orders of magnitude (at fixed fraction), we actually compare two values of $f_{\text{PBH}}$ such that each case shows the 95% C.L. exclusion.

<figure><img src="image_8.png" /><figcaption>FIG. 5. Same as Fig. 4, in the spherical accretion case.</figcaption></figure>

<footnote n="3">³We checked that varying $p \in [0.50, 0.75]$ affects our results to below 10% level.</footnote>
Interestingly, because of different time dependence of the energy injection, the shape of the EE and TE CMB power spectra residuals is significantly different (see, e.g., Ref. [11] for a review of the effect of e.m. energy injection). This shows that the effect of DM halos could potentially be distinguished (and thus should be taken into account) if a signal were detected. It might even tell us something about the nature of DM as halos may not form depending on DM properties (e.g., if DM is warm or fuzzy). Further work will be required to accurately characterize the signal from PBH accretion given the sensitivity of future CMB experiments to polarization anisotropy.

## B. Analysis

We run a Markov-chain Monte Carlo (MCMC) using the public code MONTEPYTHON-V3.4 [42,43], interfaced with our modified version of CLASS. We perform the analysis with a Metropolis-Hasting algorithm, assuming flat priors on {$\omega_b, \omega_{cdm}, \theta_s, A_s, n_s, \tau_{reio}, f_{PBH}$} at fixed PBH mass $M/M_\odot = [10, 20, 30, 40, 50, 10^2, 10^3, 10^4, 10^5]$. We additionally perform a MCMC run at fixed $f_{PBH} = 1$, with $M/M_\odot$ free to vary to determine the minimal PBH mass probed by cosmological data through accretion. In practice we find that $M_{min}^{95\%} = 15M_\odot > 10M_\odot$ in the spherical case, we therefore did not run with $M = 10M_\odot$ in that accretion scenario. We adopt the Planck collaboration convention and model free-streaming neutrinos as two massless species and one massive with $M_\nu = 0.06$ eV. Our data set includes Planck 2018 high-$\ell$ and low-$\ell$ TT, EE and lensing likelihood [44,45]; the isotropic BAO measurements from 6dFGS at $z = 0.106$ [46] and from the MGS galaxy sample of SDSS at $z = 0.15$ [47]; the anisotropic BAO and the growth function $f_{\sigma_8}(z)$ measurements from the CMASS and LOWZ galaxy samples of BOSS DR12 at $z = 0.38, 0.51$, and 0.61 [48]. Additionally, we use the Pantheon<fnref n="5" /> supernovae dataset [49], which includes measurements of the luminosity distances of 1048 SNe Ia in the redshift range $0.01 < z < 2.3$. As usual, we use a Choleski decomposition [50] to deal with the numerous nuisance parameters associated with the likelihoods (not recalled here for brevity). We consider chains to be converged using the Gelman-Rubin [51] criterion $R - 1 < 0.05$. We perform four sets of runs, assuming either spherical or disk accretion, and absence or presence of a DM halo around the PBHs. For simplicity, we adopt a monochromatic PBH mass function, keeping in mind that for extended mass functions (which are to be generically expected from single-field inflationary models [52]) bounds typically tighten [53,54], as we explicitly checked for the CMB ones in our previous article [11].

As a warm-up, we derive the limit in case PBH are accreting at Eddington luminosity. As argued, in this case one obtains a mass-independent bound, which reads

$$f_{\text{PBH}} < 2.9 \times 10^{-9} \quad (L_{\text{acc}} = L_{\text{E}}). \qquad (25)$$

This is an optimistic benchmark for what is presumably the best limit that CMB can yield to. Our more realistic con-

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 6. Bounds on the abundance of PBH assuming disk accretion [panel (a)] or spherical accretion [panel (b)]. We show the results with (light-shaded) and without (dark-shaded) the formation of a DM halo. The horizontal line shows the limiting bound of Eq. (25). See text for details.</figcaption></figure>

straints at 95% C.L. are shown in Fig. 6. Note that the bounds on PBH in the absence of DM halos (dark shaded regions) are themselves stronger than bounds previously derived in [11] by a factor ~4. This improvement is due roughly equally to the new Planck 2018 low multipole polarization data and the additional use of BAO and Pantheon data, as well as to the improvements in the treatment of energy deposition, now implemented in EXOCLASS. Accounting for the halo (light shaded regions) does not lead to significant differences unless $f_{PBH}$ is sufficiently small, i.e., there is sufficient material for growing a sizable DM halo. The threshold to see significant improvements is $f_{PBH} \lesssim 0.01$ for the disk accretion case but already at $f_{PBH} \lesssim 0.2$ for the spherical accretion case. For the latter case, the steep improvement of the bound around $M \sim 30 M_\odot$ in the presence of a halo is only indicative, since for $f_{PBH} \gtrsim 0.01$ a non-negligible fraction of the DM can be gravitationally bound to two or more PBH, and the radial profile derived in the isolated-PBH approximation breaks down [37]. At higher masses, sensitive to lower $f_{PBH}$, the approximation is however robust: The formation of a DM halo around the PBH can strikingly improve the bound by up to ~2 orders of magnitude in the covered mass range. Our results also show that the bounds eventually flatten when

<footnote n="4">⁴https://github.com/brinckmann/montepython_public.</footnote>

<footnote n="5">⁵https://github.com/dscolnic/Pantheon.</footnote>
<figure><img src="image_11.png" /><figcaption>FIG. 7. Bounds on the abundance of PBH assuming disk accretion (light shade) or spherical accretion (dark shade), accounting for the formation of a DM halo. The most stringent bounds in the same mass regions are also shown: Second order GW ones from Ref. [3], Icarus ones from Ref. [60], LIGO ones according to Ref. [5], BBN bounds from Ref. [75], spectral CMB distortions from Ref. [76]. The arrow indicates that for masses $M \gtrsim 200 M_{\odot}$, PBH can in principle grow in mass up to $10^9 M_{\odot}$ by $z = 7.5$ by accreting baryons at Eddington luminosity with $\epsilon = 0.1$.</figcaption></figure>

$M \gtrsim 10^4 M_{\odot}$. This is a consequence of the accretion attaining the Eddington limit for longer and longer periods of time, thus converging to Eq. (25). As previously argued, in this range the bounds become shaky since the working hypotheses break-down.

It is worth commenting on the relative strength of the derived bounds with other existing ones, with the most stringent ones reported in Fig. 7. Since curvature perturbations couple to tensor perturbations at second order, PBH below the solar-mass scale are associated to GWs generated in conjunction with their formation, falling in the frequency probed by pulsar timing arrays. The nonobservation of a stochastic signal in the nHz range sets tight bounds [3,55]. Galactic microlensing constraints [56–58], roughly excluding $f_{\text{PBH}} \gtrsim \mathcal{O}(0.1)$, also apply. Other “direct” bounds come from the nonobservation of mergers by LIGO/Virgo [59]. At few solar masses, leading constraints come from caustic crossing events in giant arcs (produced by stars embedded in high magnification regions due to a Galaxy cluster) [60], but one may expect similar or tighter constraints from the extrapolation of the analysis of Ref. [59] to higher masses.

In the 10–100 $M_{\odot}$ range, the binary coalescence rate inferred by LIGO/Virgo is estimated to yield bounds at a level between $10^{-3}$ and $10^{-2}$ [4,5], a result whose robustness to a number of effects has been checked in Refs. [6,61,62]. Note however that, according to Refs. [63,64], accounting for binary disruption can relax these limits to some extent.

Other constraints at $M \sim \mathcal{O}(10) M_{\odot}$ roughly in the ballpark of $f_{\text{PBH}} \lesssim \mathcal{O}(0.1)$ come from the nonobservation of a stochastic gravitational wave (GW) background (due to the mergers of PBH binaries at high-z, in the matter dominated era) [65], quasar microlensing [66], lensing of type-Ia supernovae [67], or the orbital dynamics of halo wide binaries [68]. When approaching the ~100 $M_{\odot}$ scale, radio and x-ray observations of the Milky Way [69], the half-light radius of dwarf galaxies [70,71], or the stellar distribution of dwarf galaxies [72] take over as more and more stringent bounds. Basically, the CMB constraints surpass all these at masses $M \gtrsim 20–50 M_{\odot}$, and remain the dominant constraint until at least $10^{3.5} M_{\odot}$, when they become comparable to (or slightly better than) BBN ones [73–75], before being definitely surpassed by CMB spectral distortions (see Refs. [76,77] and Refs. therein) at $M \gtrsim 10^{4.5} M_{\odot}$. Needless to say, since different constraints are derived in different systems and are affected by different systematics, the existence of multiple arguments excluding some parameter space strengthens their credibility and robustness. In particular, over all the stellar mass range, PBH as the totality of DM are excluded by at least two arguments, often more. It is worth noting that even in the most conservative case that we consider, the CMB now provides an independent argument excluding PBH of $M \gtrsim 15 M_{\odot}$ as the totality of DM; around $30 M_{\odot}$, no more than $f_{\text{PBH}} \sim \mathcal{O}(0.1)$ is allowed. On the other hand, once accounting for uncertainties, the CMB is not capable of disproving a
primordial origin of the bulk of LIGO/Virgo merger events, estimated according to Refs. [4,5]. But it is interesting to see how the disk accretion scenario is in strong tension with this interpretation, thus providing a phenomenological motivation for reducing accretion-related uncertainties.

## V. IMPLICATIONS FOR SMBH

In this section, we discuss the implications, if any, of the CMB bounds previously derived for the still mysterious genesis of SMBH.

### A. The “SMBH problem”

Supermassive black holes (SMBH)—loosely defined as BH whose mass exceeds $10^5 M_\odot$—are believed to sit at the center of almost all galaxies, and their integrated accretion disk emission almost saturates the cosmological x-ray background, see for instance the review [78]. Additionally, SMBH—including a few very massive ones with $M \gtrsim 10^9 M_\odot$—have been observed at high redshift $z \gtrsim 6$, a fact which seriously constraints their formation mechanism [79].

One may wonder if these SMBH may have developed from lighter black holes via accretion phenomena, a process that is known to be at play at $0 \le z \le 6$. Since SMBH with $M \gtrsim 10^9 M_\odot$ were already in place by the time the universe was 1 billion years old, they raise some concerns. The standard argument goes as follows:

Based on Eqs. (3)–(5), if $\epsilon$ can be estimated, $L_E$ can be linked to what is considered an optimistic benchmark accretion value, since we expect the maximal$^6$ mass accretion rate to be the complement of mass-accretion rate converted into radiation at the Eddington limit, i.e.,

$$ \dot{M} \lesssim \dot{M}_E \equiv \frac{1 - \epsilon L_E}{\epsilon c^2} \propto \frac{1 - \epsilon}{\epsilon} M. \quad (26) $$

Note that the Eddington mass accretion rate $\dot{M}_E$ depends on the unknown quantity $\epsilon$ and is mathematically unbounded from above, when $\epsilon \to 0$, while it can be arbitrarily small, when $\epsilon \to 1$. However, the largest known values of $\epsilon$ that can be attained are $\epsilon \simeq 0.42$ for a maximally rotating Kerr BH, so that the following inequality

$$ \dot{M}_E \gtrsim \dot{M}_{E,\min} \approx \frac{L_E}{c^2} = \frac{M}{\tau_E} \quad (27) $$

seems to hold for all known systems.

In the literature, the mass accretion rate of a BH of mass $M_i$ at formation time $t_i$ is thus limited as

$$ M(t) \lesssim M_i \times \exp\left(\frac{1 - \epsilon t - t_i}{\epsilon \tau_E}\right), \quad (28) $$

where $\tau_E$ was given in Eq. (4). For a typical benchmark value $\epsilon \approx 0.1$, this implies an e-fold time for the BH of the order of 0.04 Gyr or a maximum estimated growth by accretion of 17 e folds between the epoch of first stars dying at $z \simeq 15$ and observations at $z \simeq 6$. Thus, there is barely the time for

a stellar BH of mass $O(100) M_\odot$ to grow to the size of the heaviest SMBH at $z \gtrsim 6$.

The above argument hides the loophole that, if the efficiency of conversion of accretion into luminosity drops, the actual mass growth may be significantly larger. Radiatively inefficient accretion at very high inflow rates has been discussed in the past as a way to evade the above argument and suggested by a number of simulations and theoretical arguments, see for instance Refs. [80,81]. In general, accretion in this regime may become nonstationary, see for instance the considerations in Ref. [29]. There are also some mechanisms alternative to the accretion mechanism onto stellar-mass black holes to overcome the difficulty in forming SMBH, such as invoking runaway mergers in dense clusters [82,83], or direct collapse of BH from a gas cloud [84], see Ref. [85] for a review.

Nonetheless, the difficulty of achieving the required conditions has led several authors (see for instance Ref. [86]) to speculate that SMBHs or rather of their seeds may have a primordial origin, possibly linked with the origin of galactic structures, an old idea recently reviewed in Ref. [77].

### B. Primordial SMBH?

Similarly to the Milky Way halo mass fraction in the SGR A* black hole, SMBHs currently account for about $10^{-5}$ of the DM mass density in the universe (e.g., Ref. [87], see also the gray band in Fig. 5 of Ref. [88]). However, as we just reviewed, it is known that SMBH undergo significant growth with time. In fact, at $z \simeq 6$ the overall mass density into SMBH above $10^6 M_\odot$ was only about a factor $10^{-3.5}$ of the current value, such that the difference between these figures must be accounted for via mergers, accretion, and newly formed objects. A more quantitative description of the high-redshift SMBH mass function can be given in terms of the so-called Schechter function,

$$ \frac{dn_{\text{BH}}}{d \log_{10} m} = m \ln 10 \frac{dn_{\text{BH}}}{dm} = \kappa m^{\alpha} e^{-m}, \quad (29) $$

with inferred values at $z=6$ of $\kappa = 1.23 \times 10^{-8} \text{Mpc}^{-3}$, $\alpha = -1.03$ and $m \equiv M/M_*$, with $M_* = 2.24 \times 10^9 M_\odot$ (see Ref. [89] or equivalently Fig. 2 in Ref. [90]). This is consistent with the inferred co-moving density $>1.1 \times 10^{-9} \text{Mpc}^{-3}$ above $10^9 M_\odot$ between $z=6.44$ and $z=7.44$ reported in Ref. [89]. If translated in terms of the DM fraction, Eq. (29) yields about $96 M_\odot \text{Mpc}^{-3}$ above $10^6 M_\odot$, equivalent to a fraction of the DM abundance in SMBH above $10^6 M_\odot$ of $f_{\text{PBH}} \simeq 2.9 \times 10^{-9}$. Thus, even under the extreme case of Eq. (25), the CMB angular power spectra do not exclude a primordial origin hypothesis for the SMBHs.

Are there counterarguments to this? An apparent theoretical difficulty is that one expects a direct formation of SMBHs to happen after the weak reaction freeze-out, since the horizon mass scales roughly as $M_H \simeq 10^5 (t/s) M_\odot$. However, having a very tiny fraction of matter in the form of primordial SMBHs at BBN times, or even somewhat after BBN, is not obviously excluded and only limited by theoretical creativity. A more serious concern is that, if SMBHs form from (quasi)Gaussian fluctuations, the mass $6 \times 10^4 M_\odot \lesssim M \lesssim 5 \times 10^{13} M_\odot$ is subject to tight constraints coming from

<footnote>$^6$This is probably unrealistic, since it excludes any sizable outflows.</footnote>
CMB spectral distortions [76]. No cosmologically relevant abundance is allowed in this range unless the PBH form out of highly non-Gaussian tail fluctuations [91–93].

In summary, this discussion leaves two possible (primordial) scenarios:

(1) *Primordial SMBH hypothesis*: SMBHs with a mass function similar to the inferred one, Eq. (29), are directly of primordial origin. This requires PBHs to form under rather peculiar highly non-Gaussian conditions in order to fulfill CMB spectral constraints. Also, the bulk of the SMBH population is required to undergo negligible mass growth in the period before reionization (dark ages) not to overshoot the inferred mass function, a condition which appears rather challenging to fulfill and puzzling if compared to the $10^{3.5}$ growth observationally deduced between $z \simeq 6$ and today.

(2) *Primordial SMBH seed hypothesis*: If PBHs form with masses $M \lesssim 10^4 M_\odot$ and $f_{PBH} \lesssim 10^{-9}$, they are consistent with all present bounds, also for initial Gaussian conditions. The presence of SMBHs with masses $M \gtrsim 10^9 M_\odot$ at redshift $z \gtrsim 6$ requires then a mass growth by a factor of at least $\simeq 10^5$ or $11.5$ e folds, which however appears “easily” achieved within the naive theory sketched above. A sufficient growth can be attained for $M_{PBH} \simeq 200 M_\odot$ if accreting at Eddington limit and with $\epsilon = 0.1$. A similar mass growth process is needed anyway also in astrophysical scenarios, which are usually further constrained by a later formation epoch (typically $z \lesssim 15$) and, at least for popIII (as opposed to direct collapse) scenarios, by a seed mass not usually exceeding $10^2 M_\odot$.

In the second scenario, note that although the PBH get dressed with a DM halo which is one to two orders of magnitude its mass, the DM halo does not necessarily contribute to the “inferred” SMBH mass, since the latter is typically deduced from the properties of the inner accretion disk emission (x-ray or radio data). The DM halo rather creates favorable conditions to boost the baryonic accretion. Finally, since PBH of the relevant masses can be formed before the earliest cosmological timescales probed (weak reaction freeze-out), it is conceivable that theoretical models are more easily constructed in this scenario.

## VI. CONCLUSIONS

Stellar mass or heavier primordial black holes may have interesting cosmological and astrophysical consequences, even if they only constitute a small fraction of the overall amount of dark matter. However, in such a situation their interplay with the remaining fraction of DM may have peculiar consequences, as it has been noted in several instances (see, e.g., Refs. [40,94–96] for implications for DM models). Here we have revisited the impact that the growth of DM halos around PBH has on CMB anisotropy constraints. We have elucidated the effects of these DM halos on the baryonic accretion thanks to both simple semianalytical models and dedicated numerical simulations, and derived state-of-the-art cosmological bounds. These limits are the leading ones in the window between a few tens solar mass (where a number of astrophysical bounds exist, typically at the $f_{PBH} \sim 10^{-3} \rightarrow 10^{-1}$ level) and the tight cosmological bounds from CMB spectral distortions at masses $M \gtrsim 10^4 M_\odot$ (for a summary, see Fig. 7). The CMB anisotropy bounds reach very deep

down in the $f_{PBH}$ range, becoming as stringent as $\lesssim 10^{-8}$ at the highest masses of applicability, both under the hypotheses of spherical or disk accretion. As the largest uncertainty comes from the accretion model, these bounds could be further refined via dedicated hydrodynamical simulations in a cosmological setting, which would be particularly useful to explore the role of nonstationarity in SMBH accretion and extend the bounds beyond $10^4 M_\odot$, but also to assess if the CMB bounds exclude a PBH origin of the LIGO/Virgo merger events. In particular, analytical calculations predict supersonic motion of the baryons at cosmological scales [25]; starting from those initial conditions, it would be interesting to understand the dynamics of the baryonic gas actually accreted onto the PBH and their DM halos, subject to shocks and dissipative effects.

Despite still existing uncertainties, we argued that CMB bounds do not prevent a primordial origin for the very heavy supermassive black holes observed already at $z > 6$. In particular, we find that a scenario where PBHs with $M \gtrsim 10^3 M_\odot$ act as the seeds of the SMBH easily fulfills all the known constraints. A prediction of such a scenario is that rather massive black holes are already around and accreting at $z \sim 30$. Qualitatively, we can thus expect interesting implications for the dark ages, such as nonstandard 21 cm and reionization epoch (for first explorations, see, e.g., Refs. [97,98]), which will surpass next generations CMB observations in constraining e.m. energy injection in the dark ages and therefore definitely deserve additional dedicated studies.

## ACKNOWLEDGMENTS

We acknowledge valuable discussions with Yacine Ali-Haïmoud. We thank Nagisa Hiroshima for collaboration on an early stage of the project and Nagisa Hiroshima and Guillermo Ballesteros for comments on the paper. This work is partially supported by the project *Multimessenger avenues in gravitational waves* under the program “Initiatives de Recherche Stratégique” - IDEX Univ. Grenoble-Alpes (P.D.S.), by JSPS KAKENHI Grant No. JP17H01131 (K.K.), MEXT Grant-in-Aid for Scientific Research on Innovative Areas JP15H05889. JP18H04594, JP19H05114 (K.K.), and by WPI, MEXT, Japan (K.K.), the Toshiko Yuasa France-Japan Particle Physics Laboratory “TYL-FJPPL” (V.P.), by the National Science Foundation under Grant No. 1820861 (D.I.), by the NYU IT High Performance Computing resources, services, and staff expertise (D.I.).

## APPENDIX A: SIZE OF THE REGION OF INFLUENCE OF A PBH

When computing the constraints on PBH accretion from the CMB, it is assumed that the *mean* ionized fraction is affected by the energetic radiation emitted by the PBH. We wish to check whether this assumption is valid and that one should not consider the influence of PBHs as a local *perturbation* to the homogeneous ionization fraction. It was checked in Ref. [8] that a photon can escape the very dense environment close to the PBHs; we therefore focus on the far-away region, where densities are cosmological. As an estimate of the typical size of the region influenced by PBHs, one can calculate the mean free path of emitted keV photons,
which deposit their energy mostly via Thomson scattering with nonrelativistic electrons, i.e.,

$$ \lambda_T \equiv (n_e \sigma_T)^{-1} \approx 2 \times 10^4 x_e^{-1} \left(\frac{1000}{1+z}\right)^3 \text{ pc.} \quad (A1) $$

This is to be compared with the typical distance between PBHs

$$ \bar{r} = \left( \frac{3M}{4\pi \rho_{\text{PBH}}} \right)^{1/3} \approx 2 \times 10^{-1} \left( \frac{M}{f_{\text{PBH}} M_{\odot}} \right)^{1/3} \frac{1000}{1+z} \text{ pc.} \quad (A2) $$

Hence, one finds that PBHs can influence (and ionize) all of the region separating them from another PBH as long as

$$ f_{\text{PBH}} > 10^{-15} x_e^3 \frac{M_{\odot}}{M}, \quad (A3) $$

which is always satisfied given the range of PBH masses and fractions considered in this work.

Additionally, one can calculate the number of PBHs per patch of the CMB sky as seen by Planck at the smallest multipoles, in order to make sure that these patches average out the contribution of many PBHs to the ionization fraction. In the flat sky approximation, valid at these small scales (large multipole $\ell$), the typical comoving size of the patch $\lambda \equiv d_A(z_{\text{dec}})/\ell \sim 10^{10}$ pc/$\ell$ contains

$$ N_{\text{PBH}} \sim \frac{a\lambda}{\bar{r}} \approx 5 \times 10^7 \ell^{-1} \left( \frac{f_{\text{PBH}} M_{\odot}}{M} \right)^{1/3} . \quad (A4) $$

At the highest multipole seen by Planck, $\ell \sim 2000$, there are therefore more than one PBH per patch up to masses $M \sim 10^4 M_{\odot}$ for a fraction $f_{\text{PBH}} = 10^{-9}$.

## APPENDIX B: DYNAMICAL FRICTION OF A MASSIVE PBH

Here we want to assess the cosmological relevance of the dynamical friction that a PBH experiences in the cosmological

baryonic gas of density $\rho_b$ moving relatively to it at supersonic speed $v \simeq v_L$ (the drag being suppressed at subsonic velocities [99]). A body of mass $M$ experiences the energy loss rate per unit distance with the same form as for a collisionless medium [100]

$$ -\frac{dE}{dx} = 4\pi \rho_b \frac{(GM)^2}{v^2} \ln \Lambda, \quad (B1) $$

where $\ln \Lambda \simeq \mathcal{O}(10)$ is the Coulomb logarithm, depending on the ratio of the largest and smallest linear scales involved in the process. The energy loss timescale is thus

$$ \tau_{\text{loss}} = \frac{M v^2 / 2}{-v dE/dx} = \frac{v^3}{3GM \ln \Lambda} \frac{3}{8\pi G \rho_b}. \quad (B2) $$

In order to assess its cosmological relevance, let us compute the ratio of the above scale with the Hubble time $\sim H^{-1}$. Plugging $v \simeq v_L(z)$ from Eq. (11) and using the first Friedmann equation in Eq. (B2), we obtain

$$ \tau_{\text{loss}}(z)H(z) \approx 1.8 \times 10^4 \frac{M_{\odot}}{M} \left(\frac{1+z}{100}\right)^{3/2} \frac{10}{\ln \Lambda}. \quad (B3) $$

This indicates that for stellar-mass PBH this effect is sub-leading at the cosmological epoch of interest, but that at the highest PBH masses of interest ($M \gtrsim 10^4 M_{\odot}$) the baryon-PBH motion may eventually settle to subsonic during the dark ages. Hence, for studying intermediate mass PBH accretion in pristine halos, as well as their impact on dark ages observables, this effect should be properly accounted for.

[1] G. Ballesteros, J. Rey, M. Taoso, and A. Urbano, arXiv:2001.08220.

[2] A. Kusenko, M. Sasaki, S. Sugiyama, M. Takada, V. Takhistov, and E. Vitagliano, arXiv:2001.09160.

[3] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, arXiv:2002.12778.

[4] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Phys. Rev. Lett. 117, 061101 (2016).

[5] Y. Ali-Haïmoud, E. D. Kovetz, and M. Kamionkowski, Phys. Rev. D 96, 123523 (2017).

[6] B. J. Kavanagh, D. Gaggero, and G. Bertone, Phys. Rev. D 98, 023536 (2018).

[7] M. Ricotti, J. P. Ostriker, and K. J. Mack, Astrophys. J. 680, 829 (2008).

[8] Y. Ali-Haïmoud and M. Kamionkowski, Phys. Rev. D 95, 043534 (2017).

[9] B. Horowitz, arXiv:1612.07264.

[10] L. Chen, Q.-G. Huang, and K. Wang, JCAP 12 (2016) 044.

[11] V. Poulin, P. D. Serpico, F. Calore, S. Clesse, and K. Kohri, Phys. Rev. D 96, 083524 (2017).

[12] J. L. Bernal, N. Bellomo, A. Raccanelli, and L. Verde, JCAP 10 (2017) 052.

[13] E. Bañados, B. P. Venemans, C. Mazzucchelli, E. P. Farina, F. Walter, F. Wang, R. Decarli, D. Stern, X. Fan, F. B. Davies, J. F. Hennawi, R. A. Simcoe, and M. L. Turner, H.-W. Rix, J. Yang, D. D. Kelson, G. C. Rudie, and J. M. Winters, Nature (London) 553, 473 (2018).

[14] K. Inayoshi, E. Visbal, and Z. Haiman, Annu. Rev. Astron. Astrophys. 58 (2020).

[15] P. Stöcker, M. Krämer, J. Lesgourgues, and V. Poulin, JCAP 03 (2018) 018.

[16] H. Liu, G. W. Ridgway, and T. R. Slatyer, Phys. Rev. D 101, 023530 (2020).

[17] T. R. Slatyer, N. Padmanabhan, and D. P. Finkbeiner, Phys. Rev. D 80, 043526 (2009).

[18] T. R. Slatyer, Phys. Rev. D 93, 023521 (2016).

[19] S. L. Shapiro, Astrophys. J. 180, 531 (1973).

[20] S. L. Shapiro, Astrophys. J. 185, 69 (1973).

[21] H. Bondi, Mon. Not. Roy. Astron. Soc. 112, 195 (1952).
[22] F. Hoyle and R. A. Lyttleton, *Math. Proc. Cambridge Philos. Soc.* **35**, 405 (1939).

[23] F. Hoyle and R. A. Lyttleton, *Math. Proc. Cambridge Philos. Soc.* **36**, 325 (1940).

[24] F. Hoyle and R. A. Lyttleton, *Math. Proc. Cambridge Philos. Soc.* **36**, 424 (1940).

[25] D. Tseliakhovich and C. Hirata, *Phys. Rev. D* **82**, 083520 (2010).

[26] G. Hütsi, M. Raidal, and H. Veermäe, *Phys. Rev. D* **100**, 083016 (2019).

[27] K. Park and M. Ricotti, *Astrophys. J.* **767**, 163 (2013).

[28] M. Ricotti (2017), URL icc.ub.edu/congress/ICCBSchool2017/documents/Slides/Ricotti/Lecture4.pdf.

[29] M. Ricotti, *Astrophys. J.* **662**, 53 (2007).

[30] E. Takeo, K. Inayoshi, K. Ohsuga, H. R. Takahashi, and S. Mineshige, *Mon. Not. Roy. Astron. Soc.* **488**, 2689 (2019).

[31] *Formation of the First Black Holes*, edited by M. Latif and D. R. G. Schleicher (World Scientific Publishing Company, 2018), pp. 195–228.

[32] E. Bertschinger, *Astrophys. J. Suppl.* **58**, 39 (1985).

[33] K. J. Mack, J. P. Ostriker, and M. Ricotti, *Astrophys. J.* **665**, 1277 (2007).

[34] V. S. Berezinsky, V. I. Dokuchaev, and Yu. N. Eroshenko, JCAP 11 (2013) 059.

[35] K. Park, M. Ricotti, P. Natarajan, T. Bogdanovic, and J. H. Wise, *Astrophys. J.* **818**, 184 (2016).

[36] J. Harnois-Deraps, U.-L. Pen, I. T. Iliev, H. Merz, J. D. Emberson, and V. Desjacques, *Mon. Not. Roy. Astron. Soc.* **436**, 540 (2013).

[37] D. Inman and Y. Ali-Haïmoud, *Phys. Rev. D* **100**, 083528 (2019).

[38] R. I. Epstein, *Mon. Not. Roy. Astron. Soc.* **205**, 207 (1983).

[39] R. K. Sheth, *Mon. Not. Roy. Astron. Soc.* **276**, 796 (1995).

[40] J. Adamek, C. T. Byrnes, M. Gosenca, and S. Hotchkiss, *Phys. Rev. D* **100**, 023506 (2019).

[41] P. Ullio, H. S. Zhao, and M. Kamionkowski, *Phys. Rev. D* **64**, 043504 (2001).

[42] B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, JCAP 02 (2013) 001.

[43] T. Brinckmann and J. Lesgourgues, *Phys. Dark Univ.* **24**, 100260 (2019).

[44] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak et al. (Planck Collaboration), arXiv:1807.06209.

[45] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, K. Benabed et al. (Planck Collaboration), arXiv:1907.12875.

[46] F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, *Mon. Not. Roy. Astron. Soc.* **416**, 3017 (2011).

[47] A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, *Mon. Not. Roy. Astron. Soc.* **449**, 835 (2015).

[48] S. Alam, M. Ata, S. Bailey, F. Beutler, D. Bizyaev, J. A. Blazek, A. S. Bolton, J. R. Brownstein, A. Burden, C.-H. Chuang et al. (BOSS), *Mon. Not. Roy. Astron. Soc.* **470**, 2617 (2017).

[49] D. M. Scolnic, D. O. Jones, A. Rest, Y. C. Pan, R. Chornock, R. J. Foley, M. E. Huber, R. Kessler, G. Narayan, A. G. Riess et al., *Astrophys. J.* **859**, 101 (2018).

[50] A. Lewis, *Phys. Rev. D* **87**, 103529 (2013).

[51] A. Gelman and D. B. Rubin, *Statist. Sci.* **7**, 457 (1992).

[52] C. T. Byrnes, P. S. Cole, and S. P. Patil, JCAP 06 (2019) 028.

[53] F. Kühnel and K. Freese, *Phys. Rev. D* **95**, 083508 (2017).

[54] B. Carr, M. Raidal, T. Tenkanen, V. Vaskonen, and H. Veermäe, *Phys. Rev. D* **96**, 023514 (2017).

[55] Z.-C. Chen, C. Yuan, and Q.-G. Huang, arXiv:1910.12239.

[56] P. Tisserand et al. (EROS-2), *Astron. Astrophys.* **469**, 387 (2007).

[57] L. Wyrzykowski et al., *Mon. Not. Roy. Astron. Soc.* **416**, 2949 (2011).

[58] A. M. Green, *Phys. Rev. D* **94**, 063530 (2016).

[59] B. P. Abbott et al. (LIGO Scientific, Virgo Collaboration), *Phys. Rev. Lett.* **121**, 231103 (2018).

[60] M. Oguri, J. M. Diego, N. Kaiser, P. L. Kelly, and T. Broadhurst, *Phys. Rev. D* **97**, 023518 (2018).

[61] M. Raidal, V. Vaskonen, and H. Veermäe, JCAP 09 (2017) 037.

[62] G. Ballesteros, P. D. Serpico, and M. Taoso, JCAP 10 (2018) 043.

[63] M. Raidal, C. Spethmann, V. Vaskonen, and H. Veermäe, JCAP 02 (2019) 018.

[64] V. Vaskonen and H. Veermäe, *Phys. Rev. D* **101**, 043015 (2020).

[65] S. Wang, Y.-F. Wang, Q.-G. Huang, and T. G. F. Li, *Phys. Rev. Lett.* **120**, 191102 (2018).

[66] E. Mediavilla, J. Jiménez-Vicente, J. A. Muñoz, H. Vives-Arias, and J. Calderón-Infante, *Astrophys. J.* **836**, L18 (2017).

[67] M. Zumalacarregui and U. Seljak, *Phys. Rev. Lett.* **121**, 141101 (2018).

[68] M. A. Monroy-Rodríguez and C. Allen, *Astrophys. J.* **790**, 159 (2014).

[69] J. Manshanden, D. Gaggero, G. Bertone, R. M. T. Connors, and M. Ricotti, JCAP 06 (2019) 026.

[70] T. D. Brandt, *Astrophys. J.* **824**, L31 (2016).

[71] T. S. Li et al. (DES), *Astrophys. J.* **838**, 8 (2017).

[72] S. M. Koushiappas and A. Loeb, *Phys. Rev. Lett.* **119**, 041102 (2017).

[73] T. Nakama, T. Suyama, and J. Yokoyama, *Phys. Rev. Lett.* **113**, 061302 (2014).

[74] D. Jeong, J. Pradler, J. Chluba, and M. Kamionkowski, *Phys. Rev. Lett.* **113**, 061301 (2014).

[75] K. Inomata, M. Kawasaki, and Y. Tada, *Phys. Rev. D* **94**, 043527 (2016).

[76] K. Kohri, T. Nakama, and T. Suyama, *Phys. Rev. D* **90**, 083514 (2014).

[77] B. Carr and J. Silk, *Mon. Not. Roy. Astron. Soc.* **478**, 3756 (2018).

[78] A. Merloni, *Lect. Notes Phys.* **905**, 101 (2016).

[79] M. Volonteri, *Astro. Astrophys. Rev.* **18**, 279 (2010).

[80] G. S. Novak, arXiv:1310.3833.

[81] M. Volonteri, J. Silk, and G. Dubus, *Astrophys. J.* **804**, 148 (2015).

[82] G. D. Quinlan and S. L. Shapiro, *Astrophys. J.* **356**, 483 (1990).
[83] S. F. Portegies Zwart, H. Baumgardt, P. Hut, J. Makino, and S. L. W. McMillan, *Nature (London)* **428**, 724 (2004).

[84] A. Loeb and F. A. Rasio, *Astrophys. J.* **432**, 52 (1994).

[85] M. C. Begelman, M. Volonteri, and M. J. Rees, *Mon. Not. Roy. Astron. Soc.* **370**, 289 (2006).

[86] N. Düchting, *Phys. Rev. D* **70**, 064015 (2004).

[87] Q.-J. Yu and S. Tremaine, *Mon. Not. Roy. Astron. Soc.* **335**, 965 (2002).

[88] F. Shankar, D. H. Weinberg, and J. Miralda-Escude, *Astrophys. J.* **690**, 20 (2009).

[89] C. J. Willott, L. Albert, D. Arzoumanian, J. Bergeron, D. Crampton, P. Delorme, J. B. Hutchings, A. Omont, C. Reylé, and D. Schade, *Astronomi. J.* **140**, 546 (2010).

[90] M. Volonteri and D. P. Stark, *MNRAS* **417**, 2085 (2011).

[91] T. Nakama, B. Carr, and J. Silk, *Phys. Rev. D* **97**, 043525 (2018).

[92] T. Nakama, T. Suyama, and J. Yokoyama, *Phys. Rev. D* **94**, 103522 (2016).

[93] J. Garcia-Bellido, M. Peloso, and C. Unal, JCAP 09 (2017) 013.

[94] B. C. Lacki and J. F. Beacom, *Astrophys. J.* **720**, L67 (2010).

[95] G. Bertone, A. M. Coogan, D. Gaggero, B. J. Kavanagh, and C. Weniger, *Phys. Rev. D* **100**, 123013 (2019).

[96] M. P. Hertzberg, E. D. Schiappacasse, and T. T. Yanagida, arXiv:2001.07476.

[97] A. Hektor, G. Hütsi, L. Marzola, M. Raidal, V. Vaskonen, and H. Veermäe, *Phys. Rev. D* **98**, 023503 (2018).

[98] O. Mena, S. Palomares-Ruiz, P. Villanueva-Domingo, and S. J. Witte, *Phys. Rev. D* **100**, 043540 (2019).

[99] E. C. Ostriker, *Astrophys. J.* **513**, 252 (1999).

[100] S. Chandrasekhar, *Astrophys. J.* **97**, 255 (1943).