CMB bounds on primordial black holes including dark matter halo accretion

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 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_{\rm 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.


I. INTRODUCTION
Despite several decades of direct, indirect or collider searches, we still ignore the nature of the dark 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 < ∼ M < ∼ 10 −11 M (see e.g. [1]). 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 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 PBH ∼ few × 10 −3 to the DM [2][3][4].
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 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 [5] to the amended and more recent calculations in Ref. [6][7][8][9][10].
In this article, we extend our previous calculation [9] to account for the accretion of the dominant, non-PBH * Electronic address: serpico@lapth.cnrs.fr DM particles onto PBH, thus altering the baryonic accretion and eventually the PBH luminosity constrained by the CMB. 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 < ∼ M < ∼ 10 4 M ) may provide the seeds leading to the super-massive black holes (SMBHs) observed at redshift z > ∼ 6 (with the current record-holder of 10 8.9 M at z = 7.54 [11]) whose origin constitutes a long-standing mystery, see e.g. [12]. To the best of our knowledge, nobody has assessed the viability of this hypothesis in the light of current CMB anisotropy bounds.
This article 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 in 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 accreting around PBH, treating them both via semi-analytical toy models and dedicated numerical simulations. Our CMB bounds are reported in Sec. IV, providing an upgrade to the the bounds presented in Ref. [9] in the light of the last Planck data release and the updated treatment of the energy release via the new ExoCLASS package [13]. Also, we present an extension of the bounds up to masses of M few × 10 4 M , beyond which some working hypotheses break down. Sec. 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: where L 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 < ∼ 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 1 from Ref. [15,16] as implemented in the ExoCLASS package [13]. As in Ref. [9], we assume that the energy-differential spectrum of L acc , L ω , is dominated by Bremsstrahlung emission (see [6,17,18]), 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 where T s ∼ O(m e ) (we use 200 keV in the following, like in Ref. [9]) and |a| < ∼ 0.5, again like in Ref. [9] (a = 0 was used in Ref. [6]). Before addressing the question of the effective PBH mass evolution, let us review typical benchmarks for L acc . Also, note that if the time-dependence in M (t) is not negligible, then f PBH may become timedependent 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 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 1 More recent tools have been developed for the computation of the energy deposited in the cosmological plasma during the cosmic dark ages [14], but it has been shown that the effect of energy injection onto the CMB bounds is accurately described by the tabulated transfer functions.
spherical system, and can be simply computed as with σ T the Thomson cross-section, m p the proton mass, and µ a chemical composition dependent parameter, equal to unity for pure hydrogen. Note that L E is a quantity linear in the accreting object mass, and allows one to introduce a fundamental scale (independent of the accreting object) known as Salpeter timescale, Under the rather extreme hypothesis that L 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 . Also, if the PBH mass accretion timescale is of the same order of magnitude of eq. (4), then the PBH mass can be considered constant and equal to its initial value for z 10, which is indeed the regime relevant for deriving CMB bounds. This remains true a fortiori for less extreme accretion rates.

B. More realistic accretion scenarios
The function L acc in Eq. (1) is usually parameterized in terms of the two following quantities: i)Ṁ , the matter accreted per unit time onto the PBH.
ii) , the overall efficiency of conversion of accreted matter into radiation, in terms of which one writes: ConcerningṀ , an analytical theory exists in two limiting cases, both applying to a homogeneous gas of mass density ρ ∞ : The hypothesis of a stationary, spherical symmetric accretion of a body at rest [19] (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 rel in the gas [20][21][22] (Hoyle-Lyttleton). These limiting cases justify the following "Bondi-Hoyle-Lyttleton" (BHL) pa-rameterizationṀ where c s is the speed of sound in the homogeneous matter of density ρ ∞ , and λ is a dimensionless coefficient dependent upon environmental parameters, of O(1) for v rel c s (Hoyle-Lyttleton), and a calculable function assuming values of O(0.1-1) (see e.g. [6]) in the limit c s v rel (Bondi). For more realistic situations, eq. (6) is often used, but with λ now intended as an adjustable parameter or function fitted e.g. to simulation results. We can also define a "generalized" Bondi radius, in terms of which Eq. (6) writeṡ and which, just like Eq. (8), reduces to Bondi's results for v eff → c s , while smoothly interpolating to the Hoyle-Lyttleton regime for larger velocities. Concerning , it can be computed in spherical symmetry under some assumptions for the radiative processes, the most up-to-date treatment being provided in Ref. [6]. 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 λ and the same prescriptions used in Ref. [9], 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 non-perturbative) theory [23], 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 [9]. 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 ∼ 2 − 5) and, as argued for instance in Ref. [24,25], the accretion should then be disk-like 2 . In this case, our fiducial parameterization closely follows [9] where, at z < ∼ 1000, we adopt [6,23] c s 6 km s ii) If, on the other hand, the relevant small-scale motions are subsonic ( i.e. M < ∼ 1), for the low values of f 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 eff 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 [6].

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: As already argued in Ref. [26], this leads to the requirement M < ∼ few × 10 4 M , which is the upper limit for which we present our results.
ii) How to deal with super-Eddington accretion, i.e. L 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 multi-dimensional effects and realistic disk radiation spectra (for a recent study, see [27]). 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 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 acc < ∼ L E . While we will comment again on this regime in Sec. V, we address the reader to reviews such as Ref. [28] 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 acc > L E .
Interestingly enough, super-Eddington accretion is attained in our formalism for M > ∼ 10 4 M , 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 App. 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 [29][30][31] 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 essentially constant down to the redshifts of interest for CMB bounds. Hence, we can safely consider f 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 [32]: where r B,eff , the effective Bondi radius, is the unknown, M PBH is the initial PBH mass and Φ h the (timedependent) gravitational potential associated to the DM halo. Our treatment of the problem consists in adopting Eq. (9), but with r B replaced by r B,eff , 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, where ρ 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 ρ PBH = 3 M PBH /(4π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 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. [30]. Under the above-mentioned approximations, we find good agreement with the results reported in Ref. [30], namely: • A time evolution given by • A density profile proportional to ∝ 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. Eq. (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 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 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. [26], Sec. 4, that other scaling solutions like the ones described in the seminal paper by Bertschinger [29] 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 [33] modified to include PBHs as a separate particle species co-evolved with a generic collisionless DM candidate [34]. We opt for homogeneous DM initial conditions at a = 10 −6 and select cosmological parameters consistent with Planck: Ω c = 0.26, Ω b = 0.05 and z eq = 3374, with Ω PBH = f PBH Ω c . Baryons are not evolved and instead assumed to be homogeneous inducing errors of order Ω b /Ω m 15%. With this setup, simulations are invariant to the numerical masses of the particles, instead being sensitive to their ratio: After running a simulation, we can set M PBH to a physical value which then fixes the volume, L 3 , of the simulation: whereρ cr is the comoving critical density. 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 PBH and cosmological parameters, this is achieved by minimizing the number of PBHs in the simulation, i.e. N PBH = 1. By doing this we no longer accurately follow Poisson fluctuations in the PBH density field; however, from the Epstein mass function [35,36] describing the Poisson distribution we can deduce that PBHs rarely interact when f PBH (1 + z) × 10 −4 . We performed an explicit test of this by running a simulation with f PBH = 10 −5 and N 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 PBH : where the last equality utilises our maximum number of particles: 2 × 512 3 . We therefore consider f PBH = 10 −5 , 10 −6 and 10 −7 with N DM = 2 × 384 3 , 2 × 512 3 and 2 × 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 PBH apart for numerical resolution effects. We show the density profile as a function of f PBH at a = 10 −3 and 10 −2 in Fig. 2. As expected, we find that the profiles are independent of f PBH . We also find that the best resolution is obtained for the f 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 [29]; 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 simulations by Adamek et al. [37], who also find a smooth transition to standard NFW-like profile at large radii.
Our chief goal is determining r B,eff . To do this, we first interpolate the particles to a grid using the Cloudin-Cell method and then solve Poisson's equation for the gravitational potential: where i can indicate PBH and CDM separately and ∇ −2 ρ is evaluated in Fourier space. We then find r B,eff by plugging the potential above in Eq. (13). We show the obtained potentials at z = 99 in Fig. 3, with r B,eff defined as the intersection of the total potential (solid black) with the horizontal grey line. We remark some numerical artifacts: For the PBH potential we know the exact result, −φ PBH (r)/(GM 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 B,eff . At high redshifts the PBH is much more relevant, r 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 semi-analytical model described in the following section, which we later use to cover PBH masses whose Bondi radii are not resolved by our simulations.

C. Semi-analytical model
In the specific case of a point-like potential due to the PBH plus the power-law matter distribution around it, with density ρ(r) ∝ r −α up to a distance r h and total mass M h , Eq. (13) rewrites where p = 3−α, and M h and r h depend from {M PBH , z}. We adopt Eq. (16) for the halo mass within the turnaround radius, where the turnaround radius is (see We identify r h = r t.a. , in order to have a self-consistent normalization of the mass. Eq. (21) admits either the solution Note that Eq. (24) tends to r B,h when p → 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 B,eff ≥ r 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 3 . We thus use this model with p = 0.75 to compute the impact of PBH accretion onto the CMB.

A. Impact of accretion onto DM halos
We modify the branch ExoCLASS [13] 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 [15]. We extrapolate with constant values of the lower (higher) boundary at lower (higher) redshifts. We checked that this has sub-percent 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 PBH /M = 100, 1000. Since we find that including the halos increase 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 PBH such that each case shows the 95% C.L. exclusion. 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. [9] 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 [38,39], interfaced with our modified version of CLASS. We perform the analysis with a Metropolis-Hasting algorithm, assuming flat priors on {ω b , ω cdm , θ s , A s , n s , τ reio , f PBH } at fixed PBH mass M/M = [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 free to vary to determine the minimal PBH mass probed by cosmological data through accretion. In practice we find that M 95% min = 15M > 10M in the spherical case, we therefore did not run with M = 10M in that accretion scenario. We adopt the Planck collaboration convention and model free-streaming neutrinos as two massless species and one massive with M ν = 0.06 eV. Our data set includes Planck 2018 high-and low-TT, EE and lensing likelihood [40,41]; the isotropic BAO measurements from 6dFGS at z = 0.106 [42] and from the MGS galaxy sample of SDSS at z = 0.15 [43]; the anisotropic BAO and the growth function f σ 8 (z) measurements from the CMASS and LOWZ galaxy samples of BOSS DR12 at z = 0.38, 0.51, and 0.61 [44]. Additionally, we use the Pantheon 5 supernovae dataset [45], 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 [46] 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 [47] 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 [48]) bounds typically tighten [49,50], as we explicitly checked for the CMB ones in our previous article [9].
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 PBH < 2.9 × 10 −9 (L acc = L E ). (25) This is an optimistic benchmark for what is presumably the best limit that CMB can yield to. Our more realistic constraints at 95% C.L. are shown in Fig. 6. Strikingly, the formation of a DM halo around the PBH can increase the bound by up to ∼ 2 orders of magnitude in the covered mass range. The improvement in the constraints degrade as the fraction of PBH making up DM increases, with the bound becoming less and less affected by the presence of a halo at f PBH > 1%. It is also worth noting that for high f PBH , a non-negligible fraction of the DM can be gravitationally bound to two or more PBH, so that the isolated-PBH approximation breaks down [34], and the bounds in presence of halos should be considered as indicative. Our results also show that the bounds flatten when M > ∼ 10 4 M . 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. Note that the bounds on PBH in absence of DM halos (dark shaded regions) are themselves stronger than bounds previously derived in [9] 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.
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 "secondary" GWs generated in conjunction with their formation, falling in the frequency probed by pulsar timing arrays. The non observation of a stochastic signal in the nHz range sets tight bounds [51,52]. Galactic microlensing constraints [53][54][55], roughly excluding f PBH > ∼ O(0.1), also apply. Other "direct" bounds come from the non-observation of mergers by LIGO/Virgo [56]. 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) [57], but one may expect similar or tighter constraints from the extrapolation of the analysis of Ref. [56] to higher masses.
In the 10-100 M range, the binary coalescence rate inferred by LIGO/Virgo is estimated to yield bounds at a level between 10 −3 and 10 −2 [2,3]. Other constraints at M ∼ O(10) M roughly in the ballpark of f PBH < ∼ O(0.1) come from the non-observation of a stochastic gravitational wave (GW) background (due to the mergers of PBH binaries at high-z, in the matter dominated era) [58], quasar microlensing [59], lensing of type-Ia supernovae [60], or the orbital dynamics of halo wide binaries [61]. When approaching the ∼ 100 M scale, radio and X-ray observations of the Milky Way [62], the half-light radius of dwarf galaxies [63,64] or the stellar distribution of dwarf galaxies [65] take over as more and more stringent bounds. Basically, the CMB constraints surpass all these at masses M > ∼ 20 − 50 M , and remain the dominant constraint until at least 10 3.5 M , when they become comparable to (or slightly worse than) BBN ones [66][67][68], before being definitely surprassed by CMB spectral distortions (see [69,70] and Refs. therein) at M > ∼ 10 4.5 M . 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 > ∼ 15 M as the totality of DM; around 30 M , no more than f PBH ∼ 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 [2,3]. 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 -are believed to sit at  [52], Icarus ones from [57], LIGO ones according to [3], BBN bounds from [68], spectral CMB distortions from [69]. The arrow indicates that for masses M > ∼ 200 M , PBH can in principle grow in mass up to 10 9 M by z = 7.5 by accreting baryons at Eddington luminosity with = 0.1.
the center of almost all galaxies, and their integrated accretion disk emission almost saturates the cosmological X-ray background, see for instance the review [71]. Additionally, SMBH-including a few very massive ones with M > ∼ 10 9 M -have been observed at high redshift z > ∼ 6, a fact which seriously constraints their formation mechanism [72].
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 ≤ z < ∼ 6. Since SMBH with M > ∼ 10 9 M 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,4,5), if 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 6 This is probably unrealistic, since it excludes any sizable outflows.
converted into radiation at Eddington limit, i.e.
Note that the Eddington mass accretion rateṀ E depends on the unknown quantity , and is mathematically unbounded from above, when → 0, while it can be arbitrarily small, when → 1. However, the largest known values of that can be attained are 0.42 for a maximally rotating Kerr BH, so that the following inequalitẏ 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 where τ E was given in Eq. (4). For a typical benchmark value ≈ 0.1, this implies a 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 15 and observations at z 6. Thus, there is barely the time for a stellar BH of mass O(100) M to grow to the size of the heaviest SMBH at z > ∼ 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 Ref. [73,74]. In general, accretion in this regime may become nonstationary, see for instance the considerations in Ref. [26]. 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 [75,76], or direct collapse of BH from a gas cloud [77], see Ref. [78] for a review.
Nonetheless, the difficulty of achieving the required conditions has led several authors (see for instance Ref. [79]) 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. [70].

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. [80], see also the gray band in Fig. 5 of Ref. [81]). However, as we just reviewed, it is known that SMBH undergo significant growth with time. In fact, at z 6 the overall mass density into SMBH above 10 6 M 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, with inferred values at z = 6 of κ = 1.23 × 10 −8 Mpc −3 , α = −1.03 and m ≡ M/M * , with M * = 2.24 × 10 9 M (see Ref. [82] or equivalently Fig. 2 in Ref. [83]). This is consistent with the inferred co-moving density > 1.1 × 10 −9 Mpc −3 above 10 9 M between z = 6.44 and z = 7.44 reported in Ref. [82]. If translated in terms of the DM fraction, Eq. (29) yields about 96M Mpc −3 above 10 6 M , equivalent to a fraction of the DM abundance in SMBH above 10 6 M of f PBH 2.9 × 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 counter-arguments 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 10 5 (t/s) M . 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 × 10 4 M < ∼ M < ∼ 5 × 10 13 M is subject to tight constraints coming from CMB spectral distortions [69]. No cosmologically relevant abundance is allowed in this range unless the PBH form out of highly non-Gaussian tail fluctuations [84][85][86].
In summary, this discussion leaves out 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 6 and today.
2 Primordial SMBH seed hypothesis: If PBHs form with masses M < ∼ 10 4 M and f PBH < ∼ 10 −9 , they are consistent with all present bounds, also for initial Gaussian conditions. The presence of SMBHs with masses M > ∼ 10 9 M at redshift z > ∼ 6 requires then a mass growth by a factor of at least 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 200 M if accreting at Eddington limit and with = 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 < ∼ 15) and, at least for popIII (as opposed to direct collapse) scenarios, by a seed mass not usually exceeding 10 2 M .
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 [37,[87][88][89] 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 semi-analytical 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 ∼ 10 −3 → 10 −1 level) and the tight cosmological bounds from primordial nucleosynthesis and CMB spectral distortions at masses M > 10 3.5 M (for a summary, see Fig. 7). The CMB anisotropy bounds reach very deep down in the f PBH range, becoming as stringent as < ∼ 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 hydro-dynamical simulations, which would be particularly useful to explore the M > 10 3 M range, but also to assess if the CMB bounds exclude a PBH origin of the LIGO/Virgo merger events.
Nonetheless, 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 > ∼ 10 3 M 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 ∼ 30. Qualitatively, we can thus expect interesting implications for the dark ages, such as non-standard 21 cm and reionization epoch, which will tremendously surpass next generations CMB observations in constraining e.m. energy injection in the dark ages and therefore definitely deserve dedicated studies.