Constraining the Primordial Black Holes as Dark Matter at JUNO

As an attractive candidate for dark matter, the primordial black holes (PBHs) in the mass range ($10^{15} \sim 10^{16}$)$\mathrm{g}$ could be detected via their Hawking radiation, including neutrinos and antineutrinos of three flavors. In this paper, we investigate the possibility to constrain the PBH as dark matter by measuring (anti)neutrino signals at the large liquid-scintillator detector of Jiangmen Underground Neutrino Observatory (JUNO). Among six available detection channels, the inverse beta decay $\overline{\nu}^{}_e + p \to e^+ + n$ is shown to be most sensitive to the fraction $f^{}_{\rm PBH}$ of PBHs contributing to the dark matter abundance. Given the PBH mass $M^{}_{\rm PBH} = 10^{15}~{\rm g}$, we find that JUNO will be able to place an upper bound $f^{}_{\rm PBH} \lesssim 3\times 10^{-5}$, which is twenty times better than the current best limit $f^{}_{\rm PBH} \lesssim 6\times 10^{-4}$ from Super-Kamiokande. For heavier PBHs with a lower Hawking temperature, the (anti)neutrinos become less energetic, leading to a relatively weaker bound.


I. INTRODUCTION
The observations of gravitational waves (GWs) by Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) [1] have stimulated very active discussions on the origin of binary black holes (BBHs). In particular, great interest in the PBHs as the seeds for the formation of BBHs has been revived [2][3][4][5][6][7][8][9]. Recent studies [2,10,11] have shown that the PBH scenario predicts a merger rate that is well consistent with the local observations. On the other hand, produced in the early Universe, PBHs could be a viable candidate for cold dark matter [12][13][14]. The fraction of dark matter in the form of PBHs is usually defined as f PBH = Ω PBH /Ω DM , where Ω PBH and Ω DM denote the energy density fraction of PBHs and that of dark matter in the present Universe, respectively. Various constraints on f PBH within a broad range of PBH masses have been derived in the literature. See recent reviews in Refs. [15,16] and references therein. The Hawking radiation [17] of the PBHs has been suggested for experimental observations and can be * wangsai@ihep.ac.cn † xiadm@cqu.edu.cn ‡ zhangxukun@ihep.ac.cn § zhoush@ihep.ac.cn ¶ changz@ihep.ac.cn used to explore the intrinsic properties of PBHs [18][19][20][21][22][23][24]. As is well known, the black holes could emit particles near their event horizons, and the black hole evaporates faster as its mass decreases [25][26][27]. For a PBH, depending on its formation time, it could have a much smaller mass than ordinary stars [12,13]. The PBHs with masses M PBH O(10 15 ) g would have evaporated over at the present age of the Universe [17]. Based on the experimental searches for the Hawking radiation, stringent constraints on f PBH within a mass range (10 14 ∼ 10 17 ) g have been obtained (See, e.g., Refs. [15,16], for recent reviews). For example, there exists an upper limit f PBH 10 −7 for M PBH 10 15 g through an observation of the e ± flux from Voyager 1 [18]. By observing the γ-ray lines from the e − -e + annihilation and the isotropic diffuse γ-ray background, respectively, International Gamma-Ray Astrophysics Laboratory (INTEGRAL) [28] and Fermi Large Area Telescope (FERMI-LAT) [29] have restricted f PBH to be less than 10 −3 for M PBH = 10 16 g [19,20]. In Ref. [19], it has also been shown that an upper limit f PBH 6 × 10 −4 for M PBH 10 15 g can be derived via the neutrino observation from Super-Kamiokande neutrino observatory (Super-K) [30].
With a perfect neutron-tagging efficiency of the LS detector, JUNO is expected to have a great potential to discover the diffuse supernova neutrino background (DSNB) [31]. As mentioned above, the upper limit on the DSNB flux from Super-K has been used to constrain f PBH for small-mass PBHs [19]. Therefore, we are well motivated to investigate how restrictive the bound on f PBH from JUNO will be. Compared to Super-K, we expect that JUNO could improve the constraint on f PBH due to its more powerful neutron tagging, lower energy threshold and better energy resolution. A similar analysis can be performed for the future Gadolinium-doping upgrade of Super-K.
The remaining part of this paper is organized as follows. In Sec. II, we calculate the neutrino fluxes from PBHs in the Galactic and extragalactic dark matter halos. Then, Sec. III is devoted to the event spectra for all six neutrino detection channels at JUNO and the estimation of the relevant backgrounds. In Sec. IV, by comparing between the signals and backgrounds, we draw the upper limit on f PBH from JUNO. Our main conclusions are finally summarized in Sec. V.

II. NEUTRINO FLUXES FROM PBHs
Generally speaking, there are two different contributions, i.e., the primary and secondary components, to the neutrino fluxes from the evaporating PBHs. The former arises directly from the Hawking radiation, while the latter stems from the decays of the secondary particles produced in the Hawking radiation [21,32]. For an evaporating PBH, the number of (anti)neutrinos (N ) in unit energy (E) and time (t) is given by where the first and second terms on the right-hand side denote the primary and secondary components, respectively. In our calculations, both components are evaluated by using blackhawk [32]. It is worthwhile to notice that the spins of PBHs in the present work are assumed to be negligible, as predicted by a class of theories for the PBH production [33,34]. Then the differential fluxes of (anti)neutrinos radiated from PBHs can be calculated by taking into account of the distribution and cosmological evolution. In units of cm −2 s −1 MeV −1 , they can be decomposed as follows where the first and second terms at the right-hand side are contributed by PBHs located in the Galactic and extragalactic dark halos, respectively. For the extragalactic PBHs, the differential (anti)neutrino flux can be written as [19,21] According to the Planck 2018 results [35], we set the average energy density of dark matter in the present Universe to be ρ DM = 2.35 × 10 −30 g cm −3 . The (anti)neutrino energy at the source has been denoted as E s , while that in the observer's frame as E. The connection between them is given by E s = [1+z(t)]E, where the redshift z(t) is a function of the cosmic time t and encodes the cosmological evolution. For the lower and upper limits of the integration over the cosmic time in Eq. (3), we take the following values. First, the (anti)neutrinos emitted from the PBHs in the very early time will be significantly redshifted such that they with extremely-low energies today cannot be detected. There is an exception for high-energy (anti)neutrinos, whose fluxes however are highly suppressed for the PBH masses of our interest. Therefore, we fix t min = 10 11 s close to the epoch of radiation-matter equality and numerically confirm that changing t min to be smaller has essentially no impact on the final results. Second, we choose t max = min{τ 0 , τ PBH }, where τ 0 is the age of the Universe, τ PBH is the lifetime of PBHs, and the function min{x 1 , x 2 } singles out the smaller one from x 1 and x 2 . For the PBHs in the Galactic dark halo, the differential (anti)neutrino flux can be calculated as [19] dF where r(l, ψ) ≡ r 2 − 2lr cos ψ + l 2 is the galactocentric distance calculated from the distance of the Earth to the Galactic center r = 8.5 kpc and the line-of-sight distance l to the PBH, with ψ being the angle between these two directions. In addition, the angular integration is defined as dΩ = 2π 0 dφ π 0 dψ sin ψ with φ being the azimuthal angle, and ρ MW (r) is the local energy density of dark matter. The maximal value of l is determined by l max = (r 2 h − r 2 sin 2 ψ) 1/2 + r cos ψ with the halo radius r h = 200 kpc. For illustration, we implement the Navarro-Frenk-White (NFW) profile [36] ρ NFW (r) = 0.4 × 8.5 r which is in units of GeV · cm −3 , and the distance r is in units of kpc. In Fig. 1, we show the differential fluxes of ν e emitted from PBHs for four benchmark masses within the range (10 15 ∼ 10 16 ) g. The results for neutrinos and antineutrinos of other flavors are quite similar, as a consequence of the thermal Hawking radiation. In the numerical calculations, we have assumed a monochromatic mass function of PBHs and taken the existing upper limit of f PBH from Super-K [19]. As one can observe from Fig. 1, neutrinos and antineutrinos from smaller-mass PBHs have higher energies, mainly due to the higher Hawking temperature, while the magnitude of the fluxes is limited by the existing bound on f PBH . However, at low energies, the fluxes converge to a similar value and their differences become small.

III. EVENT SPECTRA AND BACKGROUNDS
Given the (anti)neutrino fluxes radiated by PBHs, we can figure out the event spectra for six available detection channels in the LS detector of JUNO. In this section, we present the results of the event spectra and estimate the relevant backgrounds at JUNO.

A. IBD channel
For the IBD channel ν e + p → e + + n, both the final-state e + and n can be perfectly observed in the LS detector. At JUNO, the event spectrum of IBD is given by [37] where the energy threshold for the IBD reaction is E thr is a Gaussian function of E 0 with an expectation value E v and the standard deviation δ E = 3%/ E 0 /MeV, i.e., the energy resolution of JUNO [31,37]. The visible energy in the detector E v = m e + E e + arises from the annihilation of the final-state positrons with ambient electrons, where m e = 0.511 MeV is the electron mass. The IBD cross section σ IBD νe and the positron energy E e + de- pend on the antineutrino energy E in the differential flux dF/dE defined in Eq. (2). We obtain the dependency relation in Tab. 1 of Ref. [38]. The total number N p of target protons in the LS and the effective running time T in Eq. (6) can be found in Ref. [37].
With the differential fluxes of ν e in Fig. 1, we compute the IBD event spectra by using Eq. (6) and present the final results in Fig. 2, where different curves correspond to four benchmark PBH masses in Fig. 1. In our calculations, the fiducial mass of the JUNO detector is taken to be 20 kiloton and the operation time is set to T = 10 yrs. In Fig. 2, we can observe that the peak of the event spectrum moves towards lower energies as the PBH mass increases. It should be noticed that the peak rate is limited by the existing upper bound on f PBH .
The main backgrounds for the IBD channel can be divided into four categories. The first one is the irreducible background from reactor antineutrinos, which dominate over others in the energy range below 12 MeV [31,39]. From Fig. 2, we notice that when the PBH mass is 8 × 10 15 g or larger, the peaks of the event spectra shift to the region E 0 10 MeV.
Consequently, for M PBH 8 × 10 15 g, the IBD signal will be severely contaminated by the reactor antineutrino background. For this reason, we mainly focus on the PBHs with smaller masses or set an energy cut E 0 12 MeV to get rid of this background. The second one is DSNB, which is one the primary goals of Super-K and JUNO [31]. For illustration, we take into account the DSNB ν e flux with an average energy Eν e = 15 MeV. The third type of backgrounds arises from the atmospheric neutrino charged-current (atm.CC) background caused by the atmosphericν e , and the atmospheric neutrino neutral-current (atm.NC) background by the neutral-current (NC) reactions of high-energy atmospheric (anti)neutrinos in the liquid scintillator (LS) [31,40]. Both kinds of event spectra have been studied in Ref. [31]. The fourth background is the fast neutrons (FN), which are produced by the cosmic high-energy muons arriving at the detector and induce the IBD-like events [31]. For the last two categories of backgrounds, the dedicated method of pulse-shape discrimination (PSD) has to be utilized, as we shall explain below.
Similar to the DSNB searches in the LS detectors, the PSD approach is critically important to enhance the signal-to-noise ratios (SNRs). Both signals and backgrounds in the LS will be finally converted into light, which can be observed by the photomultiplier tubes (PMTs). For the light signals generated in different processes, the probability density functions of the photon emission time can be described as [41] where i refers to the fast, slow and slower components, and N i and τ i denote the fraction and time constant of i-th component, respectively. Then one can in principle identify different particles by the corresponding distribution functions F (t). To make use of the tail-to-total ratio [41], we define the working parameter PSD as where t cut denotes a cutoff time. For the signal and one kind of the relevant backgrounds, one can prop- erly choose the cutoff time in order to reduce the number of background events more than that of the signal events. For JUNO, Ref. [31] has estimated the PSD efficiency, which is the fraction of events surviving the PSD procedure, ν = 0.5, NC = 0.011, and FN = 0.017.
In Table I, we summarize the total number of IBD events from PBHs and that of relevant backgrounds, where both the results before and after applying the PSD approach have been shown for comparison. The IBD events from PBHs are calculated as in Fig. 2, and the backgrounds are categorized as in previous discussions. Comparing the total backgrounds before and after applying the PSD approach, we find a significant reduction of the backgrounds. One can notice that the dominant background is due to atm.NC, and for M PBH > 10 15 g, the SNRs without the PSD procedure are too small for any realistic observations.

B. eES channel
For the eES channel ν + e − → ν + e − , all species of neutrinos and antineutrinos could contribute and the event spectrum can be calculated as [37]  where the eES cross section is σ eES ν α with α standing for neutrinos and antineutrinos of three flavors, and T e is the electron kinetic energy. The minimal neutrino energy to produce a recoil electron energy T e is determined by E min e T e /2 + T e (T e + 2m e )/2. In addition, we should set N e = N p due to the electric neutrality of the LS materials.
The eES event spectra from PBHs are shown in Fig. 3, where four benchmark masses for the PBHs corresponding to those in Fig. 2 have been taken and the effective running time of 10 years and the fiducial LS mass of 20 kilotons have been assumed. As one can observe from Fig. 3, the maximal rate reaches 1.5 MeV −1 for M PBH = 8 × 10 15 g at the low-energy end. However, for the electron recoil events, there are numerous backgrounds from the 8 B solar neutrino and other three types of cosmogenic isotopes (i.e., 11 C, 10 C and 11 Be). According to Ref. [31], the total background events are estimated to 10 4 MeV −1 per 10 years per 20 kilotons at JUNO. Compared with this background rate, the eES event rates shown in Fig. 3 are much smaller. Therefore, it seems impossible to observe the eES signals from PBHs at the LS detectors.

C. pES channel
As the advantages of the low energy threshold of the LS detectors, the pES channel ν + p → ν + p can also be used to observe astrophysical neutrinos. The pES event spectrum is computed as follows [37] where dσ pES ν α /dT p denotes the differential cross section with T p being the proton recoil energy. The minimal neutrino energy required to produce the finalstate proton with a recoil energy of T p is approximately determined by E min p = (T p m p /2) 1/2 with m p = 938.27 MeV being the proton mass. However, the proton recoil energy will be quenched in the LS and observed as T p , which is related to T p by the Birks law [37,42]. The quenching effects on recoiled protons in the LS have been shown numerically in Fig. 30 of Ref. [31] and will be taken into account in our calculations.
Since the recoil energies of protons are small and the observed energies after quenching effects become even smaller, we consider only the (anti)neutrinos from the PBHs with a mass of M PBH = 10 15 g, for which the (anti)neutrino energies are much higher than those for larger PBH masses. The final result has been presented in Fig. 4, from which one can see that the observed energy turns out to be located in the range of (0.2 ∼ 1.8) MeV. Below 1 MeV, the background is mainly caused by the radioactive decays of the LS materials and surroundings, whose rate has been estimated to be about 80 events per 10 seconds [31]. In contrast, the pES event rate produced by PBHs is small, as one can observe from Fig. 4. It should be noticed that the operation time for the LS in Fig. 4 has been taken to be 10 years, together with the fiducial LS mass of 20 kilotons. For the observed energy higher than 1MeV, the rate for possible backgrounds has not been given in Ref. [31], but the signal rate is highly suppressed, rendering the realistic observation to be difficult.

D. 12 C channel
There are three different reaction channels for the interaction between astrophysical neutrinos with the 12 C nuclei in the LS detectors. The first one is the NC reaction ν + 12 C → ν + 12 C * , where the final-state nucleus resides in the excited state (15.1 MeV, 1 + ) and its deexcitation leads to a gamma ray of 15.1 MeV that can be registered in the detector. The NC event spectrum is given by where N C = 3N e /23 is the total number of 12 C nuclei in the LS target [43] and σ NC ν α denotes the total cross section. The other two are charged-current (CC) reactions, namely, ν e + 12 C → e − + 12 N and ν e + 12 C → e + + 12 B. The CC event spectrum is found to be where T e is the recoil energy of the final-state electron (or positron) and σ CC ν e and σ CC ν e denotes respectively the cross section of ν e -12 C and ν e -12 C reaction. Moreover, we have T e = E ν e − 16.827 MeV for the 12 C(ν e , e − ) 12 N channel while T e = E ν e −13.880 MeV for the 12 C(ν e , e + ) 12 B channel [44]. The cross sections of all three neutrino-carbon reaction channels can be found in Table 1 of Ref. [44].
Similar to the pES channel, only high-energy neutrinos and antineutrinos can produce the signals with energies observable in the LS detectors. Therefore, we consider only the PBHs with the lightest mass, i.e., M PBH = 10 15 g. The event spectra of NC and CC reactions are shown in Fig. 5 and Fig. 6, respectively, where the effective operation time of 10 years and the fiducial LS mass of 20 kilotons have been taken as before. In comparison with the inverse beta decay (IBD) channel, the event rates in the 12 C channel are much smaller and will not contribute much to the detection of PBHs.

IV. CONSTRAINTS FROM JUNO
In the previous discussions, we have demonstrated that the IBD channel is most sensitive to the antineutrino signals from the PBHs as dark matter. For this reason, only the IBD signal will be implemented to derive the constraints on the PBHs at JUNO in this section. By requiring the signal-to-noise ratio SNR ≥ 1.6/ √ N (i.e. 90% C.L.), where N is the total event number of backgrounds for 10 years, we have shown the upper limit on f PBH for a given M PBH as the red curve in Fig. 7. For comparison, the currently best limit on f PBH from Super-K has been plotted the blue curve and the shaded area in blue has been excluded [19]. Two important observations can be made. First, after running for 10 years, JUNO will have the capability to explore the pink area in Fig. 7, which has not been constrained by current neutrino observatories. To be specific, given M PBH = 10 15 g, the upper limit f PBH 3 × 10 −5 can be obtained from JUNO, which is about twenty times better than that f PBH 6 × 10 −4 from Super-K. Second, JUNO will be able to constrain the PBH dark matter in the mass range M PBH = (5 ∼ 8) × 10 15 g, for which f PBH = 1 is still allowed by Super-K. As we have mentioned before, although there exist other observational limits on f PBH from cosmic gamma rays, the neutrino observations will not only provide an independent limit on f PBH but also a novel way to probe the mechanism of Hawking radiation of PBHs.

V. CONCLUDING REMARKS
Motivated by the attractive scenario of PBH dark matter, we have calculated the PBH-induced neutrino and antineutrino fluxes due to the Hawking radiation and the corresponding neutrino event spectra at JUNO. The main conclusion is that the current limit on the PBH abundance f PBH from Super-K will be improved by one order of magnitude at JUNO.
One should notice that PBHs have been assumed to follow the monochromatic mass distribution and to be spinless. However, there are indeed different theoretical scenarios, in which the predicted mass distribution of PBHs is broad [45][46][47][48] or the PBH spins are high [49,50]. For the PBHs with a log-normal mass distribution, the shaded regions in Fig. 7 would be shallower but broader (e.g., see Fig. 2 in Ref. [19]). For the spinning PBHs, the constraints would be more stringent because of more intense Hawking radiation [26,27]. Therefore, the exclusion limits obtained in the present work can be regarded to be conservative. Our approach can be easily generalized to study more generic scenarios of PBHs, which is however beyond the scope of this paper. As shown in Eq. (4), a different choice of the density profile of dark halo may affect the neutrino fluxes from the Galaxy. Besides the NFW profile, we have also considered two other typical scenarios, i.e., the cored isothermal (ISO) profile and the Einasto (EIN) profile [36]. Define the relative error as err(X) ≡ |dF Gal /dE(X) − dF Gal /dE(NFW)| dF Gal /dE(NFW) .
(13) We find that its magnitude is about 1.6%, 0.3% for X = ISO and EIN, respectively. These errors are small enough for us to safely ignore the uncertainty from different choices of the density profiles. Therefore, we just take the NFW profile in the present work for illustration.
Another concern is related to the pES channel for the detection of neutrinos and antineutrinos from PBHs. Compared to the IBD channel, a considerable event rate produced from PBHs is also expected in this channel. However, we have shown that it is severely contaminated by the backgrounds, which have not been well studied so far. The derived limit from JUNO in Fig. 7 will be further improved if new analysis techniques are developed to effectively get rid of such backgrounds [31]. This exploration will be left for future works.