On the Primordial Black Hole Mass Function for Broad Spectra

We elaborate on the mass function of primordial black holes in the case in which the power spectrum of the curvature perturbation is broad. For the case of a broad and flat spectrum, we argue that such a mass function is peaked at the smallest primordial black mass which can be formed and possesses a tail decaying like $M^{-3/2}$, where $M$ is the mass of the primordial black hole.


Introduction
Following the detection of gravitational waves sourced by the merging of two ∼ 30M black holes [1], the idea that Primordial Black Holes (PBHs) might form all (or a fraction of) dark matter has recently obtained a lot of attention (see Refs. [2][3][4] and Ref. [5] for a recent review). Of course, to answer the question about how much of the dark matter is composed of PBHs, one has to take into account the observational constraints on the PBH abundance in a given mass range [6][7][8][9][10][11]. This requires accounting for the PBHs initial clustering [12][13][14][15][16][17][18][19][20][21] and other later phenomena as merging and accretion.
PBHs may form in the early universe if there is an enhancement of the primordial curvature perturbation ζ generated during inflation [22][23][24] at small-scales. Those small-scale perturbations are transferred to the radiation fluid after inflation through the re-heating process. When such perturbations re-enter the cosmological horizon, if they are large enough in amplitude, they can collapse and form a PBH of mass similar to the horizon mass (see, for example, Refs. [25,26] for more details on the criterion for formation).
One of the basic quantity to robustly assess if PBHs may or not form a given fraction of the dark matter is the so-called mass function ψ(M ), representing the fraction of PBHs with mass in (M, M + dM ), which is routinely defined by the relation [5]  While in the case in which the curvature perturbation power spectrum is peaked at a given wavelength the mass function is peaked at the corresponding horizon mass when that given wavelength re-enters the horizon [27][28][29][30], in the case in which the curvature perturbation power spectrum is broad, the calculation of the mass function is far less intuitive.
In the present note we elaborate on the mass function of PBHs for a broad power spectrum of the curvature perturbation. As a benchmark point, we take the broad and flat spectrum which is routinely discussed in the literature 2 parametrised as where Θ is the Heaviside step function 3 . We will call M s and M l the masses contained in the horizon when the wavelengths λ s ∼ k −1 s and λ l ∼ k −1 l re-enter the horizon, respectively.
2 For the production of PBHs from large and small-scale curvature fluctuations in single field models of inflation, a departure from slow-roll is needed [31]. A flat and broad power spectrum may be generated during a non-attractor phase when the inflaton potential is extremely flat. For those modes which exit the Hubble radius during such a phase the corresponding power spectrum is flat as a result of a duality symmetry which maps the non-attractor phase into a slow-roll phase [32,33]. 3 We will discuss the cases of a slightly red-and blue-tilted broad spectra with a non-vanishing spectral tilt np in the last section.
The formation of the PBHs is a rare event. One should imagine our universe rescaled at a given time in the past containing several Hubble volumes which are growing with time. At any instant of time there is a certain probability to form a PBH in each of these Hubble volumes. In the case of a broad power spectrum one expects that PBHs with different masses may form. This expectation is indeed correct and we argue that the computation of the corresponding mass function can be performed based on the following rather simple physical arguments: 1) the threshold for collapse is the same for perturbations entering the horizon (and collapsing) at different times; 2) the PBHs formation probability is independent of the time of collapse. In other words, PBHs of different masses have the same formation probability; 3) the so called "cloud-in-cloud" problem of a PBH being absorbed by a bigger one is irrelevant due to the low formation probability; 4) the time evolution of the mass fraction privileges the smaller PBHs formed earlier.
We will elaborate these points one by one and show that the mass function has a peak at the mass M s corresponding to the horizon mass when the shortest wavelength λ s re-enters the horizon and a tail which scales like M −3/2 for larger masses. Dynamically, what happens is that in each Hubble volume one may independently compute the probability of forming a PBH as a single event. This probability is the same for all masses and the final mass function has simply to account for the proper time evolution.
One should also recall that, in the absence of some primordial non-Gaussianity, PBHs are not initially clustered even if the spectrum is broad [19] and therefore the mass function is not initially modified by the clustering. Of course, it should be clear that we are not dealing here with the mass function obtained including much later and subsequent phenomena like merging and accretion. One roughly expects that merging would push the PBH mass distribution towards larger masses from the equality time to the current era. Accretion would increase both the abundance and the mass of PBHs. To the best of our knowledge, a thorough study is still missing and these issues will be solved most probably by running some dedicated N-body simulations (see Ref. [34] for a recent effort in this direction).
This note is organised as follows. In section 2 we deal with the threshold for collapse, in section 3 we show that the probability of formation is the same for all PBHs masses. In section 4 we deal with the cloud-in-cloud problem while we discuss the mass function in section 5. Finally, we conclude and add some comments in section 6.

The threshold for collapse is the same for all PBH masses
One of the most important ingredients to characterise the probability of collapse is the critical threshold. It is understood that, in order to formulate a criterion for collapse, the overdensities peaks must be described in real space taking into account contributions from all Fourier modes.
In linear approximation, the density contrast and the curvature perturbation are related as 4 where H(η) is the Hubble rate in conformal time and is the radiation linear transfer function accounting for the time evolution of the curvature field, and effectively smoothing out subhorizon modes. Assuming spherical symmetry, to determine the characteristic scale and amplitude of fluctuations that collapse to form a PBH, one defines the compaction function C(r) as the ratio of mass-excess within a sphere of radius r to the areal radius at r [26,42]. The criterion for collapse widely used in the recent literature states that a PBH forms if the maximum of the compaction function is above a certain threshold. It can be also shown that the amplitude of the fluctuations in terms of the excess mass within a spherical volume is equivalent to the local value of the compaction function. Indeed, at a given conformal time η, one can compute the average profile of the overdensities as with ξ(r, η) being the two-point correlator in real space if working in the comoving slicing, and δ 0 the amplitude of the perturbation. The compaction function corresponding to the mean profile can be computed in terms of the volume averaged density contrast δ(r, η) using which is time independent on super-Hubble scales. The typical size of the perturbation in real space is given by the scale at which the maximum of the compaction function is located, which is usually denoted as r m (η).
Given that PBHs form from large-amplitude fluctuations, one can apply the peak statistics [43], assuming spherical peaks (which is a good approximation since PBHs are rare events), to calculate r m . One can easily find that the size of the typical perturbation grows in time as r m (η) ∼ η due to the smoothing of modes done by the transfer function. The threshold for collapse δ c (α) is provided by numerical simulation and depends on the peaks profile through the parameter α = −(C (r m )r 2 m /4C(r m )) [26,44]. The dependence of the threshold on the parameter α is just a useful parametrisation of the shape dependence of the threshold. Now, the crucial point is that, at different times, the change in the mean profile is not enough to result in a significant change of the threshold which is δ c 0.5. This result holds for a tilted density power spectrum.
This can be understood from first principles. In fact, a flat power spectrum of the comoving curvature perturbation corresponds to a blue tilted density power spectrum P δρ/ρ b (k, η) ∼ k 4 Θ(H(η) − k) at each time η, see Fig. 1, where we also show that the rescaled compaction function for different times are basically the same. Therefore, all mean profiles are close to δρ/ρ b (r, η) ∼ sinc(H(η)r) which is the result for a narrow power spectrum for which δ c = 0.5, see also [45]. Notice that selecting different wavenumbers can be in principle done by a smoothing over a radius ∼ H −1 . In practice, we use the radiation transfer function to do this job.
We do so also motivated by the fact that the evolution of the perturbation should not depend upon the choice of the smoothing function as a matter of fact.

The formation probability is the same for all PBH masses
Let us sketch in this section the procedure to compute the probability of forming a PBH of mass M (r m (η)).
A black hole forms if δ(r m , η) ≥ δ c , and once the critical threshold is known, the abundance of PBHs can be computed using the Press-Schechter formalism [46], assuming Gaussian probability distribution for the curvature perturbation, as where ν c ≡ δ c /σ. For completeness we report the full expression of the variance σ which is with r m (η) 3/H(η) and One finds a constant ν c (P 0 ) 1/2 ∼ 0.3, see Fig. 2. This is the result of two ingredients. The first is the time independence of the threshold while the second is the effect given by the horizon growth which counteracts the action of the transfer function smoothing out subhorizon modes in Eq. (3.2). For a slightly red (blue) tilted power spectrum, the exponential dependence of the probability to ν c implies that the larger (smaller) PBHs are much more likely to form.
There exists an additional property of the PBHs formation that needs to be taken into account when computing the probability of formation of a PBH of mass M at a given epoch. Indeed, it was found that the masses of the PBHs respect a scaling relation near the critical threshold for collapse δ c as [27,47,48] where we take K = 4 as a reference value and γ = 0.36 in a radiation dominated universe [49][50][51][52][53]. The mass function of PBHs produced at a given epoch, see for example [27], possesses a low-mass tail scaling

Lighter PBHs are not absorbed by heavier ones
One of the key issues in determining the PBH mass function is that, in each of the horizon volumes, the probability that a PBH is formed and subsequently absorbed by a bigger and more massive one is totally negligible. This point is rather intuitive as the formation of a PBH is an extremely rare event, having two (or more) PBHs forming on the same site is even more improbable. To assess this point we can use the so-called excursion set method [54]. In this method the density perturbation performs a random walk as a function of the smoothing scale, and the PBH formation problem becomes the so-called first-passage time problem in the presence of a barrier. The cloud-in-cloud problem is thus taken care of by considering only those trajectories which cross the threshold for the first-time. The results of this section are not original and contained in Ref. [19], but we summarise them here for the reader's benefit.
As in the previous section, let us define the smoothed density contrast where, in this section, δ( is the window function and R identifies the particular smoothing scale. We also employed the fact that the convolution in real space becomes a product in momentum space. A particularly convenient choice for the window function is provided by a sharp filter in k-space as 5 with k f = 1/R. As we will appreciate in the following, with such a choice, the excursion set theory simplifies considerably. Using the Press-Schechter formula for the probability of collapse, one easily notices that such a probability depends only on the ratio δ c /σ(R), with σ(R) being the variance of the smoothed field. Therefore, in the excursion set scheme, it is useful to factorise the time evolution of the density contrast and absorb it in the threshold which conversely becomes a time dependent quantity denoted ω(a), which will be defined in the following. Thus, the density field in Eq. (4.1) becomes independent of time. We will discuss how this can be done in practice in the PBH scenario at the end of this section.
One then studies the "evolution" of the density contrast δ( x, R) with respect to the smoothing scale R, in a given position x, by defining where we set the coordinates such that x = 0 without loosing generality and designate δ(R) as the smoothed density contrast field at the origin. δ( k) is, by construction, a stochastic variable, and therefore also κ(R) inherits this property. Then Eq.
However, in the case of the filter being a top hat in momentum space, one achieves significant simplifications. If one defines Q(k f ) = −(1/k f )κ(k f ) and recalling the definition k f = 1/R, one finds Additionally, the variance of the time independent overdensities S can be employed as "time" variable where, in the second step, we introduced our choice of sharp filter in momentum space as in Eq. (4.2). The where the stochastic noise is characterised by the two-point function Therefore, a Langevin equation with a Gaussian white noise describes the evolution of the density field as a function of the variance. As a consequence, the evolution of δ(S) can be regarded as a Brownian random walk, with respect to S which can be considered a "time" variable, with no memory effects being present. As in Ref. [54], the particular realisation of the stochastic evolution of δ(S) will be referred to as a "trajectory" (see [19] and refereces therein for a schematic picture of this evolution).
It is a standard result that the distribution solving the Fokker-Planck equation (4.11) The hierarchy for the thresholds and the variances is then given by ω(a s ) ω(a l ) and S s S l . (4.12) However, one has that This condition implies and confirms that small PBHs are generated with the same abundance as large ones.
The "two-barriers" problem [58] permits the computation of the probability to form a virialized object or a PBHs, giving also the conditional probability that a particular trajectory, with an initial up-crossing of the barrier ω 1 at a certain value S 1 of the variance, will have a first up-crossing of ω 2 between S 2 and S 2 + dS 2 with S 1 S 2 and ω 1 > ω 2 . This is directly related to the cloud-in cloud problem. The probability that a large PBH incorporates an already formed smaller PBH with S(< S s ) at later times is given by [58] P (S, ω(a)|S s , ω(a s )) = 1 2 Erfc where we employed ω(a s )/ √ S s = O(6 ÷ 8) in the typical PBHs scenario, and the large enhancement which could be given by the exponential factor is teamed by the small coefficient ω(a)/ω(a s ) 1.
This result can be understood intuitively by realising that, since the formation of a PBH is a rare event, the probability that a small PBH is eaten by a larger one is very small, scaling like β(M (r m ))β(M s ) β 2 (M s ) for a flat broad spectrum [19]. One can wonder what happens if the mass hierarchy between M s and M (r m ) is not so pronounced. In such a case, the expression (4.14) may not be applied. However, one has to recall that the collapse starting with wavelength λ s entering the horizon lasts for a few Hubble times. This implies that the modes which enter immediately after λ s are involved in such a local collapse and therefore, in practice, the 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 -10 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 dynamics separates the masses and generates an effective hierarchy. This argument is confirmed by numerical simulations [25,26].

The PBH mass function for broad and flat spectra
Having established the absence of the cloud-in-cloud problem for broad flat spectra, and since the abundance β is the same for all masses, we are ready to compute the PBH mass distribution. We recall that where M eq 2.8 · 10 17 M is the horizon mass at the time of equality, a f indicates the scale factor at formation and from Eq. (3.4) one can use the approximate relation M 0.9M H . As we argued above, for the flat power spectrum the probability β(M ) is constant, thus one gets a M −3/2 tail for large masses, see also Refs. [59][60][61].
Notice that the same result is found by using the time evolution found in (5.2). Indeed, expanding for a mass M + dM one finds

Comments and conclusions
In this short note we have calculated the PBH mass function for the case of a flat broad power spectrum of the curvature perturbation. We have argued that the mass function should be peaked at the mass corresponding to the wavelength entering the horizon first, i.e. the smallest PBH mass, with a tail decaying as M −3/2 . This finding originates from realising that PBHs with different masses have the same probability to form, that there is no cloud-in-cloud problem and that the tail is originated by the time evolution.
Let us close with some comments. One can ask what might happen if the power spectrum is broad, but slightly red-or blue-tilted. Consider for instance a broad power spectrum of the form where positive (negative) values of the index n p give rise to a blue (red) spectrum. First of all, in both cases, the cloud-in-cloud problem is absent [19]. Second of all, we find that the threshold satisfies the relation ν c (P 0 ) 1/2 ∼ 0.3(η/η s ) np/2 . It is rather intuitive to understand that in the case of a blue spectrum, the generation of PBHs will be dominated by the smallest mass possible M s . One might think that when wavelengths λ s re-enter the horizon, PBHs are generated provided that the corresponding variance is sufficiently sizeable. This would imply however that the variance of the density contrast at the shortest scale λ s would be much larger than unity and therefore outside of the perturbative regime. We do expect therefore that broad spectra with a blue tilt predict a mass function similar to those generated by narrow spectra. For a red spectrum, the variances computed at horizon crossing for the longer wavelengths are bigger than those for the small-scale perturbations.
The generation of PBHs is therefore dominated by the largest possible mass in the power spectrum M l . One could think that the generation of the low mass PBHs could be increased at later times due to the influence of the longer modes perturbations on the small ones after they have entered the horizon [18]. However, such longer modes re-enter the horizon when the short ones have have already done so by some time and the radiation pressure has already caused the dispersion of the peaks at small scales. This effect is described by the radiation transfer function which decays as 1/(kη) 2 1 on subhorizon scales. One could also argue that PBHs might be also generated at small-scales. However, this will require the hierarchy σ s σ l and therefore the non-linear regime would have been achieved much earlier the long mode has re-entered the horizon, thus making all the computations again not reliable. It would be interesting to investigate the impact of the mass function we have derived onto the PBH constraints and the implications, e.g., for the LIGO events, see for instance Ref. [62].
Finally, let us draw the reader's attention that recently a discussion of the mass function has appeared in Refs. [41,63]. Their approach is similar and inspired by a peak theory-like approach. In particular, their condition that the compaction function has a maximum at a given smoothing scale R should correspond to the first-time passage requirement in the excursion set. The latter has also the advantage of making a well-defined prediction about the cloud-in-cloud problem. Furthermore, while Refs. [41,63] do not have an explicit expression for the mass function for a broad spectrum, we expect our results, based on simple physical arguments, to be reproducible by their formalism once the time evolution is accounted for in the linear variances they adopted.