ournal of C osmology and A stroparticle hysics J Constraining the local burst rate density of primordial black holes with HAWC

. Primordial Black Holes (PBHs) may have been created by density ﬂuctuations in the early Universe and could be as massive as > 10 9 solar masses or as small as the Planck mass. It has been postulated that a black hole has a temperature inversely-proportional to its mass and will thermally emit all species of fundamental particles via Hawking Radiation. PBHs with initial masses of ∼ 5 × 10 14 g (approximately one gigaton) should be expiring today with bursts of high-energy gamma radiation in the GeV-TeV energy range. The High Altitude Water Cherenkov (HAWC) Observatory is sensitive to gamma rays with energies of ∼ 300 GeV to past 100 TeV, which corresponds to the high end of the PBH burst spectrum. With its large instantaneous ﬁeld-of-view of ∼ 2 sr and a duty cycle over 95%, the HAWC Observatory is well suited to perform an all-sky search for PBH bursts. We conducted a search using 959 days of HAWC data and exclude the local PBH burst rate density above 3400 pc − 3 yr − 1 at 99% conﬁdence, the strongest limit on the local PBH burst rate density from any existing electromagnetic measurement.


Introduction
While there are no known processes in the current Universe that can create black holes with masses less than ∼ 1 M , conditions in the early Universe were conducive to the formation of black holes with a wide range of masses [1]. These black holes, with masses ranging from the Planck mass to supermassive black holes (> 10 9 M ), are called Primordial Black Holes (PBHs). PBH production in the early Universe would have broad observable consequences spanning the largest distance scales (including influencing the development of large-scale structure in the Universe and the primordial power spectrum [2][3][4]), to the smallest scales (including enhancing local dark matter clustering [5,6]). In the present Universe, PBHs in certain mass ranges may constitute a non-negligible fraction of dark matter [1,7,8]. Since the existence of stellar-mass black holes was recently confirmed during the first observational run of Advanced LIGO [9], there has been a resurgence in support for a PBH component of the total dark matter energy density (e.g., refs. [3,10,11]). Limits placed thus far indicate that f (m), the fraction of dark matter that is made up of PBHs, is 10% over a range of masses [8].
The prediction that a black hole will thermally radiate, or 'evaporate,' with a blackbody temperature inversely proportional to its mass was first developed by Hawking in a calculation that convolved quantum field theory, General Relativity, and thermodynamics [12]. The emitted radiation consists of all fundamental particles with masses less than approximately the black hole temperature [13]. For black holes in the stellar mass range and above, Hawking radiation is nearly negligible. However, for lower-mass PBHs, this process dominates their evolution over time [14]. PBHs with initial masses of ∼ 5 × 10 14 g should be expiring today producing short bursts lasting a few seconds of high-energy gamma radiation in the GeV-TeV energy range [15,16], making their final moments an ideal phenomenon to observe with HAWC.
Confirmed detection of a PBH burst -beyond proving their existence and allowing the determination of their relic density and rate-density of evaporation -would provide valuable insights into many areas of physics, including fundamental processes in the very early Universe and particle physics at energies higher than currently achievable by terrestrial -1 -

JCAP04(2020)026
accelerators [17]. Even the non-detection of PBH burst events in dedicated searches would yield important constraints about the early Universe [1]. One of the most important reasons to search for PBHs is to constrain the cosmological density fluctuation spectrum in the early Universe on scales smaller than those constrained by the cosmic microwave background [18]. A particularly interesting question is whether or not PBHs were formed from the quantum fluctuations associated with many different types of inflationary scenarios [1]. Detection or upper limits on the number density of PBHs can thus also inform inflationary models.
Numerous detectors have searched for PBH burst events using both direct (e.g., [14,19,20]) and indirect (e.g., [1,21,22]) methods. These methods explore the PBH distribution at various distance scales, from cosmological scales down to within a parsec. At cosmological distances, the PBH density can be probed using the 100 MeV extragalactic gamma-ray background [1,12,15,[23][24][25][26]. On the Galactic scale, the local PBH density and anisotropy can be studied using the 100 MeV gamma-ray measurements [21]; on kiloparsec scales, a PBH burst limit can be placed using the Galactic antiproton background [22]; and on the parsec scale, the PBH burst limits can be set directly by searching for the detection of individual evaporating PBHs [14,19,20,27]. In this work (an extension of the analysis presented at the 36th International Cosmic Ray Conference in Madison, WI, U.S.A. [28]), the High Altitude Water Cherenkov (HAWC) Observatory searches for PBH bursts at the parsec scale.
This paper is structured as follows. The HAWC Observatory is described in section 2. Section 3 provides a description of the model of PBH bursts. The data collection and analysis procedure is presented in section 4. Section 5 provides a discussion of the statistical and systematic uncertainties and concludes the paper.

HAWC observatory
The HAWC Observatory, a successor to the Milagro Observatory [29], is a very-high-energy (VHE) ground-based air shower array located on the side of the Sierra Negra volcano in Mexico at an altitude of 4,100 m above sea level. It has a wide field-of-view of ∼ 2 sr and an operational energy range of ∼ 300 GeV-100 TeV. HAWC consists of 300 cylindrical water tanks in the main array covering a total area of 22,000 m 2 . Each tank in the main array is 7.3 m in diameter and 4.5 m deep, and is equipped with four (three 8 and one 10 in diameter) upward-facing photomultiplier tubes (PMTs) anchored to the bottom of the tank. HAWC has completed installing and integrating an additional "outrigger" array [30] composed of 345 cylindrical tanks 1.55 m in diameter and 1.65 m deep, each containing a single 8 PMT. The outriggers are arranged in a concentric, circular, symmetric pattern around the main array, covering an additional instrumented area of ∼ 4.5 times that of the main array. However, the analysis presented in this work includes only data from the main HAWC array.
In both the main array and the outriggers, the PMTs detect Cherenkov light from secondary particles created in extensive air showers induced by VHE gamma rays incident on Earth's atmosphere. The main data acquisition system measures the arrival time of secondary particles on the ground and the amplitude of the PMT signals. This information is used to reconstruct the arrival direction and energy of the primary particle. The angular resolution of HAWC ranges from ∼ 0.2 • (68% containment) to 1.0 • , depending on the fraction of PMTs hit by the resulting shower [31]. With these features and a high duty cycle of greater than 95%, HAWC is ideally suited to continuously monitor the Northern Hemisphere sky for high-energy emission from gamma-ray transients such as PBHs. HAWC's wide field-of-view and continuous operations are advantageous for detecting burst transients with emission durations shorter than the slewing times of Imaging Atmospheric Cherenkov telescopes (IACTs). These features are also key for this analysis as the sensitivity to a PBH burst rate density is determined by the total observable volume and the exposure time.

Primordial black hole burst spectrum
As a PBH radiates, it continually loses mass and its temperature increases to very high energies [12]. The manner in which the PBH expires depends on the physics at this energy scale. The chosen high-energy particle physics model determines the energy spectrum. In this work, as in ref. [27], we assume the Standard Evaporation Model (SEM) [13,32] as our emission and particle physics model. In the SEM -based on the Standard Model of particle physics -a PBH should directly radiate those fundamental particles whose wavelengths (Compton wavelength λ c = h/mc for massive particles) are comparable to the size of the black hole. When the black hole temperature exceeds the Quantum Chromodynamics confinement scale of ∼250-300 MeV, quarks and gluons should be radiated [13,26,33]. On astrophysical timescales, the final products will decay to the lightest Standard Model particles: photons, neutrinos, electrons, positrons, protons, and anti-protons [13].
In the SEM, the black hole temperature (T ) can be expressed in terms of the remaining lifetime (τ ) of the black hole (that is, the time left until the black hole finishes evaporating) as [32,34], for T T p , where T p is the Planck temperature (1.22 × 10 19 GeV). Note that here we are using units such that the Boltzmann constant k = 1. The emission rate increases as the black hole shrinks [26]. For black holes with temperatures greater than several GeV at the start of the observation, the time-integrated photon flux can be parameterized as [18,34,35], for gamma-ray energies E 10 GeV. This parameterization includes both directly radiated photons and those produced by the decay of other directly radiated species. Figure 1 shows the total gamma-ray spectrum for the PBH remaining lifetimes (τ ) examined in this work, τ = 0.2 s, 1 s, and 10 s.

Creating the model
The first step in modeling PBH bursts is to simulate burst source points in HAWC's fieldof-view. We used a Monte Carlo simulation to generate a random set of PBH burst events assuming a local burst rate density,ρ, of 10 4 pc −3 yr −1 . We chose this burst rate density merely because it is close to HAWC's predicted sensitivity [27], but it could have been set to be any positive non-zero value. The events were distributed uniformly in a 50 • cone centered in the HAWC field-of-view, out to a distance of 0.5 pc. Beyond 0.5 pc, even PBH bursts simulated for 10 s at HAWC's zenith did not produce an observable signal within HAWC. The parameterization of the time-integrated photon flux given in eq. (3.2) can be used to calculate the expected number of photons detectable by an observatory on the Earth's surface. For a PBH burst of duration τ at a non-cosmological distance r and zenith angle θ, the number of expected photons may be expressed as, where dN/dE is the PBH gamma-ray spectrum integrated from remaining lifetime τ to 0. The energies E 1 and E 2 correspond to the energy range of the detector, A(E, θ) is the effective area of the detector as a function of photon energy and zenith angle, and f is the dead time fraction of the detector. While eq. (3.3) schematically illustrates how to approach this portion of the analysis, the actual calculations used a forward-folding approach [36]described below -to better characterize the detector.
To determine the expected signal at HAWC from each simulated PBH in the fieldof-view, we utilized a forward-folding approach based on a functionality within the HAWC software called ZEBRA, which stands for ZEnith Band Response Analysis [37]. ZEBRA uses simulation to characterize the response of the HAWC detector as a function of zenith angle, which is then convolved with the expected spectrum of the source to estimate the counts observed from that source during an arbitrary period of time. Our source model assumed the timescale of the emission was short enough to be considered as all coming from a single zenith angle, with the PBH spectrum following eq. (3.2) [18]. This procedure was completed for each of the remaining lifetimes we consider. The number of estimated counts given by ZEBRA for a simulated PBH with remaining lifetime of 0.2 s at a distance of 0.05 pc is shown in figure 2. The sample PBH is assumed to have a remaining lifetime of 0.2 s and be located at a distance of 0.05 pc. The slope plateaus above 45 • due to HAWC's low sensitivity at those zenith angles.
Using the estimated signal, we calculated the Poisson probability (p-value) of obtaining N or more counts given the background, B, for each burst event,

HAWC transient search data
The data set used for the PBH burst search, H data , consists of data from a self-triggered all-sky transient search conducted from March 2015 through May 2018 [38]. This program, originally designed for use in a gamma-ray burst (GRB) analysis, continuously searches for transients at energies above a few hundred GeV with sliding time windows of lengths 0.2, 1, and 10 seconds -corresponding to typical timescales of peak structures within GRBs. It detects clusters of events above the cosmic-ray rate and stores the probability of observing each cluster of events under a background-only hypothesis. Events were required to have fired at least 70 PMTs and pass a loose hadron-rejection cut optimized for gamma-ray energies < 10 TeV. The background, used herein to form H bkg , is estimated using a 1.75 hour  integration [38]. To be computationally efficient, these p-values are saved only if they fall below a certain value. This feature can be seen in the righthand edge of H data in figure 3.
The search is performed by shifting each window forward in time and binning air shower events during that window using a grid of 2.1 • × 2.1 • square spatial bins covering all points within 50 • of detector zenith. The step size of each shift is less than the bin size to allow for desired overlap. The size of the overlap is optimized to allow for fine tuning on the spatial position of the events while avoiding strong correlations between pixels, and is a different amount for each search duration. Points outside a zenith angle of 50 • are excluded from the spatial search because most photons at the energies expected from a transient burst signal with θ > 50 • do not have sufficient energy to reach HAWC due to attenuation of air showers at larger atmospheric depth. The fixed-width time windows and square bins were utilized -rather than attempting to fit a light curve profile (which would improve sensitivity) -to make a full search of the HAWC field-of-view computationally tractable [38]. The analysis threshold for each remaining lifetime searched in this PBH analysis was chosen to be the -6 -JCAP04(2020)026 p-value where there are 10 events in the background histogram to ensure the limit placed is conservative and not sensitive to fluctuations in the data. A vertical black line indicates this threshold in figure 3. Above this analysis threshold, where p > p thr , we have a fiducial region in which we can verify -independently of the analysis -that the background and data are statistically equivalent; that is, to ensure that the overall rate normalization seen in data matches the background estimate from simulation. As we can be sure the data are welldescribed by a background-only distribution above p thr , we are confident in our extrapolation of the background-only model to the low-statistics area.

Calculating an upper limit
Finding no statistically significant PBH signal in the data, we computed upper limits on the local burst rate density of PBHs. We began by choosing an approximate value for the local burst rate density of PBHs,ρ ∼ 10 4 pc −3 yr −1 -the PBH rate thrown in our Monte Carlo simulation as described in section 3.2. We then scanned overρ until we found the most probable value ofρ given the data, as well as the 99% upper limit.
To determine the 99% limit, we began by calculating the test statistic forρ, where L 0 is the background Poisson likelihood and L 1 is our model Poisson likelihood. The TS follows a χ 2 distribution with one degree-of-freedom based on Wilks' theorem [39]. In general for binned data, the log-likelihood of a Poissonian variable x, with mean λ for an independently and identically distributed sample of size n, is given by, such that we can write, summing over bins, p, instead of sample size, where H data is a histogram of the HAWC data p-values produced by HAWC's transient burst search [38,40] and H bkg is a histogram of p-values generated by Monte Carlo using HAWC background distributions, also from the transient burst search. Similarly, the model log-likelihood can be written as, where H model (p) = H bkg (p) + H PBH (p). Note that both expressions neglect the factorial term in the likelihood eq. (4.2) as it will cancel when evaluating eq. (4.1). We then found the value ofρ that yielded the largest T S value from eq. (4.1), T S max . Scanning over increasing values ofρ, we stopped when the change in T S from the maximum reached ∆T S = 5.41, which corresponds to a one-sided 99% confidence interval. Although our strictest limit was placed for a remaining lifetime of 0.2 s, in an effort to place a conservative limit we report the remaining lifetime for which HAWC's sensitivity was predicted to be the strongest -corresponding to τ = 10 s, ρ < 3400 +400 −100 pc −3 yr −1 , (4.5) the strictest limit yet placed on the local PBH burst rate density. The uncertainties in the limit are systematic only and are described in section 5.1.

Discussion
The 99% upper limits on the PBH burst rate density for each of the three remaining lifetimes searched are listed in table 1 and shown in figure 4. The expected limits, as well as the 68% and 95% containment for the null hypothesis, are also shown in figure 4. The containment bands and expected limits are calculated using 1000 simulations containing no PBHs. Note that the bands are purely statistical and indicate that the upper limits are compatible with the expected sensitivity of HAWC assuming only background events. Our reported limit, corresponding to a remaining lifetime of 10 s, is also presented in table 2 for direct comparison with the results of other direct searches for PBHs. All limits shown in figure 4 and tables 1 & 2 were obtained based on the PBH Standard Emission Model [13,25].

Systematic uncertainties
The main source of systematic uncertainties within HAWC analyses comes from discrepancies between the data and the simulated Monte Carlo events, which stem from uncertainties in the modeling of the detector. The dominant sources of these discrepancies are PMT efficiency and late light simulation. The PMT efficiency cannot be precisely determined using the calibration system; instead, an event selection based on charge and timing cuts is implemented to identify incident vertical muons. Since vertical muons provide a mono-energetic source of light, they can be used to measure the relative efficiency of each PMT by matching the muon peak position to the expected position from simulations, which is used to measure the range of uncertainties. Uncertainty also arises from mis-modeling of late light in air showers stemming from a discrepancy between the time width of the laser pulse used for calibration and the time structure of actual showers. Less dominant contributors include PMT thresholds, angular resolution discrepancy, and charge uncertainty. The PMT threshold in simulations may be, based on observations of the cosmic-ray rate, deviating ±0.05 PE from the actual value; the 68% containment radius of the reconstructed gamma-ray direction around the true direction in the HAWC   [41], CYGNUS [42], the Tibet Air Shower Array [43], H.E.S.S. [14], VER-ITAS [19], Fermi -LAT [20], and Milagro [27]. Also shown are the expected limits, 68% containment, and 95% containment for the null hypothesis based on HAWC sensitivity. Results displayed between the three explicitly evaluated burst durations are interpolations. While the statistical fluctuations of the limits show possible turn-over in the inferred limit, the sensitivity line indicates a smooth trend toward worse limits for shorter remaining lifetimes. detector Monte Carlo model is underestimated by ∼ 5%; and the charge uncertainty stems from relative differences in photon detection efficiency from PMT to PMT, encapsulating how much of a PMT measurement will vary for a fixed amount of light. These effects are described in more detail in ref. [44] and have been evaluated for correlations with none found. To account for these uncertainties in this analysis, we made a new model histogram for each source of uncertainty based off of simulations using different detector models and then repeated the analysis. We assume these variables are independent and add each source of systematic uncertainty in quadrature with the others, resulting in approximately +10.6%, −4.3% uncertainty in our results (shown in figure 5). These uncertainties, being significantly smaller than the 1σ and 2σ containment bands of the expected limit, are not shown in figure 4.

Conclusions
We evaluated three years of HAWC transient search data for PBH bursts. No significant signal was found, so we used the log-likelihood method described in section 4.2 to set upper limits on the local burst rate density at 99% confidence. We set a limit of ∼ 3400 pc −3 yr −1 using a burst duration of 10 s, the most constraining limit placed to date for very nearby PBHs. Note that the burst duration is a search parameter, not a physical parameter, thus differences between the limits placed for each burst duration searched are due to differences in the signal-to-background ratio within HAWC. Planned future work to improve this analysis includes a dedicated PBH search working directly with HAWC data rather than results from a previous blind transient search on HAWC data. This would also incorporate more durations to search, energy estimators, and more days of HAWC data. Statistical improvements for such a study would include full likelihood profiles in both time and space, as well as stacking of these likelihoods to ensure the signal is well-defined for any zenith angle. Also, data using the newly-built outrigger array will increase the sensitivity at > 10 TeV [30], relevant for hard-spectrum sources such as PBHs. We anticipate significant enhancement in our ability to search for PBH bursts with the application of such improvements.