# Bounds from Primordial Black Holes with a Near Critical Collapse Initial Mass Function

## I Introduction

Primordial black holes (PBHs) arise naturally in most cosmologies. Perhaps the least speculative mechanism for PBH formation comes from the collapse of overdense regions in the primordial density fluctuations that gave rise to structure in the universe . Thus PBHs carry information of an epoch about which we know comparatively little, and are a very useful tool for restricting theories of the very early universe, especially within the context of an inflationary scenario. In this paradigm the density spectrum is a consequence of the quantum fluctuations of the inflaton field, and can in principle be calculated given an underlying model. Thus, the non-observation of the by-products or of the effects of the energy density of these PBHs constrains the underlying microscopic theory.

The simplest bound that can be extracted from PBH formation is generated by insisting that $`\mathrm{\Omega }_{PBH}1`$. Other bounds may be derived by studying the consequences of their evaporation. Given that black holes evaporate at a rate proportional to their inverse mass squared , the phenomenological relevance of a PBH will depend upon its initial mass. Smaller mass PBHs ($`10^9<M_{bh}<10^{13}`$ g) will alter the heavy elements abundances  as well a distort the microwave background , whereas PBHs with larger masses will affect the diffuse gamma-ray background .

The net evaporation spectrum from a collection of PBHs will depend on the initial mass distribution, which in turn depends upon the probability distribution for the density fluctuations. Prior to the emergence of the COBE data, bounds from PBHs were calculated assuming a Harrison-Zel’dovich spectrum $`|\delta _k|^2k`$, with an unknown normalization. This lead to the famous “Page-Hawking” bound and all its subsequent improvements . However, assuming that the distribution is Gaussian, we can now relate the mass variance at the time of formation with the mass variance at large scales today, if we know the (model-dependent) power spectrum. None the less, bounds on $`n`$ have been derived within the class of models where the power spectrum is given by a power law $`|\delta _k|^2k^n`$ over the scales of interest . Indeed, in this way it follows that for the scale-invariant Harrison-Zel’dovich spectrum PBHs have too small a number density to be of any astrophysical significance. However, observations (see references in ) now seem to favor a tilted blue spectrum with $`n>1`$ (in CDM models) with more power at smaller scales, and thus the bounds derived from black hole evaporation can be used to constrain the tilt of the spectrum.

In this paper we revisit the aforementioned bounds in light of recent calculations which indicate that a spectrum of primordial black hole masses are produced through near critical gravitational collapse. As was pointed out by Jedamzik and Niemeyer , if PBH formation is a result of a critical phenomena, then the initial mass function will be quite different then what was expected from the classic calculation of Carr . In particular, the PBH mass formed at a given epoch is no longer necessarily proportional to the horizon mass. The resulting difference in the initial mass function leads to new bounds, which is the main thrust of this paper. In particular, we revisit the density bounds and the bounds derived from the diffuse gamma-ray observations. We will first derive bounds on $`n`$ in the class of models where the power spectrum is a power law over the scales of interest. We then relax this assumption and instead place bounds on the mass variance at the formation epoch.<sup>*</sup><sup>*</sup>*As will be discussed later, PBH formation is dominated by the earliest formation epoch if they form via critical collapse. These bounds can then be applied to a chosen model by extrapolating the variance today to the formation epoch according to the appropriate power spectrum.

Before continuing to the body of this work, we would like to point out that primordial black holes have also played a large role in attempting to explain various data. PBHs can serve as significant cosmological flux sources for all particle species via Hawking radiation. Thus it is very tempting to postulate that present-day observed particle fluxes of unknown origin are a consequence of PBH evaporation. However, to predict these fluxes, or model them realistically, we need to know the mass distribution of black holes, since their emission spectra are determined by their temperature (or inverse mass). If we assume that the black holes of interest were formed from initial density inhomogeneities generated in an inflationary scenario (which is usually assumed), then the black holes are either tremendously over-abundant or completely negligible. To get a phenomenologically interesting quantity of PBHs thus requires an extreme fine-tuning, as will be demonstrated below. Succinctly, this fine-tuning arises because the PBH number density is an extremely rapidly varying function of the spectral index $`n`$. Thus, without even analyzing the details of the spectral profile, explaining unknown fluxes via PBH evaporation is far from compelling.

## II The Initial Mass Function

Carr  first calculated the PBH spectrum resulting from a scale-invariant Harrison-Zel’dovich spectrum up to an overall normalization. Subsequently Page and Hawking calculated a bound on the normalization by calculating the expected diffuse gamma-ray spectrum from these PBHs . However, using the COBE measurements of the temperature anisotropies translated into density fluctuations (within a CDM), the overall normalization can be determined. For a scale-invariant spectrum, no significant number density of PBHs is generated. However, for a tilted blue power spectrum with more power on smaller scales, a larger number density of PBHs is expected.

Given an initially overdense region with density contrast $`\delta _i`$ and radius $`R`$ at time $`t_i`$ (using the usual comoving coordinates), analytic arguments predict  a black hole will form if

$$1/3\delta _i1.$$
(1)   
The lower bound on the density contrast comes from insisting that the size of the region at the time of collapse be greater than the Jeans length, while the upper bounds come from the consistency of the initial data with the assumption of a connected topology.

The first calculations of the PBH mass distribution assumed the relation

$$M_{bh}\gamma ^{3/2}M_h,$$
(2)   
where $`\gamma `$ determines the equation of state $`p=\gamma \rho `$ and $`M_h`$ is the horizon mass when the scale of interest crossed the horizon. Recently, numerical evidence suggests that near the threshold of black hole formation, gravitational collapse behaves as a critical phenomena with scaling and self-similarity . A scaling relation of the following form was found

$$M_{bh}(\delta )=kM_h(\delta \delta _c)^\rho ,$$
(3)   
where $`\rho `$ is a universal scaling exponent which is independent of the initial shape of the density fluctuation. It was later shown  that such scaling should be relevant for PBH formation. Indeed, in Ref.  the authors found relation (3) to hold for PBH formation when the initial conditions are adjusted to be nearly critical. They found the exponent to be $`\rho 0.37`$. They also found that for several different initial density shapes, $`\delta _c0.7`$, which is significantly larger than the analytic prediction of $`1/3`$ found by requiring that the initial overdensity be larger than the Jeans mass.

Given Eq. (3), calculating the initial PBH mass distribution becomes analytically cumbersome, since in principle one needs to sum over all epochs of PBH formation. However, we would expect that the initial mass function would be dominated by the earliest epoch of formation if we assume a Gaussian distribution with a blue power spectrum, since for larger scales the formation probability should be suppressed. This expectation was tested in Ref. , where the authors used the excursion set formalism  to calculate the initial mass function allowing for PBH formation at all epochs. The authors found that it was approximately true the that earliest epoch dominates, for the conditions of interest to us.

This simplification allows us to derive quite easily the initial mass distribution . We assume Gaussian fluctuations (the effects of non-Gaussianity will be briefly discussed in Sec. V) and define the usual smoothed density contrast

$$\delta _R(x)=d^3y\delta (x+y)W_R(y),$$
(4)   
where $`\delta (x)=(\rho (x)\rho _b)/\rho _b`$, and $`\rho _b`$ is background energy density. $`W_R`$ is the window function with support in a region of size $`R`$. The probability that a region of size $`R`$ has density contrast between $`\delta `$ and $`\delta +d\delta `$ is given by

$$P(R,\delta )d\delta =\frac{1}{\sqrt{2\pi }\sigma _R}\mathrm{exp}\left(\frac{\delta ^2}{2\sigma _R^2}\right)d\delta ,$$
(5)   
where $`\sigma _R`$ is the mass variance for a region of size $`R`$, $`\sigma _R^2=\delta _R^2(x)/R^3`$.

Then using Eq. (3), the physical number density of PBHs within the horizon per logarithmic mass interval, at the formation epoch, can be written as

$$\frac{dn_{bh}}{d\mathrm{log}M_{bh}}=V_h^1P[\delta (M_{bh})]\frac{d\delta }{d\mathrm{log}M_{bh}}=\frac{V_h^1}{\sqrt{2\pi }\sigma \rho }\left(\frac{M_{bh}}{kM_H}\right)^{1/\rho }\mathrm{exp}\left[\frac{1}{2\sigma ^2}\left[\delta _c+\left(\frac{M_{bh}}{kM_H}\right)^{1/\rho }\right]^2\right],$$
(6)   
where $`V_h`$ is the horizon volume the the epoch of PBH formation. We assume prompt reheating, and therefore take the formation epoch to be the time of reheating, which corresponds to the minimum horizon mass .

Assuming that the power law spectrum holds down to the scales of the reheat temperature, then $`\sigma _H^2R^{(n+3)}`$. We can then relate the mass variance today $`\sigma _0`$ to the mass variance at the epoch of PBH formation $`\sigma (M_H)`$, using 

$$\sigma ^2(M_H)=\sigma _0^2\left(\frac{M_{eq}}{M_0}\right)^{(1n)/3}\left(\frac{M_H}{M_{eq}}\right)^{(1n)/2},$$
(7)   
where

$$M_H=M_0\left(\frac{T_{eq}}{T_{RH}}\right)^2\left(\frac{T_0}{T_{eq}}\right)^{3/2},$$
(8)   
$`M_0`$ is the mass inside the horizon today, $`T_{RH}`$ is the reheat temperature, and $`T_{eq}`$ is the temperature at radiation-matter equality. This relation is essential to connect present-day density fluctuations to those of much earlier times. Two important assumptions underly this useful relation: First, $`n`$ is taken to be constant over all the scales of interest. Second, the universe was assumed to be radiation dominated until the temperature dropped below $`T_{eq}5`$ eV (matter-radiation equality), and then matter dominated thereafter.

From the COBE anisotropy data, the mass variance can be calculated 

$$\sigma _0=9.5\times 10^5.$$
(9)   
Using this result we can then calculate the physical number density per unit mass interval at $`T=T_{RH}`$

$$\frac{dn_{bh}}{dM_{bh}}=\frac{V_h^1}{\sqrt{2\pi }\sigma (M_H)M_{bh}\rho }y^{1/\rho }\mathrm{exp}\left[\frac{\left(\delta _c+y^{1/\rho }\right)^2}{2\sigma ^2(M_H)}\right],$$
(10)   
where

$$y=\frac{M_{bh}}{kM_H}=\frac{M_{bh}T_{RH}^2}{0.301kg_{}^{1/2}M_{Pl}^3},$$
(11)   
and $`M_{Pl}`$ is the Planck mass. The physical number density at time $`t`$ is simply Eq. (10) rescaled by a ratio of scale factors, that can be written as

$$\frac{dn_{bh}}{dM_{bh}}(t)=\left(\frac{T(t)}{T_{RH}}\right)^3\frac{dn_{bh}}{dM_{bh}}.$$
(12)   
Finally, we should note that relation (3) is only valid for $`\delta \delta _c`$. Thus we expect that we may integrate over $`\delta `$ with small errors, as long as the width of the Gaussian is sufficiently small. In particular, our results should be trustworthy provided $`\sigma 1`$, which implies $`n`$ should not exceed the maximum

$$n_{max}1+\frac{2\mathrm{log}(\sigma _0^2)}{\mathrm{log}(T_{eq}T_0/T_{RH}^2)}.$$
(13)   
We will see that $`n`$ does not exceed this maximum value for all of the bounds we consider.

## III Bounds from $`\mathrm{\Omega }1`$

Let us now calculate the total energy density in PBHs. We assume the “standard cosmology” where the universe began in an inflationary phase, reheated, was radiation dominated from the reheating period until matter-radiation equality, and then has been matter dominated. The contribution of a PBH with a given initial mass, $`M_{bh}`$, to the energy density today will depend upon its lifetime. The time-dependent PBH mass $`M(t)`$ is given by 

$$M(t)=M_{}\left[\left(\frac{M_{bh}}{M_{}}\right)^3\frac{t}{t_0}\right]^{1/3},$$
(14)   
where $`M_{}`$ is the initial mass of a PBH which would be decaying today, $`M_{}5\times 10^{14}`$ g. It is a good approximation to assume that the black hole decays instantaneously at a fixed decay time, $`t_d`$, which we use in the following.

There are two components to the PBH density bounds that we can calculate. The first is the total energy density of the PBHs that have not decayed by a given time $`t`$. The second is the total energy density of the products of PBH evaporation. The sum of these components must be less than the critical density $`\mathrm{\Omega }_{pbh,evap}+\mathrm{\Omega }_{pbh}<1`$, at any time. The evaporated products of PBHs, in particular photons, could break up elements during nucleosynthesis, disrupting the well-measured elemental abundances. This and other processes during nucleosynthesis provide additional bounds on the density of PBHs  that we do not discuss here.

The simplest bound comes from the density of PBHs that have not decayed by time $`t`$,

$$\rho _{pbh}(t)=\left(\frac{T(t)}{T_{RH}}\right)^3\rho _{tot}(t_{RH})I_{M_{}(t)}^{M_{max}}(0)$$
(15)   
where $`T(t)`$ is the temperature of the universe at time $`t`$, and $`M_{}(t)`$ is the initial PBH mass that has just completely evaporated by time $`t`$,

$$M_{}(t)M_{}\left(\frac{t}{t_0}\right)^{1/3}.$$
(16)   
$`I_{M_1}^{M_2}(\xi )`$ is a dimensionless weighted integral over the PBH mass spectrum between $`M_1`$ to $`M_2`$, normalized to the total density $`\rho _{tot}(t_{RH})`$,

$$I_{M_1}^{M_2}(\xi )=\frac{1}{\rho _{tot}(t_{RH})}_{M_1}^{M_2}𝑑M_{bh}M_{bh}\frac{dn_{bh}}{dM_{bh}}\left(\frac{M_{bh}}{M_{}}\right)^\xi ,$$
(17)   
where $`dn_{bh}/dM_{bh}`$ is given by Eq. (10). We can use the above to trivially compute the ratio of the PBH density to the critical density, $`\mathrm{\Omega }(t)`$. In particular, we need only compute the density ratio at three relevant epochs: immediately after reheating $`t=t_{RH}`$, at matter-radiation equality $`t=t_{eq}`$, and present-day $`t=t_0`$. The density ratios areNote that $`\mathrm{\Omega }(t_{RH})=I_0^{M_{max}}(0)`$ is often denoted by $`\beta (t_{RH})`$.

$`\mathrm{\Omega }(t_{RH})`$ $`=`$ $`I_0^{M_{max}}(0)`$ (18)   
$`\mathrm{\Omega }(t_{eq})`$ $`=`$ $`{\displaystyle \frac{T_{RH}}{T_{eq}}}I_{M_{eq}}^{M_{max}}(0)`$ (19)   
$`\mathrm{\Omega }(t_0)`$ $`=`$ $`{\displaystyle \frac{T_{RH}}{T_{eq}}}I_M_{}^{M_{max}}(0).`$ (20)   
where $`M_{max}`$ corresponds to $`\delta =1`$, which is of order the horizon mass at reheating. The integral should be independent of the upper limit $`M_{max}`$ if we are to trust our results.

Making the conservative approximation that all the PBH decay products are relativistic, the contribution to the density ratio of the products of PBH evaporation that has occurred up until today can be written as

$$\mathrm{\Omega }_{pbh,evap}(t_0)=\frac{T_{RH}}{T_{eq}}\left(\frac{T_0}{T_{eq}}\right)^{1/4}I_0^{M_{eq}}(3/2)+\frac{T_{RH}}{T_{eq}}I_{M_{eq}}^M_{}(2).$$
(21)   
In Fig. 1, we show the upper limit on $`n`$ as a function of $`T_{RH}`$ coming from bounding $`\mathrm{\Omega }_{pbh,evap}(t_0)+\mathrm{\Omega }_{pbh}(t_0)<1`$ (solid line). For larger values of the reheat temperature we get a more stringent bound by imposing the constraint $`\mathrm{\Omega }_{pbh}(t_{RH})1`$, simply because as $`T_{RH}`$ is increased more of the black holes will have decayed at an earlier epoch. Given that most of the energy of the decay products resides in radiation, the effect on $`\mathrm{\Omega }_{pbh}(t_0)`$ is diminished due to the redshifting. This new bound is given by the dotted line in Fig. 1.

If we assume that the PBH leaves behind a Planck mass remnant, then we have additional bounds which become important for very large reheat temperature. The best bound in this case comes from calculating $`\mathrm{\Omega }_{remnant}(t_{eq})`$ which is given by

$$\mathrm{\Omega }_{remnant}(t_{eq})=\frac{T_{RH}}{T_{eq}}\frac{M_{Pl}}{M_{}}I_0^{M_{eq}}(1).$$
(22)   
The bound in this case is shown as the dashed line in Fig. 1, and is the best bound at the largest values of the reheat temperature.

## IV Bounds from Diffuse Gamma-Rays

For a certain range of $`T_{RH}`$ we can improve our bounds on $`n`$ from diffuse gamma-ray constraints. The present day flux is determined by convoluting the initial mass function with the black hole emission spectrum

$$f(x)=\frac{1}{2\pi }\frac{\mathrm{\Gamma }_s(x)}{\mathrm{exp}(8\pi x)(1)^{2s}},$$
(23)   
where $`s`$ is the spin of the emitted particle, $`x=\omega (t)M(t)/M_{Pl}^2`$, $`\omega (t)`$ is the frequency and $`M(t)`$ is the PBH mass at the time $`t`$ of emission. $`\mathrm{\Gamma }_s(x)`$ is the absorption coefficient and may be written as $`[\omega (t)]^2\sigma _s/\pi `$. $`\sigma _s`$ is the absorption cross section and is calculated using the principle of detailed balance.

The values for $`\sigma _s`$ were calculated some time ago by Page. Let us consider how $`\sigma _s`$ behaves for massless particles. At large values of $`x`$, $`\sigma _s`$ performs small oscillations about the geometric optics limit of $`\sigma _g=27\pi M^2/M_{Pl}^4`$. As $`x`$ approaches zero, $`\sigma _s`$ goes to zero for $`s=1/2,1`$ but goes to a constant value for $`s=0`$. We will use the approximation

$$\mathrm{\Gamma }_s(x)=(56.7,20.4)x^2/\pi \text{for}s=(\frac{1}{2},1).$$
(24)   
This approximation is poor at low energies, as it is in error by $`50\%`$ at $`x=0.05`$. However, as we shall see, the contribution to the spectrum of interest is greatly peaked at $`x0.2`$. The case of strongly interacting particles is complicated by the hadronization process. There is a large contribution coming from pion decay, however, given the extreme sensitivity of the flux to the value $`n`$, the effect on the bound is negligible.

It has been recently suggested  that the self-interactions of the emitted particles will induce a photosphere, thus distorting the spectrum considerably from Eq. (23). It was suggested that two types of photospheres should form. A QCD photosphereIn the case of QCD what is meant by “photosphere” is a quark-gluon cloud. generated by parton-parton interactions as well as a QED photosphere generated by electron-positron-photon interactions. This is idea has been tested more quantitatively via a numerical solution of the Boltzmann equation . Again, while this effect may change the spectrum, especially at higher energies, it is irrelevant as far as the bound on the spectral index is concerned.

The flux measured today is given by

$$\frac{dJ}{d\omega _0}=\frac{1}{4\pi }_{t_i}^{t_0}𝑑t(1+z)𝑑M_{bh}\frac{dn_{bh}}{dM_{bh}}(t)f(x),$$
(25)   
where $`dn_{bh}/dM_{bh}`$ is evaluated at time $`t`$ using Eq. (12), $`t_0`$ is the age of the universe, $`t_i`$ is the time of last scatter, and $`f(x)`$ is the instantaneous emission spectrum given above with

$$x=\frac{\omega (t)M(t)}{M_{Pl}^2}=\frac{\omega _0(1+z)}{M_{Pl}^2}M_{}\left[\left(\frac{M_{bh}}{M_{}}\right)^3\frac{t}{t_0}\right]^{1/3}.$$
(26)   
The integral over $`t`$ is cut off at early times, since at redshifts above $`z=z_0700`$ the optical depth will be larger than unity due to either pair production off of matter or ionized matter . Those processes will degrade the energy below the window we are interested in.

This integral may be rewritten in the more illuminating form

$$\frac{dJ}{d\omega _0}=\frac{1}{4\pi }\frac{M_{Pl}^6}{(\omega _0M_{})^3}_0^{z_0}\frac{dz}{H_0(1+z)^{5/2}}_0^{\mathrm{}}𝑑xx^2\alpha ^2f(x)\frac{dn_{bh}(x,z)}{dM_{bh}},$$
(27)   
where

$$\alpha =\frac{M(t)}{M_{}}=\left\{(1+z)^{3/2}+\frac{x^3M_{Pl}^6}{[(1+z)\omega _0M_{}]^3}\right\}^{1/3}.$$
(28)   
Let us study the qualitative behavior of the above integral as a function of $`\omega `$ at fixed $`n`$ and $`T_{RH}`$. The $`x`$ integration is controlled by the Boltzmann factor in $`f(x)`$. Indeed, a little manipulation shows the the integrand is highly peaked near $`x0.2`$. Furthermore, the $`\omega `$ dependence in $`\alpha ^2`$ is almost completely canceled by the $`\omega `$ dependence in the factor $`M_{bh}^1y^{1/\rho }\alpha ^{(1/\rho 1)}`$ in $`dn_{bh}/dM_{bh}`$. Thus the $`\omega `$ dependent part of the integrand may be written as

$$\frac{dJ}{d\omega _0}\omega _0^3\mathrm{exp}\left[\frac{\left(\delta _c+aT_{RH}^{2/\rho }\alpha ^{1/\rho }\right)^2}{2\sigma ^2(M_H)}\right],$$
(29)   
where $`a^\rho =M_{}g_{}^{1/2}/(0.301kM_{Pl}^3)`$, and the only $`\omega `$ dependence in the exponential is through $`\alpha `$. If for now we assume that the dominant contribution the higher energy photons comes from recent decays ($`z0`$), and most of the support for the $`x`$ integral comes with $`x0.2`$, then $`\alpha `$ simplifies to

$$\alpha \left[1+\left(\frac{0.2M_{Pl}^2}{w_0M_{}}\right)^3\right]^{1/3}.$$
(30)   
As $`\omega _0`$ gets larger than $`0.2M_{Pl}^2/M_{}100\mathrm{MeV}`$, $`\alpha `$ becomes independent of $`\omega _0`$ and therefore the flux behaves as $`dJ/d\omega _0\omega _0^3`$. For lower energies we can make the approximation $`\alpha 0.2M_{Pl}^2/(w_0M_{})`$, and we would expect that at some point the $`\omega `$ dependence in the exponential will begin to dominate such that the flux should begin to rapidly decrease as we go to lower photon energies. The energy at which the flux turns over is determined by the competition between the two terms $`\delta _c`$ and $`aT_{RH}^{2/\rho }\alpha ^{1/\rho }`$ in the exponential, which is set by the reheat temperature. As we lower the reheat temperature the position of the kink moves to lower energies. If the reheat temperature is higher than $`T_{RH}10^9\mathrm{GeV}`$, however, the peak will stay around 100 MeV, since at these temperatures the second term in the exponential will always dominate. Indeed, we expect the position of the fall off to be near

$$\omega _{\mathrm{kink}}\mathrm{min}(100\mathrm{MeV},\frac{0.2g_{}^{1/2}T_{RH}^2}{0.301kM_{Pl}\delta _c^\rho }).$$
(31)   
Figure 2 shows the flux for fixed $`n`$ for a few different reheat temperatures. The position of the kink is well tracked by Eq. (31). Note however that the flux does not fall off exponentially at energies below the kink. This is because as we go to lower energies we pick up more of a contribution from higher redshifts.

It is interesting to contrast this behavior with the flux calculated assuming that the mass of a PBH formed at a given epoch is proportional to the horizon mass at the time of collapse. In Refs.  the authors calculated an initial mass function following the Press-Schecter formalism, summing over all epochs and assuming the relation $`M_{bh}\gamma ^{3/2}M_H`$ at each epoch. They found

$$\frac{dn_{bh}}{dM_{bh}}=\frac{n+3}{4}\sqrt{\frac{2}{\pi }}\gamma ^{7/4}\rho _iM_{H_i}^{1/2}M_{bh}^{5/2}\sigma _H^1\mathrm{exp}\left(\frac{\gamma ^2}{2\sigma _H^2}\right),$$
(32)   
where $`\rho _i`$ and $`M_{H_i}`$ are the energy density and horizon mass at $`T_{RH}`$ and

$$\sigma _H=\sigma _0\left(\frac{M_{bh}}{\gamma ^{3/2}M_0}\right)^{(1n)/4}.$$
(33)   
This result reduces to the initial mass function first computed by Page  for the Harrison-Zel’dovich spectrum with $`n=1`$ and $`dn_{bh}/dM_{bh}M_{bh}^{5/2}`$. The $`\omega _0`$ dependence of this result arises only through $`M(t)`$ given by Eq. (28). Using this initial mass distribution in Eq. (25), we expect, as in the previous case, $`dJ/d\omega _0\omega _0^3`$ for larger energies, and exponential decay into the lower energies (which will again be mollified from photons descending from higher redshifts). However, for this case the position of the kink will be fixed at around 100 MeV, independent of the reheat temperature.

Let us compare the above prediction with the recent COMPTEL and EGRET data. The EGRET collaboration found that the flux in the energy range $`30\mathrm{MeV}100\mathrm{GeV}`$ is well fit by the single power law 

$$\frac{dJ}{d\omega _0}=7.32\pm 0.34\times 10^9\left(\frac{\omega _0}{451\mathrm{MeV}}\right)^{2.10\pm 0.03}(\mathrm{cm}^2\mathrm{sec}\mathrm{sr}\mathrm{MeV})^1,$$
(34)   
while the COMPTEL data  in the range $`0.8`$$`30\mathrm{MeV}`$ can be fit  to the power law

$$\frac{dJ}{d\omega _0}=6.40\times 10^3\left(\frac{\omega _0}{1\mathrm{MeV}}\right)^{2.38}(\mathrm{cm}^2\mathrm{sec}\mathrm{sr}\mathrm{MeV})^1.$$
(35)   
Below $`0.8\mathrm{MeV}`$ there is large increase in the measured flux. Thus, the best bounds are found by comparing the measured flux to predicted flux at $`\omega _{\mathrm{kink}}`$ or at $`0.8\mathrm{MeV}`$, whichever is larger. Because of the rapid rise of the predicted spectrum relative to the measured spectrum, a change in the kink position can change the bound on $`n`$ on the order of $`0.01`$, which we consider within the accuracy of our calculation. The bounds on $`n`$ from the diffuse gamma-rays are specified by the dot-dashed line in Fig. 1. The bound terminates when all but the exponential tail of the PBHs decay prior to a redshift of $`700`$, since the optical depth at such early times exceeds unity, as discussed above.

We may compare our results to those derived by Yokoyama ,<sup>§</sup><sup>§</sup>§After this work was completed, we became aware of Ref.  that also utilized the critical collapse initial mass function to derive bounds on the PBH mass density by requiring that LSPs (in supersymmetric models) not be overproduced. where the author placed bounds on mass fraction of PBHs at $`t_{RH}`$, $`\beta (t_{RH})=\mathrm{\Omega }(t_{RH})`$, using the initial mass function, Eq. (10). He found that the bounds on $`\beta `$ did not differ significantly from the previous bounds derived using the standard initial mass functions, except for the bounds coming from diffuse gamma-rays. In the latter case, applicable for horizon masses in the range $`M_H5\times 10^{14}`$ g, he found more stringent constraints. Our bounds on $`n`$, translated into bounds on $`\beta `$, agree with his bounds coming from energy density constraints except for the case of larger reheat temperature, since we included the proper scaling of the energy density of photons emitted after PBH decay. Thus our bounds on $`\beta `$ can differ by many orders of magnitude. Our bounds coming from diffuse gamma-rays can also differ by orders of magnitude, but in this case it is for a different reason. Yokoyama determined his bound on $`\beta `$ by imposing the constraint on $`\mathrm{\Omega }_{pbh}(t_0)`$ derived in Ref. . However, when we change the initial mass function we also change the diffuse gamma-ray spectrum significantly in both shape and normalization, as discussed above. Thus, it is inappropriate to directly take the bounds from Ref.  and apply them to the case with the new initial mass function, Eq. (10). We find that our bounds on $`\beta `$ from diffuse gamma-rays are more stringent than those determined in Ref.  in the range $`M_H>5\times 10^{15}`$ g by several orders of magnitude.

## V Robustness of the Bounds

Let us consider the robustness of the bounds. We might worry that the bounds are highly sensitive to the choice of parameters given the sharpness of the initial mass function. Indeed, in the case where it is assumed that the PBH mass is given by Eq. (2), the bounds are $`n`$ are exceptionally sensitive to the exactness of this relation. This is clear from the exponential factor in Eq. (32). Given the initial mass function calculated by Jedamzik and Niemeyer, Eq. (10), we must check the sensitivity to the parameters $`\delta _c,k`$ and $`\sigma _0`$. In Ref. , the authors tested the scaling relation (3) using several different initial shapes density perturbations shapes. They found $`(\delta _c=0.70,k=11.9)`$, $`(\delta _c=0.67,k=2.85)`$, $`(\delta _c=0.71,k=2.39)`$, for Gaussian, Mexican Hat and fourth order polynomial fluctuations, respectively.

We varied the value of $`\delta _c`$ between $`0.600.80`$ and found that the bounds changed by at most $`0.01`$. The sensitivity to the variation being maximal at the smaller values of $`T_{RH}`$. Given that the initial mass function is peaked at a number smaller than $`kM_h`$, the sensitivity is increased at smaller $`T_{RH}`$ because $`\sigma `$ is an decreasing function of $`T_{RH}`$. Variations in $`k`$ are equivalent to a scaling in $`T_{RH}`$. Thus, varying $`k`$ by an order of magnitude has essentially no effect on the bound. Lastly, let us consider the sensitivity to the parameter $`\sigma _0`$. The value we used for $`\sigma _0`$ in Eq. (9) was calculated in Ref.  using the result It should be emphasized that this result assumed the spectra can be approximated as a power law over the range of $`k`$ that COBE probes. We are then making the further assumption that $`n`$ is constant down to the mass scales of relevance for PBHs.

$$\delta _0=1.91\times 10^5\frac{\mathrm{exp}[1.01(1n)]}{\sqrt{1+0.75r}},$$
(36)   
where $`r`$ is a measure of the size of the tensor perturbations. The $`1\sigma `$ observational error being $`7\%`$. The fit, Eq. (36), is good to within $`1.5\%`$ everywhere within the region $`0.7n1.3`$ and $`0r2`$. The authors of Ref.  quote a $`9\%`$ uncertainty in Eq. (36) at $`1\sigma `$, once uncertainties in the systematics and variations in the cosmological parameters are taken into account. The value in Eq. (9) was determined ignoring tensor perturbations. Given that $`\sigma _0`$ scales with $`\delta _0`$ we find that varying $`\sigma _0`$ at the $`2\sigma `$ level has no effect on our bound at the level of $`0.01`$. On the other hand, including some contribution from tensor perturbation will weaken the bound. We found that taking $`r=2`$ weakened the bound by $`0.010.02`$ throughout the range in the reheat temperature.

We can also consider the effects of a non-vanishing $`\mathrm{\Omega }_\mathrm{\Lambda }`$. Bunn et al.  extended their results to this case and found

$$\delta _0|_{\mathrm{\Omega }_\mathrm{\Lambda }}=1.91\times 10^5\frac{\mathrm{exp}[1.01(1n)]}{\sqrt{1+(0.750.13\mathrm{\Omega }_\mathrm{\Lambda }^2)r}}\mathrm{\Omega }_0^{0.800.05\mathrm{log}\mathrm{\Omega }_0}\left(1+0.18(n1)0.03r\mathrm{\Omega }_\mathrm{\Lambda }\right).$$
(37)   
If $`0r2`$, we can express $`\delta _0|_{\mathrm{\Omega }_\mathrm{\Lambda }}`$ extracted assuming a nonzero cosmological constant to a very good approximation by a scaling of $`\delta _0`$ extracted without a cosmological constant

$$\delta _0|_{\mathrm{\Omega }_\mathrm{\Lambda }}\mathrm{\Omega }_0^{0.800.05\mathrm{log}\mathrm{\Omega }_0}\delta _0,$$
(38)   
(where $`\mathrm{\Omega }_0+\mathrm{\Omega }_\mathrm{\Lambda }=1`$) and thus $`\sigma _0`$ also acquires a correction. Consequently, the bound on $`n`$ is shifted by

$$\mathrm{\Delta }nnn|_{\mathrm{\Omega }_\mathrm{\Lambda }}=\frac{2(0.80.05\mathrm{ln}\mathrm{\Omega }_0)\mathrm{ln}\mathrm{\Omega }_0}{42.9+\mathrm{ln}(T_{RH}/10^8\mathrm{GeV})}.$$
(39)   
In Fig. 3 we show the above correction as a function of $`T_{RH}`$ for several choices of $`\mathrm{\Omega }_\mathrm{\Lambda }`$. If we take $`\mathrm{\Omega }_\mathrm{\Lambda }0.7`$ as recent observational data suggests, our bounds on $`n`$ strengthen by about $`0.030.06`$ for $`T_{RH}`$ between $`10^{16}10^3`$ GeV respectively, as shown in Fig. 4.

We can also calculate bounds on the mass variance at reheating  which is essentially model-independent. If the relation Eq. (7) is violated by, for example, a power spectrum with a spectral index that depends on scale, then our previous bounds on $`n`$ cannot be applied. However, given a inflationary model one could in principle calculate the power spectrum, normalize to the COBE data at our present epoch, and then match onto the mass variance at reheating. In Fig. 5 we show the bounds on the mass variance from both the density bounds as well as the bounds from the diffuse gamma-ray observations. Notice that the diffuse gamma-ray observation bounds on $`\sigma (M_H)`$ are a significant improvement over the density bounds in the applicable range of reheating temperatures.

Finally, we must address the issue of non-Gaussianity. It has been pointed out  that skewness could very well be important for PBH formation given that its effects are amplified in the tail of the distribution $`P[\delta ]`$, which contributes to PBH formation. In general, the amount of non-Gaussianity expected is highly model dependent. Bullock and Primack investigated several inflationary models to study the amount of non-Gaussianity one would expect at larger values of $`\delta `$. They calculated $`P[\delta ]`$ for three toy models, and found in one case no deviation from Gaussianity and in the other two found a significant suppression in the probability of of large perturbations. However, as was pointed out in Ref. , while these effects can drastically effect the PBH mass fraction $`\beta `$, we expect that, even in the most extreme case considered in Ref. , the effect on the bound on $`n`$ is only at the level of $`0.05`$. For hybrid inflation, where the approximation that $`n`$ is constant actually holds, the perturbations are in fact Gaussian due to the linear dynamics of the inflaton field . Therefore, these bounds should be applied to specific models, with the roughness of the bound determined by the deviations away from Gaussianity.

## VI Conclusions

We have calculated the density of primordial black holes using the the near critical collapse mass function that results in a spectrum of PBH masses for a given horizon mass. The normalization of the PBH mass spectrum was determined using the COBE anisotropy data that allowed us to set bounds on the spectral index $`n`$ as a function on the reheat temperature. We find that restricting the density of PBHs to be less than the critical density corresponds to the restriction that the spectral index $`n`$ be less than about $`1.45`$ to $`1.2`$, throughout the range of reheating temperatures resulting after inflation, $`10^3`$ to $`10^{16}`$ GeV respectively. (The precise limits are shown in Fig. 1.) For a smaller range of reheating temperatures, between about $`10^7`$ to $`10^{10}`$ GeV, significant PBH evaporation occurs when the optical depth of the universe is less than one. Hence, we found a slightly stronger bound on the spectral index by restricting the cosmological PBH evaporation into photons to be less than the present-day observed diffuse gamma-ray flux. Due to the extreme sensitivity of the PBH mass density to the spectral index, effects such as the indirect photon flux from PBH evaporation into quarks and gluons which fragment into pions or the formation of a QCD photosphere are completely negligible when calculating the bound on $`n`$. We should also remark that slightly stronger bounds on $`n`$ for larger reheating temperatures $`10^{10}`$ GeV are expected from PBHs that decay during the epoch of nucleosynthesis.

If the universe is vacuum-energy dominated, there are corrections to our bounds on $`n`$ that can be substantial. We calculated these corrections for a range of $`\mathrm{\Omega }_\mathrm{\Lambda }`$ and applied them to our bounds on $`n`$ for the case of $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. The improvement is apparent by contrasting Fig. 1 with Fig. 4. Finally, we calculated bounds on the mass variance at reheating. These bounds in principle could be used to constrain any given inflationary model, once the power spectrum is calculated.

###### Acknowledgements.

This work was supported in part by the Department of Energy under grant number DOE-ER-40682-143. We thank Rich Holman and Jane MacGibbon for useful discussions. We also thank Andrew Liddle useful discussions and comments on the manuscript.