Small-scale structure of primordial black hole dark matter and its implications for accretion

Primordial black hole (PBH) dark matter (DM) non-linear small-scale structure formation begins before the epoch of recombination due to large Poisson density fluctuations. Those small-scale effects survive until today, distinguishing physics of PBH DM structure formation from the one involving WIMP DM. We construct an analytic model for the small-scale PBH velocities which reproduces the velocity floor seen in numerical simulations, and investigate how these motions impact PBH accretion bounds at different redshifts. We find that the effect is small at the time of recombination, leaving the CMB bounds on PBH abundance unchanged. However, already at $z=20$ the PBH internal motion significantly reduces their accretion due to the additional $1/v^6$ suppression, affecting the 21 cm bounds. Today the accretion bounds arising from dwarf galaxies or smaller PBH sub-structures are all reduced by the PBH velocity floor. We also investigate the feasibility for the PBH clusters to coherently accrete gas leading to a possible enhancement proportional to the cluster's occupation number but find this effect to be insignificant for PBH around $10 M_{\odot}$ or lighter. Those results should be reconsidered if the initial PBH distribution is not Poisson, for example, in the case of large initial PBH clustering.


INTRODUCTION
Primordial black holes (PBHs) can make up the entirety or a fraction of dark matter (DM) providing a viable alternative to particle DM [1,2]. PBHs form from the gravitational collapse of large curvature fluctuations [3,4] and can thus open a window into the very early universe. Even when PBHs make up a small fraction of DM they may serve as seeds of supermassive BHs [5] or provide an origin [6][7][8][9] for the recently observed binary black hole mergers [10][11][12][13][14][15].
The abundance of PBHs is, however, constrained by several experimental observations (see e.g. [2,16,17]). Recent revisions of the femtolensing [18] and the HSC/Subaru microlensing [19] surveys have opened the mass window 10 −16 −10 −11 M for PBH DM [20,21], that may be extended to even lower masses if their radiation would be modified [22]. For higher masses PBH abundance is constrained by microlensing [23][24][25][26], survival of stars in dwarf galaxies [27,28], survival of wide binaries [29], gravitational wave observations [7][8][9]30] and the modification of CMB [31][32][33][34] or 21cm physics [35,36] due to accreting PBHs. Accretion bounds on PBHs heavier than 0.1M using CMB data were first obtained in Ref. [31]. These constraints were revised later by [33], who showed that the bounds are significantly weaker constraining PBHs heavier than 100M . These bounds depend sensitively on gas temperature, its ionization fraction and on PBH motions with respect to the gas. The motions can be broken up into several components: (i) large-scale streaming motions of the gas with respect to the DM distribution, (ii) thermal motions of the gas particles, characterized by the sound speed, (ii) small-scale PBH motions induced by the initial discreteness noise of the PBH population. The constraints from CMB [33] as well as the constraints from 21cm observations [35,36] include components (i) & (ii) in their estimates but omit the small-scale contribution (iii). The purpose of this paper is to investigate the potential impact of component (iii) on the allowed PBH mass fraction, f PBH ≡ Ω PBH /Ω DM .
After matter-radiation equality the discreteness noise of the PBH distribution drives early small-scale structure formation, leading to the formation of binaries and an early buildup of systems with multiple PBHs. Sufficiently compact PBH systems could begin accreting coherently which will lead to an enhanced accretion rate when compared to the case where the PBHs are treated as independent accretors. In particular, once their accretion radii start to overlap significantly, a system consisting of N PBHs might start to accrete as a coherent whole, leading to an enhancement by a factor of N compared to the situation with N independent accretors [37,38]. The investigation of this coherent boost factor is another task for this paper.
In this study we make use of the small-scale N -body simulations of [9] which investigated the evolution of 30 M PBHs up to redshift z 1100. To extrapolate the results towards lower redshifts we build a simple analytic model for the small-scale PBH motions, which is checked against the N -body results at z = 1100.
A monochromatic PBH mass function is assumed throughout the paper.

ACCRETION BASICS
The motion of PBHs due to the Poisson enhanced small-scale structure affects mainly the constraints on PBH abundance arising from accretion of baryons during the cosmic dark ages -the period from recombination up to the epoch of reionization -which is the period we will focus on in this paper. Due to complexities involved with low redshift structure formation -nonlinearities and baryonic feedback -we will not consider the period after reionization.
The accretion rate of a BH of mass M PBH can be cast asṀ where r a ≡ GM/v 2 a is the accretion radius, v a is a characteristic velocity and ρ is the average energy density of gas. For subluminal spherical accretion v a is the speed of sound, c s , at infinity [39]. The coefficient λ is a complicated function of redshift, the characteristic velocity and the PBH mass. For redshifts z 1100 it can decrease by an order of magnitude from 1.12 to at most 0.12 with the extremes corresponding to isothermal and adiabatic accretion, respectively [33]. When the motion of the BH is supersonic, then v a is taken to be the relative velocity between the BH and the gas, i.e. v a = v rel , and λ ≈ 0.5 [40]. An order of magnitude estimate interpolating between these two regimes can be obtained by setting The luminosity of an accreting BH, L = ηṀ , is characterized by the radiative efficiency η, which scales roughly asṀ . Thus, naively we would obtain L ∝Ṁ 2 ∝ v −6 a . However, the coefficients, e.g. λ, do also depend on v a , which can modify the velocity dependence. This dependence, obtained from the analytic model of Ref. [33], is illustrated in fig. 1. Assuming a power law dependence at a given redshift characterized by the parameter κ, i.e L ∝ v κ a , fig. 1 shows that for the naive estimate κ = −6 works relatively well when z 1000. If z 1000 the scaling is qualitatively different as the velocity dependence will asymptote to κ → 1. This transition is due to the ionization of the universe which results in the scaling λ ∝ v −3 a , characteristic of a high viscosity [41], as well as an extra linear dependence in the radiative efficiency.
We remark that in Ref. [33] the coefficients λ and η were derived under the assumption of subluminal spherical accretion. However, when z 10 4 , accretion is mostly superluminal. By relying on the analysis in Ref. [33] we thus implicitly assume that the velocity dependence of λ and η is roughly the same in the subluminal and superluminal accretion regime. The accretion bounds are sensitive to the injected power density, which we express as where f PBH is the fraction of DM in the form of PBHs, M PBH is the PBH mass and v eff is an effective populationaveraged velocity, Here f (v) is the PBH velocity distribution with respect to the baryonic medium. Unless stated otherwise, we will assume κ = −6 which agrees well with the radiatively inefficient advection dominated accretion flow (ADAF) model under the limit of low accretion speeds [42,43]. Since κ < −6 during the dark ages, the choice κ = −6 is conservative as it slightly underestimates the suppression of luminosity effect by PBH motion. The scaling of the injected power density with the PBH mass M −κ/3 PBH is also consistent with Ref. [33] as the parameters λ and η depend on the mass and velocity through combination M PBH v −3 a . We will investigate how the peculiar motion of PBHs affects accretion of gas. To obtain a relative enhancement/suppression factor, it is sufficient to know the scaling relations given in (3), but not the overall factor. Thus, our conclusions will remain valid if the PBH form an accretion disk in which case the luminosity would be enhanced but the κ ≈ −6 scaling remains intact [34].

MOTION OF THE PBH
The velocity distribution function f (v) is given by two independent components: (i) large-scale baryon-DM streaming motions described by a Maxwell-Boltzmann (MB) distribution with 3D rms velocity ∼ 30 km/s at recombination [44] with the usual (1 + z) redshift dependence, i.e., v DM−b rms min[1, (1 + z)/10 3 ] 30 km/s [45], (ii) small-scale PBH motions driven by the discreteness noise of the PBH distribution. Fig. 2 shows the PBH velocity distributions at redshift z 1100 for M PBH = 30 M PBH and for different values of f PBH obtained from the numerical simulations of Ref. [9]. The numerical results are compared to MB distributions (dashed lines) with the 1D velocity dispersions obtained from linear perturbation theory (see Appendix) They provide a decent fit to the low velocity tail, which will give the dominant contribution to luminosity during the dark ages when κ ≈ −6. The yellow dotted line shows the MB distribution for the baryon-DM stream- The low velocity tail of the PBH velocity distribution will thus dominate over baryon-DM streaming when 1 + z 50 f The extended tail of the PBH velocity distribution can be explained by the Poisson enhanced small scale structure, especially by the formation of PBH binaries from random close PBH pairs [46]. The simulation in Ref. [9] gives an overall velocity dispersion of about f 2/3 PBH 6.0 km/s for 30M PBH that remains constant after matter-radiation equality, i.e., once the early binaries have formed. Since it is possible that κ > 0 when M PBH 10 2 M , the presence of a high velocity tail may enhance the luminosity of these PBH. However, as seen in Fig. 2, the effect is present for a small fraction of PBH and, moreover, given the CMB constraint [33] and, given the velocity dispersion scales as M 1/3 PBH , the latter will not exceed 10 km/s. In comparison to 30 km/s streaming velocities, this leads to a less than 5% correction for v eff . In conclusion, before recombination, the motion of PBHs can be safely neglected for approximate estimates even when the high velocity tail is included.
From Eq. (5) one can see that the linearized continuity equation (see Appendix) dictates that the PBH velocities to grow as (1 + z) −1/2 , and thus, even though at the recombination epoch σ DM−b σ PBH , at low enough redshifts σ PBH is expected to dominate over σ DM−b , which decays as (1 + z).
In the following we assume that the velocity distribution function f (v) in Eq. (4) has a MB form, i.e. v dP/dv The colored histograms show the PBH velocity distributions at redshift z 1100 obtained numerically in Ref. [9]. Here MPBH = 30 M and from left to right fPBH = 0.01, 0.1 and 1, respectively. The dashed lines show MB distributions with dispersion (5) aiming to describe the low-velocity tails of the distributions. The yellow dotted MB curve corresponds to the large-scale baryon-DM streaming motions.
is also a good approximation, since the integral in Eq. (4) is determined by the low velocity tail, which is well approximated by the MB distribution. For σ PBH ∼ σ DM−b this approximation somewhat underestimates the effective σ. Our results for the accretion bounds are therefore conservative. A more realistic and precise treatment demands a full model for the PBH velocity distribution along with its temporal evolution.
For the sound speed c s we use the approximation where T k is the gas kinetic temperature, a the scale factor, γ = 5/3 the adiabatic index and µ 1.22 is the mean atomic weight of the neutral gas in units of the proton mass m p . T k (a) has the analytic fitting form suggested in [44]. Here T CMB = 2.725 K is the CMB temperature at z = 0 [47] and a 1 = 1/136 and a 2 = 1/181. 1 The above approximation for T k is accurate 4% when compared against the numerical results from the RECFAST Power-law index, α, for the coherent accretion boost as a function of the ratio between the mean intracluster PBH distance, d , and an effective accretion radius, ra. The points with error bars show the results from numerical simulations [38]. The green solid line represents our analytic fit used throughout this paper.
code [48]. During the dark ages σ DM−b > 2c s , that is, the streaming velocity dominates for most PBH. When c s σ, then v eff ≈ 1.03 √ c s σ.

COHERENT ACCRETION BOOST
The total accretion rate of a cluster comprising N PBHs scales asṀ tot ∝ N when all cluster members accrete independently. This approximation is valid when the average distance between the PBHs is significantly larger than their effective accretion radius. At the other extreme, that is, once the accretion radii start to overlap significantly, the cluster begins to accrete as a coherent whole. In this caseṀ tot ∝ N 2 , that is, the accretion is enhanced by a factor of N [37,38]. Between these extremal cases, the enhancement factor over the standard incoherent feeding can be approximated as N α , where 0 < α < 1. The exponent α depends on the ratio of mean PBH distance inside the clusters l to the accretion radius r a . We use the functional form for α shown in Fig. 3, which is motivated by the simulation results of Ref. [38]. 2 To obtain the mean PBH distance inside clusters we make a standard assumption that clusters are identified as objects with overdensities ∆ times over the background matter density. For the redshifts considered in this paper ∆ is well approximated by its standard Einstein-de Sitter value, ∆ = 18π 2 . Under these assump-tions, the average proper PBH distance inside clusters is The accretion radius will vary for different clusters depending on the relative motion between the PBH and the gas. To estimate the effect of motion consider the luminosity-weighted average, with weight w which for κ = −6 is half the accretion radius of a PBH at rest. The approximation is valid when c s σ and κ < −3 and can thus be used during the dark ages. The dependence on PBH velocity distribution drops out since the luminosity is dominated by the slowest PBHs.
As the accretion radius decreases faster with redshift than the average distance, the largest coherent boost is expected at smaller redshifts. The speed of sound (6) is approximated by c s ≈ (1 + z) × 15m/s around z = 10. As r a L ≈ z −1 . The condition r a L d is satisfied, when So, as the abundance of the PBH with masses over 100M are strongly constrained, coherent accretion of PBH is suppressed in viable mass ranges. An exception is possible in the central regions of halos where the density can be much above the average or for the small clusters where the average distance can drop by an order of magnitude (see Fig. 6). For binaries, the large peculiar velocity will suppress accretion when compared to the individual BHs. This effect is, however, milder for highly eccentric binaries [49]. The effect of coherent accretion on injected energy is shown in Fig. 4. At z = 10 it can be sizable already for 100M PBH. This is partly due to the large N of PBH clusters at z = 10.

The modification of the energy input by accretion is given by the factor
where the superscript 'std' represents the "standard" case with an uniform PBH distribution and without PBH peculiar motions. The average is also taken over the cluster mass distribution for which we use a discretized Press-Schechter like halo mass function (see Appendix, Eq. (16)). This is justified because the coherent accretion boost is active during the late dark ages where the cluster size can be relatively large and when the smallest PBH clusters have been evaporated due to their relatively short dynamical timescales. The results for the factor A are presented in Fig. 4 for different values of f PBH and M PBH , with and without coherent accretion boost factor N α . When the latter is omitted, then (10) is approximated by A ≈ (1 + (σ PBH /σ DM−b ) 2 ) −3/2 . The effect of the N α -term can be neglected for order of magnitude estimates when f PBH 0.1 or M PBH 10 and the dominant contribution is driven by the v eff dependence.
Consider now the effect on constraints on PBH abundance. If the accretion bound without discretenessinduced effects is f PBH,max , the corrected bound can be obtained by comparing the injected energy densities.

This gives
assuming the injected power density in (3) scales linearly with f PBH when the PBH motion is omitted. The modified boundary of the allowed region, f new PBH , is plotted in Fig. 5 for redshifts 10 and 50 for various values of M PBH . The larger the value of M PBH the more f new PBH deviates from the bound without PBH motions. It is interesting to note the non-monotonic behavior of the f new PBH -f PBH relation. This is easily understood with Eq. (3), according to which an increase of f PBH , when it is small, leads to an increase in the injected power density j. However, this increase in j will saturate once f PBH becomes sufficiently large for the PBH motions to become noticeable.
The modification on the PBH abundance bounds is shown is Fig. 5. For example, according to Fig. 5, for M PBH = 100 M and z = 10, if, neglecting PBH motions, the bound on f PBH is f PBH 0.001, then, after accounting for the corrections, the unconstrained region consists of two regions f PBH 0.001 and f PBH 0.01, i.e. a new allowed region emerges. For the slightly weaker non-corrected constraint, e.g. f PBH 0.002, however, all the values of f new PBH are allowed.

DISCUSSION AND SUMMARY
We have investigated (i) the effect of PBH motions on energy injection from the gas accretion along with (ii) possible boost due to coherent accretion inside PBH clusters. We find that the first effect has a dominant impact while the second one can be neglected for PBH masses and abundances allowed by the constraints.
At redshifts z ∼ 1000 the PBH motions are still small compared to the dominant baryon-DM streaming motions and thus the impact on existing CMB bounds, e.g. [33], is negligible. However, at redshifts z 100 the effect of PBH motions should certainly be included. Without it one ends up with unrealistically (up to several orders of magnitude) tight bounds as can be seen in Figs. 4 and 5.
We should point out that in this paper the discretenessinduced PBH motions were estimated via a linearized continuity equation, which was shown to capture the early evolution quite well. In particular, the low-velocity tail of the PBH velocity distribution, where most of the energy input from accretion is released, turned out to be well described through this simplified treatment. It is clear that nonlinear evolution adds extra small-scale motions, resulting in even stronger impact on deducible accretion bounds. Thus, our treatment here is certainly on the conservative side. A complete treatment here calls for a dedicated numerical simulation, which is beyond the scope of this paper. We assumed a monochromatic PBH mass distribution and that the primordial spatial distribution is Poisson. In case the PBH are clustered, i.e. there are non-trivial primordial correlations between PBHs, both the coherent enhancement and the correction due to PBH motion are expected to increase as the formation of structure begins earlier. In this case even the CMB bounds may be tightened, especially for heavier PBH for which the luminosity can grow when the velocity is increased. An extended mass function will non-trivially affect the evolution of PBH clusters because the heavier PBH tend to migrate towards the center of such clusters and are more likely to form hard binaries while disrupting the lighter ones. In particular, once the small scale effects become relevant, one cannot use the method of e.g. [16] to obtain constraints for extended mass functions.
The main message of this work is that discretenessinduced PBH motions must be accounted for in case one wishes to derive reliable accretion bounds in the late universe. Thus, e.g. the 21cm PBH bounds derived in [35,36], which neglect the above motions, need to be appropriately adjusted. In the wider context these discreteness driven motions provide an inescapable lower velocity floor which throughout can only grow as it evolves. Therefore, e.g. compared to the standard CDM case the halos cannot have very cool central regions, also the early stages of structure formation (z ∼ few×10), during which the shot-noise driven PBH motions have considerable impact, differ significantly. To investigate these issues in greater detail one must rely on dedicated numerical simulations. These will certainly be quite demanding, ideally requiring mass resolutions many orders of magnitude below the resolutions of typical cosmological simulations, since only then can one properly account for the physical PBH shot-noise level. A corresponding CDM simulation against which to compare the PBH run needs an even higher resolution.
Here δ c 1.686 is the critical density for spherical collapse. At the recombination epoch, z ∼ 1100, this gives M * ∼ 6M PBH f 2 PBH , thus for f PBH ≤ 0.6 we expect that the PBH will not form structures beyond binaries.
In the above calculations the halo mass was a continuous parameter. Here we use a simple method to convert to a corresponding discrete case, where the halo mass can take values M N = N M PBH , N ∈ N. In particular, we are interested in calculating a total mass fraction in halos with mass M N . Integrating Eq. (14) over halo mass M gives us a total mass density in PBHs, ρ PBH ≡ n PBH M PBH , and so the corresponding mass probability distribution function is given by Eq. (14) divided with ρ P BH . To obtain a discrete version of this probability distribution, f N ,we approximate a contribution from halos with mass M N by integrating from (N − 1)M PBH to N M PBH . This gives By construction, this probability distribution is correctly normalized, i.e.
∞ N =1 f N = 1, and that it reduces to the usual PS formalism for large N , The comparison of this analytic approximation against the N -body results of [9] are shown in Fig. 6 for f PBH = 1 and z = 1100. Although the analytic cluster mass function shows a good agreement with the numerical data, we see in Fig. 6 that the average PBH separation in small clusters can be much larger than expected from the PS formalism.
When estimating the coherent accretion boost we will use (15) as the approximation that gives This expression works well when PBH form larger clusters, i.e. when N * (z) 1, which is satisfied for viable f PBH and M PBH when z 1000.

Linear peculiar velocity field
To estimate the evolution of the low velocity tail of the velocity distribution, as characterized by (5), we start from the linear continuity equatioṅ which can be rewritten as v = −iaf (a)H(a)g(a) δ k (a = 1) k , the peculiar velocity power spectrum can be recast as Above, v is the peculiar velocity component parallel to the wavevector k, δ k is the Fourier component of the density fluctuation, a is the scale factor (normalized such that a = 1 at z = 0), H(a) ≡ȧ/a is the Hubble parameter, g(a) the linear growth factor, and f (a) is the dimensionless linear growth rate f (a) ≡ d ln g(a)/d ln a.
The dispersion of the velocity fluctuation field smoothed over the comoving scale R with filter W can be expressed as The spectrum Eq. (13) with a top-hat spatial filter, i.e. W (x) = 3(sin x − x cos x)/x 3 , gives the velocity field dispersion Taking a smoothing scale equal to the average PBH comoving distancē the corresponding 1D velocity dispersion, σ 1D v = σ v / For example, z 1100, for f PBH = 1 and M PBH = 30 M we obtain σ 1D v 0.56 km/s, which describes the low tail of the velocity distribution obtained from our simulations relatively well.