Search for Gamma-Ray Emission from Local Primordial Black Holes with the Fermi Large Area Telescope

The formation of primordial black holes (PBHs) is a prediction of some models of the early universe (Zel'dovich & Novikov 1966; Hawking 1971). In this paper, we search for evidence of gamma-rays produced by the Hawking radiation of low-mass, high-temperature PBHs in the Fermi Large Area Telescope (LAT) data.

In classical general relativity, black holes (BHs) have zero temperature and do not emit particles. The concept of a black hole temperature was formally introduced as a way to resolve the paradox of information loss for matter falling into the black hole (Bekenstein 1973, 1974). Hawking (1974, 1975) argued that black holes should emit particles due to pair creation near the horizon with gray body emission spectra with the temperature of

$$\begin{eqnarray}&&{T}_{\mathrm{BH}}=\displaystyle \frac{{\hslash }{c}^{3}}{8\pi {GMk}}\approx {10}^{-7}{\left(\displaystyle \frac{M}{{M}_{\odot }}\right)}^{-1}K,\end{eqnarray} \tag{ 1 }$$

where M is the black hole mass and M_⊙ is the solar mass. Astrophysical black holes, such as those created in a collapse of a massive star, have masses larger than M_⊙. The corresponding temperature is less than 10⁻⁷ K, which is less than the fluctuations in the cosmic microwave background (CMB) radiation; this makes the signal from these BHs extremely difficult to detect. However, PBHs with masses as small as the Planck mass (∼10⁻⁵ g) may have been created in the early universe (Zel'dovich & Novikov 1966; Hawking 1971).

The mass of a PBH created at time t after the Big Bang is proportional to the mass within the particle horizon at time t (for a review, see Halzen et al. 1991; Carr 2005; Khlopov 2010),

$$\begin{eqnarray}&&M(t)\approx {10}^{15}\left(\displaystyle \frac{t}{{10}^{-23}\,{\rm{s}}}\right)\,{\rm{g}}.\end{eqnarray} \tag{ 2 }$$

Limits on the density of PBHs at mass M constrain the magnitude of density fluctuations (or the equation of state) at time t. These constraints are complementary to the constraints obtained from observations of the CMB since they are sensitive to fluctuations on spatial scales much smaller than those of the scales observed in the CMB (e.g., Linde et al. 2013).

Since Hawking radiation is emitted in all available particle species, the total emitted power and therefore the lifetime of a PBH depends on the number of available particle states. For the Standard Model of particles, the lifetime of a PBH can be approximated as (Carr et al. 2010)

$$\begin{eqnarray}&&\tau \approx 400{\left(\displaystyle \frac{M}{{10}^{10}{\rm{g}}}\right)}^{3}{\rm{s}}.\end{eqnarray} \tag{ 3 }$$

The PBHs that were formed in the early universe with a mass around M_* = 5 × 10¹⁴ g have a lifetime close to the lifetime of the universe and reach late stages of evaporation at the present time. The PBHs with the remaining lifetime of months to years will be the subject of the search in this paper. The temperature corresponding to M_* is T_BH ≈ 20 MeV (MacGibbon 1991; Carr et al. 2010).⁶² A limit on the average cosmological density of PBHs can be obtained by integrating the flux from PBHs over the lifetime of the universe. This flux has to be smaller than the observed extragalactic gamma-ray background (Carr & MacGibbon 1998; Carr et al. 2010). This method is most sensitive to PBHs with initial masses of M ≈ M_*. Since PBHs are expected to be concentrated in galaxies similarly to dark matter, stronger constraints for initial masses slightly larger than M_* can be obtained by searching for diffuse emission from PBHs in the halo of the Milky Way Galaxy (Lehoucq et al. 2009). A comparison of extragalactic and Galactic gamma-ray background constraints with the constraints from the search for individual PBHs will be presented in Section 5.

Constraints on the local PBH evaporation rate have been obtained by looking for bursts of high-energy gamma-ray emission with a duration of a fraction of a second to several seconds. In particular, PBHs with a mass of M ∼ 10⁹ g have a temperature of T_BH ≈ (10¹⁰ g/M_BH)TeV ∼ 10 TeV (Carr et al. 2010) and a lifetime of τ ∼ 0.4 s. A search for such PBHs in the context of the Standard Model has been carried out by searching for high-energy gamma-rays (≈100 GeV < E < ≈ 50 TeV) with Cherenkov telescopes (e.g., Archambault & VERITAS 2017; Glicenstein et al. 2013; Abdo et al. 2015). In the case of the Hagedorn model (Hagedorn 1968), which was proposed before the Standard Model of particle physics was confirmed, the number of states available for the PBH evaporation increases exponentially when the temperature reaches the Hagedorn transition energy around 160 MeV. This leads to a microsecond burst of gamma-rays around a few hundred MeV (Page & Hawking 1976). A search for such bursts was carried out using EGRET data (Fichtel et al. 1994). One can also expect short radio bursts in this model (Blandford 1977; Rees 1977), which can be used to constrain the PBH evaporation rate (e.g., Cutchin et al. 2015).

Fermi-LAT is a pair conversion telescope that is sensitive to gamma-rays from ≈20 MeV to more than 300 GeV. Overviews of the Fermi-LAT design and performance can be found in Atwood et al. (2009) and Ackermann et al. (2012). PBHs with masses of M ≲ 10¹⁵ g have temperatures of T ≳ 10 MeV and emit gamma-rays (Page & Hawking 1976) that can be detected with the Fermi-LAT.

In this paper, we search for PBHs in the context of the Standard Model of particles using the Fermi-LAT gamma-ray data. In Section 2, we estimate that for the differential point source (PS) sensitivity in the 4 years of observation the Fermi-LAT is most sensitive to PBHs with a temperature of T_BH ∼ 16 GeV (corresponding to a remaining lifetime of ∼4 years). The corresponding distance to which such a PBH can be detected is ≲0.02 pc. Assuming a random relative velocity of PBHs around Earth similar to that of particle dark matter velocity dispersion, the displacement of a PBH over 4 years within 0.02 pc is, on average, greater than 1°. For comparison, the PS localization radius for a near-threshold source is approximately 01 (Ackermann et al. 2012). Thus a smoking-gun signature of a PBH in Fermi-LAT data is a moving source with a hard spectrum of gamma-rays. Moving γ-ray sources in the context of dark matter microhalos have been discussed (Koushiappas 2006), and an algorithm for statistical detection of sub-threshold moving sources using the two-point correlation function was developed (Geringer-Sameth & Koushiappas 2012). However, a search for individual bright moving sources is a novel feature of our analysis that has not been reported in the literature related to PBH searches.

In Section 3, we search for PBH candidates among the sources of the third Fermi-LAT catalog (3FGL; Acero et al. 2015). The majority of sources can be excluded as PBH candidates based on their association with known sources or on their spectra, but several sources cannot be excluded based on these criteria alone. For these sources, we analyze the Fermi-LAT data near the position of the PBH candidates and find one candidate whose proper motion is inconsistent with zero. We discuss this source in more detail in Section 3.2. In Section 4, we use Monte Carlo (MC) simulations to derive the efficiency of our PBH selection criteria, which we then use to derive an upper limit on the PBH evaporation rate. We find that PBHs can be detected up to a distance of ∼0.03 pc at the expense of a reduction of volume in phase space: the PBHs should be moving preferentially along the line of sight to be detected as PSs. The detectability distance of ≈0.03 pc derived with MC simulations and the full PS analysis pipeline is generally consistent with the simple estimate of 0.02 pc, which we obtain for nonmoving PBHs using the known differential sensitivity of Fermi-LAT to PSs. In Section 5, we present the conclusions. In the Appendix, we discuss some uncertainties on the spectrum of gamma-rays emitted by PBHs.

We start our analysis by estimating the sensitivity domain, i.e., the relevant energies and timescales of the Fermi-LAT to individual PBHs. One of the main questions is whether the Fermi-LAT is more sensitive to a population of PBHs with a low temperature (e.g., ≲10 GeV), hence, with a smaller intensity of emission but a large spatial density, or to a population with a high temperature (e.g., ≳10 GeV) but low spatial density.

In this section, we derive the range of masses, temperatures, and distance to Earth where PBHs are detectable by Fermi-LAT. In Figure 1, we compare the spectra of PBHs with the differential PS sensitivity for four years of Pass 7 reprocessed data (Ackermann et al. 2012; Bregeon et al. 2013). The same data sample was used in the derivation of the 3FGL catalog (Acero et al. 2015), which we use in the next section to search for PBH candidates. To determine the four-year equivalent flux from a PBH, we integrate the PBH spectrum either over the lifetime of the PBH or over four years, whichever is smaller, and divide by four years. The PBH spectra are derived taking into account both primary and secondary (mostly from hadronic showers of quarks and gluons) production of gamma-rays by the BHs (MacGibbon & Carr 1991). Although all Standard Model particles can be emitted by PBHs in principle, the γ-ray spectrum is dominated by the QCD degrees of freedom (MacGibbon et al. 2015). The normalization for each of the curves is chosen such that the PBH flux is equal to the differential PS sensitivity in one of the energy bins.

Zoom In Zoom Out Reset image size

Since the intrinsic luminosity of PBHs for a given temperature is fixed, we can use the condition that the flux should be larger than the PS sensitivity to estimate the maximal distance at which the PBHs can be detected as a function of the initial temperature (shown by the solid line in Figure 2 on the right). Typical values near the maximum are ≲0.02 pc for initial temperatures between 10 and 40 GeV. The corresponding remaining lifetime is between  20 years and  4 months. Over the age of the universe, we expect PBHs to virialize in the same manner as dark matter. If we take into account the orbital velocity of the solar system around the GC v_rot ∼ 250 km s⁻¹ (Wilkinson & Evans 1999; Xue et al. 2008; Sofue et al. 2009; Gnedin et al. 2010; McMillan & Binney 2010; McMillan 2011) and add a velocity dispersion similar to the expected dark matter velocity dispersion near the Sun v_disp ∼ 270 km s⁻¹ (e.g., Kuhlen et al. 2010), then the displacement of a PBH during 3 years at R = 0.02 pc from the Earth could be as large as α = v_perp t/R ∝ 3°. Here we estimated the average velocity perpendicular to the line of sight in a random direction on the sky as ${v}_{\mathrm{perp}}=\sqrt{2/3({v}_{\mathrm{rot}}^{2}+{v}_{\mathrm{disp}}^{2})}\,\sim 300\,\mathrm{km}\,{{\rm{s}}}^{-1}$. Because this displacement is ∼15 times larger than the Fermi-LAT PS localization radius for a threshold source, such a black hole will appear as a linear streak of gamma-ray emission instead of as a PS.

Zoom In Zoom Out Reset image size

Since the large majority of Fermi-LAT sources are point-like, linear spatial extension is a powerful criterium for identifying PBHs. One of the difficulties in observing a moving PS is that the flux becomes smeared out into an extended track following the trajectory of the source, but PS catalogs are optimized for sensitivity to PSs as opposed to extended sources. A simple estimate of the Fermi-LAT sensitivity to a moving source can be obtained by integrating the flux during the time when the moving source appears as a point-like source, e.g., when the displacement is less than the PS localization radius. The corresponding time is shown in Figure 2 (left).⁶³ Using the reduced integration time, we also derive the corresponding detectability distance (red dashed–dotted line in Figure 2 on the right). For large temperatures (small remaining lifetime) the detectability distance for a moving source is the same as that for a stationary one because the expected displacement is smaller than the PS localization radius.

Given the characteristic detection radius for PBHs, including the proper motion R ∼ 0.01 pc for 4 years of observations (Figure 2, right panel), we estimate that the Fermi-LAT is sensitive to a PBH evaporation rate of $\dot{\rho }=1/{Vt}\propto 6\,\times {10}^{4}\,{\mathrm{pc}}^{-3}\,{\mathrm{yr}}^{-1}$, where V = 4π R³/3 is the detectability volume and t = 4 years is the observation time. The main result of this section is that the Fermi-LAT sensitivity is potentially competitive with the sensitivity to individual PBHs of Cherenkov observatories (Archambault & VERITAS 2017; Glicenstein et al. 2013; Abdo et al. 2015; Ukwatta et al. 2016), e.g., $\dot{\rho }=1.4\times {10}^{4}\,{\mathrm{pc}}^{-3}\,{\mathrm{yr}}^{-1}$ derived by the H.E.S.S. collaboration (Glicenstein et al. 2013), but one has to take the proper motion of PBHs into account. In the following sections, we use the 3FGL catalog, which employs the full PS sensitivity of the Fermi-LAT rather than the differential sensitivity. We also treat the proper motion more rigorously. As a result, the derived Fermi-LAT sensitivity is a factor of a few better than the simple estimate presented in this section.

The Fermi-LAT surveys the entire sky approximately every three hours, and has relatively uniform exposure over long timescales (at 1 GeV over 4 years the exposure varies by approximately 40% over the entire sky). Combined with a large effective area (approximately 1 m² between 1 GeV and 1 TeV), this makes it an ideal instrument for detecting a large number of gamma-ray PSs. The most complete catalog of PSs is currently the third Fermi PS catalog (3FGL, Acero et al. 2015). It contains sources that are significantly detected above 100 MeV and spans the first 4 years of the Fermi mission. We used the 3FGL catalog to search for PBH candidates and to constrain the local PBH evaporation rate.

To find PBH candidates in the 3FGL catalog, we first excluded from further consideration PSs associated with known astrophysical sources (such as blazars). Of 3033 sources in the 3FGL catalog, 1010 are unassociated sources. We also excluded sources that are within 10° of the Galactic plane, which removed a further 468 sources. Our analysis was restricted to high-latitude sources because (a) detectable PBHs are expected be distributed isotropically given the detectability distances estimated in Section 2, while astrophysical sources are concentrated along the Galactic plane, and (b) association of extragalactic sources such as blazars is easier at high latitudes (see, e.g., Ackermann et al. 2015).

To test the remaining unassociated 3FGL sources as PBH candidates, we fit their spectra with the time-integrated gamma-ray spectra emitted by a PBH. For this analysis, we used the fluxes of the candidate sources as reported in the 3FGL, which are provided in five energy bands (0.1–0.3, 0.3–1, 1–3, 3–10, and 10–100 GeV). The time-integrated PBH spectrum depends on two parameters: initial mass (or temperature) and the distance to the Earth (equivalent to an overall normalization). We varied these two parameters to obtain the best fit to the PS spectrum in the five energy bands. The quality of the spectral fit to the reported spectrum was determined by calculating the value of the χ² over the five energy bands:

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=1}^{5}\displaystyle \frac{{({{\rm{\Phi }}}_{i}-{{\rm{\Psi }}}_{i})}^{2}}{{\sigma }_{i}^{2}},\end{eqnarray} \tag{ 4 }$$

where Φ_i is the PBH spectrum integrated over the width of bin i and Ψ_i is the flux in bin i from the 3FGL. Here σ_i represents the uncertainty on the flux in bin i. In the 3FGL, flux uncertainty is represented by a 68% confidence interval; the value of the uncertainty can then be written as the difference between the best-fit flux and either the upper or lower bound. We set σ_i to be the larger of the two in order to be conservative. We require the value of the best-fit χ² to be below the critical value of 11.3, which corresponds to 99% exclusion for five degrees of freedom and two parameters. In other words, sources with χ² values greater than 11.3 have only a 1% chance of being spectrally consistent with a PBH. After the spectral consistency was computed, 318 sources out of the 542 unassociated candidates remained as PBH candidates.

The candidate sources next underwent a check for proper motion. We use the following algorithm to determine the magnitude and significance of proper motion:

• 1.  
  All source-class photons above 1 GeV within 5° of the source's reported 3FGL location were collected. The time range (2008 August–2012 July) and data reconstruction (P7REP_SOURCE_V15) were consistent with that of the data used to construct the 3FGL. Since the angular resolution of the Fermi-LAT decreases quickly below 1 GeV, including photons below 1 GeV did not have a significant impact on the final results.
• 2.  
  Data covering a longer time range (2008 August–2017 July) and a more recent event reconstruction (P8R2_SOURCE_V6) were held in reserve for validation, and were used to test the PBH hypothesis for any sources that passed the proper motion cut.
• 3.  
  The expected number of photons N from the source of interest was calculated by multiplying the flux in each energy bin by the Fermi-LAT exposure at the bin's midpoint energy, and summing over the three relevant bins (1–3, 3–10, 10–100 GeV).
• 4.  
  In order to estimate the velocity of a PS, we compare the maxima of the likelihood function ${ \mathcal L }({{\boldsymbol{x}}}_{i},{t}_{i},{{\boldsymbol{x}}}_{0},{{\boldsymbol{v}}}_{0})$ in two cases: fixed ${{\boldsymbol{v}}}_{0}=0$ and free ${{\boldsymbol{v}}}_{0}$. We approximate the point-spread function of Fermi-LAT by a Gaussian for simplicity. The likelihood function is given by multiplying over N photons around the initial position of the source:

  $$\begin{eqnarray}&&{ \mathcal L }=\displaystyle \prod _{i=1}^{N}{w}_{i}\times {\rm{\exp }}\left\{\displaystyle \frac{-{({{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{0}-{{\boldsymbol{v}}}_{0}{t}_{i})}^{2}}{{\sigma }_{i}^{2}}\right\},\end{eqnarray} \tag{ 5 }$$

  where ${{\boldsymbol{x}}}_{i}$ is the coordinate of the photon, ${{\boldsymbol{v}}}_{0}$ is the proper motion of the source, t_i is the photon arrival time, ${{\boldsymbol{x}}}_{0}$ is the source location at the beginning of the observation time, and σ_i is the 68% angular containment radius for a photon at energy E_i, which is ∼07 at 1 GeV (Bregeon et al. 2013). Here w_i is the weight assigned to each photon, the calculation of which is defined in Section 3.1. In practice, we use the natural logarithm of the likelihood:

  $$\begin{eqnarray}&&\mathrm{log}{ \mathcal L }=\displaystyle \sum _{i=1}^{N}\mathrm{log}{w}_{i}-\displaystyle \frac{{({{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{0}-{{\boldsymbol{v}}}_{0}{t}_{i})}^{2}}{{\sigma }_{i}^{2}}.\end{eqnarray} \tag{ 6 }$$

  The main difficulty is separating the N photons attributed to the source from background photons. Our algorithm chooses a four-dimensional grid of points around an initial value of ${{\boldsymbol{x}}}_{0}$ and ${{\boldsymbol{v}}}_{0}=0$, and for each grid point finds the N photons inside the 5° ROI that have the highest contribution to $\mathrm{log}{ \mathcal L }$, i.e., the photons that most likely belong to the source given a particular position and velocity. Therefore, the weights w_i do not appear as a prefactor in Equations (5) and (6) because the N best-fit photons change given different assumptions of ${{\boldsymbol{x}}}_{0}$ and ${{\boldsymbol{v}}}_{0}$. The best-fit ${{\boldsymbol{x}}}_{0}$ and ${{\boldsymbol{v}}}_{0}$ are found by maximizing ${\rm{\Delta }}\mathrm{log}{ \mathcal L }=\mathrm{log}{ \mathcal L }-\mathrm{log}{ \mathcal L }({{\boldsymbol{v}}}_{0}=0)$ on the grid. With the additional degrees of freedom from allowing ${{\boldsymbol{v}}}_{0}$ to float, the value of ${\rm{\Delta }}\mathrm{log}{ \mathcal L }$ is always non-negative. We find in MC simulations (described in Section 4) that this algorithm tends to underestimate the input velocity by ≈25%; the best-fit velocity should therefore be considered a lower bound on the true velocity and sufficient for our purpose of separation of moving and stationary sources. The underestimation occurs because source photons that are far away from the average source position have lower weights and so are less often included in the likelihood calculation. In their stead are background photons, whose distribution in time is random, and therefore cause the algorithm to favor a slower overall velocity.
• 5.  
  The significance of ${\rm{\Delta }}\mathrm{log}{ \mathcal L }$ for each source was found by assigning random times t_i drawn from a flat distribution to each photon but fixing the positions ${{\boldsymbol{x}}}_{i}$ of all the photons, and reoptimizing ${\rm{\Delta }}\mathrm{log}{ \mathcal L }$. This process was repeated 50 times for each source, and the original value of ${\rm{\Delta }}\mathrm{log}{ \mathcal L }$ was compared with the distribution of ${\rm{\Delta }}\mathrm{log}{ \mathcal L }$ for the data sets scrambled in time to find a local significance σ:

  $$\begin{eqnarray}&&\sigma =\displaystyle \frac{{\rm{\Delta }}\mathrm{log}{{ \mathcal L }}_{0}-\overline{{\rm{\Delta }}\mathrm{log}{{ \mathcal L }}_{s}}}{{\rm{std}}({\rm{\Delta }}\mathrm{log}{{ \mathcal L }}_{s})},\end{eqnarray} \tag{ 7 }$$

  where ${\rm{\Delta }}\mathrm{log}{{ \mathcal L }}_{0}$ is the original value of the improvement in likelihood, $\overline{{\rm{\Delta }}\mathrm{log}{{ \mathcal L }}_{s}}$ is the mean of the scrambled likelihood improvements, and ${\rm{std}}({\rm{\Delta }}\mathrm{log}{{ \mathcal L }}_{s})$ is the standard deviation of the scrambled likelihood improvements.
• 6.  
  A cut on the local significance for each source was made at 3.6σ, which corresponds to a global significance of 2σ for 318 sources.

A single source (3FGL J2310.1–0557, see Section 3.2 for a discussion of this source) exceeded this cut on local significance, and the standard deviation of the local significances of the entire set of candidates was 1.03, which is consistent with statistical fluctuations. After examining the data held in reserve for J2310.1–0557 (described in Section 3.2), we concluded that no likely PBH candidates exist in the 3FGL catalog.

3.1. Calculating Photon Weights

The photon weights w_i in Equations (5) and (6) are defined as the probability that a given photon originated from the candidate PS, and are calculated by performing a standard likelihood optimization with the Fermi Science Tool gtlike.⁶⁴ The model used includes all 3FGL PS within 5° of the candidate source, as well as the standard Pass 7 models for Galactic and isotropic diffuse emission. The candidate source is modeled as an extended source with a radial Gaussian profile with σ = 025 instead of a PS, in order to account for the possibility of proper motion. The data were binned into three logarithmically spaced energy bands between 1 and 100 GeV and in 01 × 01 spatial pixels. After the model was optimized by gtlike, weights were assigned to each photon (described its coordinates x, y, and energy E) by calculating the fraction of the flux belonging to the candidate source in each pixel:

$$\begin{eqnarray}&&{w}_{x,y,E}=\displaystyle \frac{{\rm{\Phi }}^{\prime} (x,y,E)}{{\displaystyle \sum }_{i}{{\rm{\Phi }}}_{i}(x,y,E)}\end{eqnarray} \tag{ 8 }$$

where Φ'(x, y, E) is the predicted flux from the candidate source in the pixel and Φ_i(x, y, E) are the fluxes from all the sources in the model. In addition, the 3FGL PS were masked by assigning a weight of zero to all the photons that fell in a pixel more than 1° from the candidate source position where the summed contribution of the non-candidate PS fluxes exceeded 10% of the total flux in that pixel. This meant that all the photons in the calculation had a high probability of originating either from the candidate source or the diffuse background.

The weighting has little impact on the reconstruction of proper motion because individual photon weights do not change as the likelihood maximization from step 4 optimizes ${{\boldsymbol{x}}}_{0}$ and ${{\boldsymbol{v}}}_{0}$. However, weighting the photons in this way prevents the algorithm from interpreting photons from nearby sources as originating from the candidate source. Without weighting, we found that flaring nearby sources could mimic a moving source and therefore lead to false positive results.

3.2. J2310.1–0557

We used MC simulations to derive the efficiency for detecting PBHs, and used the efficiency to place upper limits on the local PBH evaporation rate. We generated a sample of PBHs within 0.08 pc of the Earth with uniform spatial density and random velocities with an average speed of 250 km s⁻¹, which is close to an upper bound on orbital velocity of the Sun around the Galactic center (Wilkinson & Evans 1999; Xue et al. 2008; Sofue et al. 2009; Gnedin et al. 2010; McMillan & Binney 2010; McMillan 2011), and three-dimensional velocity dispersion equal to the local velocity dispersion of dark matter, 270 km s⁻¹ (Kuhlen et al. 2010). At the end of this section, we also derive the limits for different assumptions about the PBH distribution, such as the relative velocity and velocity dispersion, to estimate the corresponding systematic uncertainty.

We assume that the PBH population has a constant rate of PBH evaporations, ${\dot{\rho }}_{\mathrm{PBH}}=\mathrm{const}.$ We also assume a uniform PBH density distribution in the vicinity of the Earth. A constant rate of evaporation implies that the derivative of the PBH density is related to the PBH temperature as

$$\begin{eqnarray}&&\displaystyle \frac{d{\rho }_{\mathrm{PBH}}}{{dT}}\propto {T}^{-4}.\end{eqnarray} \tag{ 9 }$$

The following steps were performed in the derivation of the PBH evaporation rate limit:

• 1.  
  A sample of PBHs $({T}_{i},{{\boldsymbol{x}}}_{i},{{\boldsymbol{v}}}_{i})$ was simulated with temperatures T_i > 5 GeV and T_i < 60 GeV distributed according to Equation (9), and distances R_i within R < 0.08 pc around the Earth. The velocities v_i of the sample PBHs were distributed with a mean equal to the orbital velocity of the Sun, v_rot = 250 km s⁻¹, and a dispersion v_disp = 270 km s⁻¹.
• 2.  
  For each PBH, we simulated the detection of the photons emitted over the 4 year 3FGL time period, consistent with the PBH evolution. The energies were distributed according to the instantaneous PBH spectrum of the appropriate temperature, and the positions of the photons were smeared according to the Fermi-LAT point-spread function (modeled as a Gaussian distribution). The emission spectra of PBHs Φ(E, t) are discussed in the Appendix. We used a time step of Δt = 1 day in modeling the evolution of the PBH position and temperature, and the number of photons detected by the Fermi-LAT each day was given by a Poisson random value with a mean of

  $$\begin{eqnarray}&&\overline{N(t)}=\displaystyle \frac{{\rm{\Delta }}t}{4\pi {R}^{2}}{\int }_{E=100\,\mathrm{MeV}}^{E=500\,\mathrm{GeV}}{\rm{\Phi }}(E,t)A(E){dE},\end{eqnarray} \tag{ 10 }$$

  where A is defined as the average Fermi-LAT exposure per unit time at the position of the simulated PBH and R is the distance from the Earth. The energy of each photon was found by random sampling of Φ(E, t) × A(E). Fermi-LAT has relatively uniform exposure on time periods longer than 1 day.
• 3.  
  The list of simulated PBH photons was concatenated to the real photons present within 5° of the final location of the PBH, with the same data selection as the 3FGL. A likelihood fit using the Fermi Science Tool gtlike was performed in a 7° × 7° ROI centered at the same location, using a model of the sky that included all 3FGL sources within 5° of the ROI center, as well as models of the isotropic diffuse and Galactic diffuse emission.⁶⁶ The PBH was modeled as a source with a LogParabola spectrum, with fitting parameters restricted to the ranges 1.2 < α < 3.0 and 0.0 < β < 1.0. Once the likelihood maximization was complete, the PBH was considered detected if its TS value was greater than 25, which is consistent with the 3FGL cutoff.
• 4.  
  If the PBH source was detected, the results from the likelihood fit were used to find the source flux in the five energy bins reported in the 3FGL catalog. The spectral consistency with a PBH spectrum was then calculated in the same way as described in Section 3.
• 5.  
  If the source was found to be spectrally consistent with a PBH, the significance of any proper motion was evaluated by the algorithm described in Section 3. The combined efficiency of steps 3–5 is displayed in Figure 3. We smoothed the results by convolving the detectability map with a 3 × 3 matrix of those that had a minor (≈8%) impact on the resulting limit. The impact of fluctuations was quantified by observing the change in the limit as the number of simulations increased; we found that an increase of the number of simulations by 100% had less than a 20% change in the resulting limit.
• 6.  
  To derive an upper limit on the number of PBH evaporations in our search region, we begin with the number of expected detections:

  $$\begin{eqnarray}&&N=\rho \epsilon V,\end{eqnarray} \tag{ 11 }$$

  where ρ is the true density of PBHs and V is the volume searched. is the average PBH detection efficiency in time t = 4 year and within the search volume V (a sphere with radius 0.08 pc, with the wedge corresponding to $| b| \lt 10^\circ$ removed); it is calculated by taking the mean over the pixels in Figure 3 with the weight R²T⁻⁴:

  $$\begin{eqnarray}&&\epsilon =\displaystyle \frac{\iint \epsilon (R,T)\tfrac{{R}^{2}}{{T}^{4}}\,{dR}\,{dT}}{\iint \tfrac{{R}^{2}}{{T}^{4}}\,{dR}\,{dT}},\end{eqnarray} \tag{ 12 }$$

  where the integrals run over the space of parameters described in step 1. Equation (11) can be inverted to find the PBH density ρ as a function of the number of detections N, or the upper limit on ρ given an upper limit on N. Given that one PBH candidate passed the selection criteria described in Section 3, we set an upper limit $N\lt 6.64$, which is the 99% confidence upper limit on the mean of a Poisson distribution with 1 observed event.
• 7.  
  We convert the upper limit on ρ to an upper limit on $\dot{\rho }$ by finding the fraction f of PBHs that would have evaporated during the search time t. Given a time of observation of 4 years, we find that all PBHs with initial temperature above 16.4 GeV would evaporate. Therefore,

  $$\begin{eqnarray}&&f=\displaystyle \frac{{\int }_{\text{16.4 GeV}}^{\text{60 GeV}}{T}^{-4}\,{dT}}{{\int }_{\text{5 GeV}}^{\text{60 GeV}}{T}^{-4}\,{dT}}.\end{eqnarray} \tag{ 13 }$$

  We calculate the 99% upper limit on ${\dot{\rho }}_{\mathrm{PBH}}$ to be:

  $$\begin{eqnarray}&&{\dot{\rho }}_{\mathrm{PBH}}\lt f\displaystyle \frac{6.64}{\epsilon {Vt}}=7.2\times {10}^{3}{\text{pc}}^{-3}\,{{\rm{}}{\rm{year}}}^{-1}.\end{eqnarray} \tag{ 14 }$$
• 8.  
  We estimated the systematic uncertainties arising from the uncertaintes in the PBH spectrum by varying the overall normalization of the PBH spectrum (see the Appendix) and varying the velocity distributions of the Milky Way disk and DM halo. We consider two scenarios ("aggressive" and "conservative") that give the best and worst sensitivity, respectively. Steps 1–7 are then repeated to find the resulting limit. The parameters of the aggressive and conservative models, as well as the resulting limits, are listed in Table 1. The limit including the systematic uncertainties is

  $$\begin{eqnarray}&&{\dot{\rho }}_{\mathrm{PBH}}\lt ({7.2}_{-2.4}^{+8.1})\times {10}^{3}{\text{pc}}^{-3}\,{{\rm{}}{\rm{yr}}}^{-1}.\end{eqnarray} \tag{ 15 }$$

Zoom In Zoom Out Reset image size

Table 1.  Parameters Used in the Estimation of Systematic Uncertainty

Model	Spectrum Normalization	Orbital Velocity (km s−1)	DM Halo Velocity (km s−1)	Limit
Aggressive	$\tfrac{0.45}{0.35}$	100	150	$4.8\times {10}^{3}\,{{\rm{}}{\rm{pc}}}^{-3}\,{{\rm{yr}}}^{-1}$
Conservative	$\tfrac{0.25}{0.35}$	300	350	$15.3\times {10}^{3}\,{{\rm{}}{\rm{pc}}}^{-3}\,{{\rm{}}{\rm{yr}}}^{-1}$

Note. To be more conservative in the estimates of the systematic uncertainties, we have tested the ranges of orbital velocities and the DM dispersion velocities that are larger than most of the values reported in the literature (Wilkinson & Evans 1999; Xue et al. 2008; Sofue et al. 2009; Gnedin et al. 2010; Kuhlen et al. 2010; McMillan & Binney 2010; McMillan 2011).

Download table as:  ASCIITypeset image

The potential existence of PBHs that emit detectable Hawking radiation is one of the most intriguing features of some theories of cosmological evolution. In addition to providing evidence for these theories, the possibility for direct observation of Hawking radiation (which would be a major discovery in its own right) was our main motivation for the search of PBHs with the Fermi-LAT.

In this paper, we searched for individual PBHs among 3FGL catalog sources. We performed calculations showing that the characteristic distance to a detectable PBH is of the order of ∼0.03 pc, in which case the proper motion of the PBH must be taken into account. We developed a new algorithm that can detect the proper motion of a PS in the presence of a known background. We found several 3FGL sources that have spectra consistent with PBHs, but none of these sources exhibit proper motion, which would be the smoking-gun signature of a PBH in the Fermi-LAT sensitivity domain. As a result, we derived upper limits on the local PBH evaporation rate.

To derive the efficiency of a PBH passing our selection criteria, we developed an MC that simulates PBHs with realistic velocity distribution and initial temperature distribution expected for the steady-state evaporation rate of the PBHs. The efficiency was then used to find the upper limits on the local PBH evaporation rate. The inferred systematic uncertainties are related to the uncertainties in the PBH gamma-ray spectrum as well as the uncertainties in the Galactic rotational velocity and the DM velocity dispersion. We calculated upper limits in scenarios for which these parameters covered a wide range of reasonable values.

In Figure 4, we compare the Fermi-LAT upper limit on the rate of PBH evaporations with the limits from Cherenkov telescopes and observe that they are similar. Although ground-based gamma-ray observatories are sensitive to timescales of a minute or less, and Fermi-LAT is sensitive to timescales of months to years, both the Cherenkov telescopes and the Fermi-LAT are probing the same quasi-stationary population of PBHs.

Zoom In Zoom Out Reset image size

The local evaporation rate can be related to the local mass density (relative to the critical density ρ_c) as (Halzen et al. 1991)

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{pbh}}^{\mathrm{loc}}=\displaystyle \frac{{M}_{* }}{(\alpha -2)}\displaystyle \frac{{\tau }_{0}}{{\rho }_{{\rm{c}}}}{\dot{\rho }}_{\mathrm{pbh}},\end{eqnarray} \tag{ 16 }$$

where α is the index of the initial distribution of PBH masses, dn/dM ∼ M^−α, and M_* is the mass of a PBH with a lifetime equal to the age of the universe. In the following, we will assume α = 2.5, which corresponds to PBHs formed in the radiation-dominated era (Halzen et al. 1991). The local density of PBHs corresponding to the evaporation rate in Equation (15) is ${{\rm{\Omega }}}_{\mathrm{pbh}}^{\mathrm{loc}}\leqslant ({3.6}_{-1.2}^{+4.1})\times {10}^{2}$. In general, the distribution of the PBH masses does not need to follow a power-law function. If there is a short period of low-pressure dust-like equation of state or a phase transition (see, e.g., the discussion in Carr et al. 2010 and references therein), then the PBH masses will be distributed around the mass in Equation (1), where t is the time of the PBH formation. In order for the PBHs to be evaporating now, we need M ∼ 10¹⁵ g, which corresponds to the formation time t ∼ 10⁻²³ s after the Big Bang (Carr et al. 2010). Thus, in general, limits on local PBH evaporation can constrain models where PBHs are formed with lifetimes close to that of the universe, i.e., the PBHs that are close to evaporation now. This requires certain "fine-tuning" of the time when the PBHs are formed. If, for example, the PBHs are formed close to the QCD phase transition with t_Q ∼ 10⁻⁴ s and T_Q ∼ 100 MeV, then the PBH masses are M ≳ 5M_⊙ ≈ 10³⁴ g (Dolgov & Blinnikov 2014). These PBHs have a temperature much smaller than 10⁻⁷ K and lifetime much longer than the lifetime of the universe, i.e., they cannot be detected by Fermi-LAT.

The local density of PBHs is expected to be enhanced compared to the average density in the universe in a similar way that the density of DM in the Galaxy is larger than the average density of DM in the universe. The enhancement factor for DM near the Sun (Bovy & Tremaine 2012) compared to the average DM density (Hinshaw et al. 2013) is k ∼ 2.2 × 10⁵. With this enhancement, the limit on the average PBH density is ${{\rm{\Omega }}}_{\mathrm{pbh}}\leqslant ({1.5}_{-0.5}^{+1.7})\times {10}^{-3}$. This limit is several orders of magnitude less constraining than the limits obtained from extragalactic and Galactic gamma-ray backgrounds ${{\rm{\Omega }}}_{\mathrm{pbh}}\,\leqslant {10}^{-8}-5\times {10}^{-10}$ (Carr & MacGibbon 1998; Lehoucq et al. 2009; Carr et al. 2010). The latter limits are calculated either by integrating the PBH evaporations inside the visible universe or inside the halo of our Galaxy, i.e., they are derived on kiloparsec to gigaparsec scales, while the limit in this paper is derived for distances less than a fraction of a parsec.

The limit on the average current density of PBHs can be translated to a limit on the density of PBHs at the time of formation, $\beta \sim {10}^{-18}{{\rm{\Omega }}}_{\mathrm{pbh}}\sqrt{M/{10}^{15}{\rm{g}}}$ (Carr 2005; Carr et al. 2010). This limit, in turn, can be used to constrain the spectrum of density fluctuations in the early universe (Josan & Green 2010; Linde et al. 2013; Torres-Lomas & Ureña-López 2013). In some cases, nonobservation of PBHs provides the only way to limit theories of inflation, especially the theories that predict large fluctuations of density at small distances (e.g., Linde et al. 2013).

The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.

The authors would like to thank Jane MacGibbon for numerous helpful discussions.

Facility: Fermi - Fermi Gamma-Ray Space Telescope (formerly GLAST).

Software: Astropy (http://www.astropy.org, Robitaille et al. 2013), matplotlib (Hunter 2007).

One of the most important factors in the systematic uncertainty of the limit on PBH evaporation is the spectrum of emitted gamma-rays. The spectrum of fundamental particles emitted by the black hole was computed by Hawking (1975). Quarks and gluons emitted by PBHs hadronize into mesons and baryons, which subsequently decay into stable particles. The spectra of the stable particles emitted by PBHs at high temperatures were computed by MacGibbon & Webber (1990).

To obtain the spectra of gamma-ray emission from PBHs, we use the values given in MacGibbon & Webber (1990) for T_PBH = 0.3, 1, 10, 50, 100 GeV and interpolate for different values of PBH temperatures. To cross-check numerical calculations, we use two different interpolation methods:

• 1.  
  We use an analytic approximation by fitting the PBH emission rate ${\dot{N}}_{\gamma }(E,T)$ with a cubic log polynomial

  $$\begin{eqnarray*}&&\mathrm{log}\dot{N}={c}_{0}(T)+{c}_{1}(T)\mathrm{log}x+{c}_{2}(T){(\mathrm{log}x)}^{2}+{c}_{3}(T){(\mathrm{log}x)}^{3},\end{eqnarray*}$$

  where x = E/T and interpolate the fit coefficients c_i(T) (we use this appriximation in Section 2).
• 2.  
  We create a table for a set of E and T values and use a two-dimensional interpolation directly from the results in MacGibbon & Webber (1990; this approximation is used in Sections 3 and 4).

We compare the first interpolation with other parametrizations available in the literature (Halzen et al. 1991; Ukwatta et al. 2013) in Figure 5 left. There seems to be a rather significant discrepancy in the total energy emitted in gamma-rays. For instance, MacGibbon & Webber (1990) found that 24%–25% of the energy is emitted in gamma-rays for a large range of initial temperatures, while analytical integration of the parametrization (Bugaev et al. 2008; Petkov et al. 2008)

$$\begin{eqnarray}\displaystyle \frac{{{dN}}_{\gamma }}{{{dE}}_{\gamma }}=9\times {10}^{35}\left\{\begin{array}{ll}{\left(\tfrac{1\mathrm{GeV}}{{T}_{\tau }}\right)}^{3/2}{\left(\tfrac{1\mathrm{GeV}}{{E}_{\gamma }}\right)}^{3/2}{\mathrm{GeV}}^{-1} & {E}_{\gamma }\lt {T}_{\tau };\\ {\left(\tfrac{1\mathrm{GeV}}{{E}_{\gamma }}\right)}^{3}{\mathrm{GeV}}^{-1} & {E}_{\gamma }\geqslant {T}_{\tau }\end{array}\right.\end{eqnarray} \tag{ 17 }$$

gives about 47% of the energy in gamma-rays. A more recent parametrization of the gamma-ray spectra at temperatures ≳1 TeV (Figure 11 of Ukwatta et al. 2016) gives about 35% of the energy in gamma-rays.

For our baseline model, we take the spectrum of MacGibbon & Webber (1990) rescaled to give 35% of energy in gamma-rays to match the more recent calculation in Ukwatta et al. (2016). For the calculations, we use the interpolation presented in Figure 5 (right). At low temperatures, the "Interpolation" curve is calculated by taking the instantaneous PBH spectrum from MacGibbon & Webber (1990) rescaled to 35% energy going to gamma-rays times 4 years. At temperatures above the temperature of a PBH with the 4 year lifetime, the "Interpolation" curve is the Bugaev et al. (2008) parametrization rescaled by 0.45. There is a good agreement between the integrated spectra of MacGibbon & Webber (1990; rescaled to 35% efficiency) and the "Interpolation" curve.

Zoom In Zoom Out Reset image size