# Primordial black holes as a source of extremely high energy cosmic rays

## 1 Introduction

Small Primordial Black Holes (PBHs), with masses well below the self-gravitational collapse limit and possibly as low as the Planck Mass ($`M_{Pl}5.5\times 10^5g`$) may have formed in the primordial universe . Numerous processes, compatible with standard cosmological scenarios, can be put forward to explain their formation . In particular, if they result from initial density perturbations (with an initial mass determined by the horizon mass at this epoch), the mass spectrum can be analytically determined  following the natural hypothesis of scale-invariant Gaussian fluctuations:

$$\frac{\mathrm{d}^2n}{\mathrm{d}M\mathrm{d}V}=(\alpha 2)M^\alpha M_{evap}^{\alpha 2}\mathrm{\Omega }_{PBH}\rho _{crit}$$
(1)   
where $`M_{evap}10^{15}`$ g is the mass of a PBH evaporating nowadays, $`\mathrm{\Omega }_{PBH}`$ is the current density of PBHs in units of critical density $`\rho _{crit}`$ and $`\alpha =(1+3\gamma )/(1+\gamma )+1`$, $`\gamma =p/\rho `$ being the pressure to density ratio.

This study is dedicated to the final-stage emission of PBHs to investigate if they can be considered as candidates for extremely high energy cosmic rays (EHECR), beyond 100 EeV ($`10^{20}`$ eV). Observational data  show that the cosmic ray flux seems to be curiously unaffected by the expected Greisen-Zatsepin-Kuz’min (GZK) cutoff (due to interaction with the 2.7K cosmological background above photoproduction threshold). The integrated emission of PBHs is estimated in the following sections for a volume of universe where predicted effects of this interaction are weak (i.e. for a radius close to the attenuation length $``$ 50 Mpc).

## 2 Individual emissions

Hawking showed  that black holes can radiate particles in a process qualitatively equivalent to $`e^+e^{}`$ pairs production in a strong electric field. When the hole temperature becomes greater that the quantum chromodymanics confinment scale ($`T>\mathrm{\Lambda }_{QCD}`$), i.e. some hundreds of MeV, emitted particles are fundamental constituents rather than composite hadrons . The EHECR production by PBHs (which are particularly affected by evaporation effects because of their low mass) has to be understood in such an approach. The emission spectrum for particles of energy $`Q`$ per unit of time $`t`$ is:

$$\frac{\mathrm{d}^2N}{\mathrm{d}Q\mathrm{d}t}=\frac{\mathrm{\Gamma }_s}{h\left(exp\left(\frac{Q}{h\kappa /4\pi ^2c}\right)(1)^{2s}\right)}$$
(2)   
where contributions of angular velocity and electric potential have been neglected since the black hole discharges and finishes its rotation much faster than it evaporates  . $`\kappa `$ is the surface gravity, $`s`$ is the spin of the emitted species and $`\mathrm{\Gamma }_s`$ is the absorption probability. In the general case, $`\mathrm{\Gamma }_s`$ is a function of $`Q`$, the particle mass $`m`$, the hole mass $`M`$, and the number of degrees of freedom of the species. Its value can only be computed by numerical approximations based on expansion in spherical harmonics of the scattering matrix . In the optical limit (i.e. $`Q\mathrm{}`$) , which is totally justified for energies considered in this work,

$$\mathrm{\Gamma }_s\frac{27Q^2}{64\pi ^2(kT)^2}$$
(3)   
where $`T`$ is the ”temperature” defined by

$$kT=\frac{hc^3}{16\pi ^2GM}10^4\left(\frac{1\mathrm{g}}{M}\right)\mathrm{EeV}$$
(4)   
In such a description, the black hole behaviour mimics a black body whose temperature increases when the mass decreases until it reaches the Planck limit where this theoretical description becomes unadapted. The time-evolution of the system depends on the emitted constituents’ degrees of freedom and is therefore based on the choice of a particle physics model. It is likely that new particles, absent from the standard model, are emitted when the black hole temperature becomes extremely high, but the general behaviour remains unchanged: all the emission above 100 EeV is nearly instantaneous. The mass loss rate of a PBH is  :

$$\frac{\mathrm{d}M}{\mathrm{d}t}\frac{(7.8d_{s=1/2}+3.1d_{s=1})10^{24}\mathrm{g}^3\mathrm{s}^1}{M^2}$$
(5)   
where $`d`$ is the mass-dependant number of degrees of freedom for the emitted particles of spin $`s`$. In the standard model, $`\mathrm{d}M/\mathrm{d}t7.9\times 10^{26}/M^2`$ above the top quark production threshold. It leads to

$$\mathrm{d}t=\frac{1}{(7.8d_{s=1/2}+3.1d_{s=1})}\frac{h^3c^6}{(4\pi )^6G^3}\frac{d(kT)}{(kT)^4}$$
(6)   
or

$$\mathrm{d}t_{}1.510^{15}\frac{\mathrm{d}(kT_{})}{(kT_{})^4}$$
(7)   
where $`t_{}=t/1\mathrm{s}`$ and $`kT_{}=kT/1\mathrm{E}\mathrm{e}\mathrm{V}`$. Since it has been checked that only particles emitted at $`kT5`$ EeV will contribute (within a few percent) to the flux of cosmic rays with energies beyond 100 EeV, the characteristic production time is $`\mathrm{\Delta }t4\times 10^{18}`$ s. As a comparison, the total evaporation time for a $`10^{15}`$ g black hole is of the order of the age of the universe.

## 3 Extremely high energy emission

Taking into account formula (6) relating the temperature to the mass, the previous emission spectrum can be rewritten in its integral form  per particle species above a threshold $`E_{th}`$:

$$N=\frac{1}{(7.8d_{s=1/2}+3.1d_{s=1})}\frac{27h^2c^9}{8^6\pi ^8G^3}_{kT_i}^{kT_{Pl}}\frac{1}{(kT)^6}_{E_{th}}^{\mathrm{}}\frac{Q^2\mathrm{d}(kT)\mathrm{d}Q}{e^{Q/(kT)}(1)^{2s}}$$
(8)   
where $`T_i`$ and $`T_{Pl}`$ are the initial and Planck temperatures. It can be numerically expressed as:

$$N=1.5610^{16}_{kT_i}^{kT_{Pl}}\frac{\mathrm{d}(kT_{})}{(kT_{})^3}_{E_{th}/(kT)}^{\mathrm{}}\frac{x^2\mathrm{d}x}{e^x(1)^{2s}}$$
(9)   
where $`x=Q/(kT)`$. Fig 1 shows that the Planck cutoff is effective for energies well beyond those of interest here.

After their production, emitted quark and gluons fragment and produce a subsequent number of hadrons. Monte Carlo simulation codes tuned to reproduce experimental data obtained on colliders cannot be used because the energies considered here are several orders of magnitude greater than those available today. The multiplicity $`n_h`$ of charged hadrons produced in a jet of energy $`Q`$ is therefore estimated by means of the leading log QCD computation :

$$n_h(Q)3\times 10^2e^{2.7\sqrt{ln(Q/\mathrm{\Lambda })}}+2$$
(10)   
To get the resulting hadron spectrum, Hill  derived the following distribution:

$$\frac{dn_h}{dz}10^1e^{2.7\sqrt{ln\frac{1}{z}}}\times \frac{(1z)^2}{z\sqrt{ln\frac{1}{z}}}$$
(11)   
The number of emitted hadrons above the threshold, by a PBH of temperature $`T`$, can then be written as:

$$N_h=1.5610^{16}_{kT_i}^{kT_{Pl}}\frac{d(kT_{})}{(kT_{})^3}_{mc^2/(kT)}^{\mathrm{}}\frac{x^2dx}{e^x(1)^{2s}}_{E_{th}/(xkT)}^1\frac{\mathrm{d}n_h}{\mathrm{d}z}dz$$
(12)   
per particle species of mass $`m`$. The numerical computation has been compared to what is given by the empirical function  $`dn_h/dz=(15/16)\times z^{3/2}(1z)^2`$, leading to a multiplicity which can be easily calculated analytically. Results are in agreement within an error of 12% which is certainly not the dominant uncertainty in this evaluation. Figure 2 illustrates the general behavior of the total hadronic multiplicity $`_{E_{th}/Q}^1\frac{dn_h}{dz}𝑑z`$ above a given threshold.

The total number of emitted particles above a detector threshold $`E_{th}=100\mathrm{EeV}`$ can then be estimated by summing the direct flux (taking into account all the standard model degrees of freedom) of fundamental stable particles and the fragmentated flux resulting from the previous computation for coloured objects. The numerical result is $`F(100\mathrm{EeV})8.510^{14}`$ particles over the lifetime of a PBH.

## 4 Resulting flux above 100 EeV

The derivation of the exact resulting spectrum on the Earth is a complete study by itself, well beyond the scope of this work. It is straightforward to demonstrate that the integrated emission goes as $`E^2`$, which seems quite difficult to conciliate with the cosmic-ray experimental data if energy-dependent confinement effects are ignored. The following section therefore aims at evaluating the orders of magnitude involved.   

To derive the resulting flux reaching the earth, PBHs have been considered as classical (non baryonic) cold dark matter clustered in galactic halos. The Milky Way mass distribution is therefore assumed to follow the simple law in spherical coordinates:

$$\rho (R)=\rho _{}\frac{R_c^2+R_{}^2}{R_c^2+R^2}$$
(13)   
where $`R`$ is the distance between the considered PBHs and the Galactic Center, $`\rho _{}`$ is the local density of exploding PBHs, and $`R_c`$ is the core radius of the halo. The particle flux becomes:

$$\left(\frac{\mathrm{d}N}{\mathrm{d}t}\right)_{galactic}=\rho _{}\times F\times J(R_H,R_c,R_{})$$
(14)   
where

$$J(R_H,R_c,R_{})=\frac{1}{8\pi }_0^\pi _0^{R_H}R^2\frac{R_c^2+R_{}^2}{R_c^2+R^2}\frac{sin\varphi \mathrm{d}R\mathrm{d}\varphi }{R_{}^22RR_{}cos\varphi +R^2}$$
(15)   
$`R_H`$ being the total radius of the halo. Table 1 gives fluxes normalized to the average for extreme values of $`R_c`$ and $`R_H`$ (for $`R_{}=8`$ kpc): it shows a quite low dependance on the halo parameters.

The extragalactic contribution is computed by assuming a standard galaxy distribution $`\rho _G0.01e^{\pm 0.4}h^3\mathrm{Mpc}^1`$  (with the Hubble parameter defined as $`H_0=h100`$ km.s<sup>-1</sup>.Mpc<sup>-1</sup>). The resulting flux is:

$$\left(\frac{\mathrm{d}N}{\mathrm{d}t}\right)_{extragalactic}=K(F,R_H,R_c,R_{})\times \rho _{}\times \rho _G\times R_{GZK}$$
(16)   
where $`K(R_H,R_c,R_{})`$ is the average emission of a single galaxy (obtained by the previous method) and $`R_{GZK}50\mathrm{Mpc}`$ is the radius of a sphere ”unaffected” by the GZK cutoff. On such distances, it is not necessary to redshift energies within the expected accuracy of a few percent.

Numerical results for average values of physical parameters show that the galactic contribution is nearly three orders of magnitude larger that the extra-galactic component, even assuming the highest galaxies number density and the upper Hubble parameter limits (h$``$1). As it only depends linearly on $`R_{GZK}`$, the accurate determination of this radius is also irrelevant. The total flux above 100 EeV is then:

$$\left(\frac{\mathrm{d}n}{\mathrm{d}t}\right)_{PBH}=3.810^{23}\times \left(\frac{\rho _{}}{1\mathrm{year}^1\mathrm{pc}^3}\right)\mathrm{m}^2\mathrm{s}^1\mathrm{sr}^1$$
(17)   
Experimental data on EHECR show an integrated flux of the order of $`\left(\frac{\mathrm{d}n}{\mathrm{d}t}\right)_{exp}10^{16}\mathrm{m}^2\mathrm{s}^1\mathrm{sr}^1`$ . The required density of exploding PBHs near the earth to reproduce such a signal is then $`\rho _{}2.610^6\mathrm{year}^1\mathrm{pc}^3`$.

## 5 Discussion

Direct observational constraints on the local PBH explosion rate $`\rho _{}`$ are quite difficult to obtain. A reliable search for short bursts of ultra high-energy gamma radiations from an arbirary direction have been performed using the CYGNUS air-shower array . No strong 1 second burst was observed and the resulting upper limit, based on the exhaustive analysis of a very fine binning of the sky, is in the range $`\rho _{}0.910^6\mathrm{year}^1\mathrm{pc}^3`$. Very similar results were derived by the Tibet  and the AIROBIC collaborations . TeV gamma-rays have also been used to search for short time-scale coincidence events, thanks to the imaging atmospheric Cherenkov technique developped by the Whipple collaboration. The very high-energy gamma-ray bursts detected are compatible with the expected background, within $`\pm 1.7\sigma `$. The resulting upper limit obtained with 5 years of data , i.e. $`\rho _{}310^6\mathrm{year}^1\mathrm{pc}^3`$, is substantially better than the previous published results in the TeV range. All those limits are roughly compatible with the density required to generate the observed EHECR spectrum.

At the opposite, low-energy ($`<0.5`$ GeV) cosmic-ray antiprotons detected by a BESS 13-hours ballon flight have been used  to put a much more severe upper limit of $`\rho _{}210^2\mathrm{year}^1\mathrm{pc}^3`$, which could exclude PBHs as serious candidates for EHECR. This analysis is particularly promising since the authors have shown that the local PHB-antiproton flux can only be due to contributions from black holes that are very close to explosion, and exist within a few kpc away from the Solar system. Nevertheless, such data suffer from an important lack of statisics and from contamination effects due to interactions with the atmosphere. Future results from the AMS  spectrometer on board the International Space Station will give a much more accurate antiproton spectrum in the 0.1-1 GeV range. Those data should allow a stringent upper limit (if not a positive detection) on nearby exploding PBHs.   

An entirely different approch is to study the diffuse gamma-ray background spectrum. The emission from PBHs over the lifetime of the Universe is integrated so as to evaluate the resulting particles and radiations. This method  leads to $`\rho _{}10\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}^1\mathrm{pc}^3`$ for clustered black holes. It should, anyway, be emphasized that such a study does not directly constrain $`\rho _{}`$. The resulting ”Page-Hawking bound”  on $`\mathrm{\Omega }_{PBH}`$, derived to match the observed spectrum at 100 MeV, is converted into an upper limit on the initial number density of holes per logarithmic mass interval $`N_{PBH}`$ at $`M=M_{evap}`$ under assumptions on the Hubble parameter, on the relative matter density ($`\mathrm{\Omega }_M`$), on the equation of state of the Universe at the formation epoch, and on the gaussian distribution of initial density perturbations. This latter point is rather controvertial. The upper limit on $`\rho _{}`$ which can then be derived has to account for the large (possibly up to 8 orders of magnitude) uncertainties associated with clustering. Recent reviews on the detection of PBHs captured around massive objects  show that, when the first astrophysical objects with masses of the order of the Jeans mass were forming, black holes haloes on the sub-galactic scale could have formed around old globular clusters, dark matter clusters or population III stars. This makes the use of 100 MeV gamma-rays a quite difficult way of ruling out an important local rate of PBH explosions, though future GLAST  data should change the situation by a dramatic improvement in statistics and resolution.   

Furthermore, the first results from the AUGER Observatory  will soon give high statistics samples of EHECR. With the PBH space distribution previously assumed, the resulting EHECR flux would be from 6.0 to 2.2 times higher in the Galactic Center direction than in the opposite direction, for an integrated observation angle from 10 to 90 degrees. After five years of operation, the AUGER observatory should collect up to 300 cosmic rays above 100 EeV. Such an anisotropy would be detectable if more than approximately 50% of them come from PBHs.   

Finally, some evidences for PBH signatures available nowadays should be noted. Studies of the BATSE 1B and 1C catalogs have shown  that some gamma-ray bursts (GRBs) were consistent with a PBH evaporation origin at the quark-gluon phase transition. Characteristics of selected events are in remarkable agreement with the ”Fireball” PBH picture. The resulting (model dependent) limit is significantly lower than what is expected in the present work, and disfavours a PBH origin for EHECR. Nevertheless, new analysis of EGRET data  gives some evidences for a gamma-ray halo ”glow” due to PBH emission. Those first tentative detections are very promising for further investigations on the subject.   

From the theoretical point of view, it should also be emphasized that results given in this paper are based on the standard particle physics model. The probable increase of degrees of freedom available when the black hole temperature exceeds energies currently available on colliders would modify the estimated fluxes, making the final explosion much more violent. This could validate PBHs as a good source candidate for a fraction of the observed high energy cosmic rays.   

Acknowledgments. I would like to thank Cecile Renault for very helpful discussions.