# Primordial black hole probes of heavy neutral leptons

Valentina De Romeri  0000-0003-3585-7437<EMAIL_ADDRESS>Instituto de
Física Corpuscular (CSIC-Universitat de València), Parc Científic UV C/
Catedrático José Beltrán, 2 E-46980 Paterna (Valencia) - Spain  Yuber F.
Perez-Gonzalez  0000-0002-2020-7223<EMAIL_ADDRESS>Institute for Particle Physics Phenomenology, Durham University, South Road,
Durham DH1 3LE United Kingdom.  Agnese Tolino  0009-0003-3278-0902
<EMAIL_ADDRESS>Instituto de Física Corpuscular (CSIC-Universitat de
València), Parc Científic UV C/ Catedrático José Beltrán, 2 E-46980 Paterna
(Valencia) - Spain

###### Abstract

Primordial black holes (PBH), while still constituting a viable dark matter
component, are expected to evaporate through Hawking radiation. Assuming the
semi-classical approximation holds up to near the Planck scale, PBHs are
expected to evaporate by the present time, emitting a significant flux of
particles in their final moments, if produced in the early Universe with an
initial mass of $\sim 10^{15}$ g. These “exploding” black holes will release a
burst of Standard Model particles alongside any additional degrees of freedom,
should they exist. We explore the possibility that heavy neutral leptons
(HNL), mixing with active neutrinos, are emitted in the final evaporation
stages. We calculate the expected number of active neutrinos from such an
event, including contributions due to the HNL decay for different assumptions
on the mixings. We infer sensitivities on the active-sterile neutrino mixing
and on the sterile neutrino mass, finding that, for instance, for the scenario
where $U_{\tau 4}\neq 0$, IceCube could improve current constraints by $\sim
2$ orders of magnitude, for HNLs masses between 0.1 - 1 GeV, for a PBH at a
distance of $\sim 10^{-4}$ pc from Earth.

primordial black holes, neutrinos, heavy neutral leptons, IceCube

††preprint: IPPP/24/23

##  I Introduction

One major shortcoming of the Standard Model (SM) of particle physics is the
lack of an explanation of the largest matter component of the Universe.
Numerous observations at different cosmological scales indicate that most
matter in the Universe is made up of dark matter (DM), whose nature and
properties are still unknown [1].

While most observations seem to prefer a particle nature for DM, thus calling
for the introduction of new particles beyond the Standard Model (BSM), null
searches at typical DM direct detection experiments [2, 3] have recently
revived the interest in alternative candidates. One such an option is that the
DM is made not by BSM particles but by Primordial Black Holes (PBHs), see e.g.
Refs. [4, 5, 6] for recent reviews on PBH-DM.

PBHs are hypothetical black holes that would have formed in the early Universe
before the Big-Bang Nucleosynthesis (BBN) epoch [7, 8, 9, 10]. Black holes
have received renewed interest in recent years after the observation of
gravitational waves in the LIGO-Virgo experiments [11, 12], which have
definitely proven their existence and opened the possibility of searching for
black holes of primordial origin. There are different constraints on PBHs
depending on their masses and whether they still exist at present. Crucially,
some of the strongest bounds rely on the particle production from the black
hole gravitational field, i. e. via Hawking radiation [13, 14]. As
demonstrated by Hawking, black holes have a temperature inversely proportional
to their masses and they are expected to evaporate by emitting a flux of
particles with thermal spectra. Clearly, the observation of PBH evaporation
would confirm this prediction and provide invaluable information about several
aspects of physics. However, non-detection of PBH evaporation allows to impose
a large number of constraints on the abundance of PBHs with masses in the
range $10^{15}\sim 10^{17}$ g, using $\gamma$ rays [15, 16, 17, 18, 19, 20],
electrons/positrons [21, 22, 23, 24] and neutrinos [22, 25, 26, 27]. Lighter
PBHs, that would have evaporated at earlier times in the Universe’s history,
can be constrained through BBN [28, 29, 30], reionization [31, 32] or CMB [28,
29] measurements. See, e.g., Ref. [33] for a comprehensive review.

Non-rotating PBHs with masses $\lesssim 5\times 10^{14}$ g would have
completely evaporated by now [34, 35, 36] and could not constitute the
totality of DM. However, they could still give rise to interesting signals.
Following an initial stage during which the evaporation proceeds slowly, as
the PBH mass decreases, its temperature and flux increase, leading to a final
stage that would resemble an explosion [13]. The luminosity of such exploding
PBHs is established by the types of degrees of freedom that they can emit,
together with their instantaneous properties (mass spectrum, charge and spin).
Clearly, photons would be one smoking-gun signature of these explosions,
having strong matches with a gamma-ray burst though with some intrinsic
differences [37]. Facilities such as H.E.S.S [38, 39], Milagro [40], VERITAS
[41], Fermi-LAT [42], and HAWC [43, 44] have already placed constraints on the
rate of PBH explosions. Notice that, in order to be detectable by current (or
even near-future) telescopes, they have to occur (at most) a few parsecs from
Earth [45, 37, 46]. Currently, the H.E.S.S array of imaging atmospheric
Cherenkov telescopes [47] sets the strongest direct limit on the rate of
exploding PBHs:
$\dot{n}_{\mathrm{PBH}}<2000~{}\mathrm{pc}^{-3}~{}\mathrm{yr}^{-1}$ for a
burst interval of 120 seconds, at the $95\%$ confidence level.

Evaporating PBHs will radiate all elementary particles with a mass below their
temperature. Among the SM states, neutrinos constitute one possible signature
of exploding PBHs [48, 49, 50, 46, 51]. On the other hand, new degrees of
freedom that may appear in BSM scenarios could also be probed through the
final stages of PBHs evaporation [52, 53, 54, 55, 56, 33, 57]. For instance,
several BSM extensions call upon the introduction of SM-singlet degrees of
freedom, dubbed sterile neutrinos or heavy neutral leptons (HNLs), to address
neutrino masses and mixings, another shortcoming of the SM and actually the
first laboratory probe of BSM physics. HNLs can have masses that vary over
several orders of magnitude, depending on the actual mechanism that originates
neutrino masses, and extended phenomenological implications. Despite intense
experimental efforts in searching for these sterile fermions, no positive
evidence has been found to date [58, 59, 60]. However, HNLs in the MeV-GeV
mass range remain appealing BSM candidates, as they emerge rather naturally in
several SM extensions (e.g. [60]) besides being accessible also at terrestrial
facilities. Interestingly, the existence of HNLs and PBHs in the early
Universe could affect leptogenesis scenarios, as noted in Refs. [61, 62, 63,
64, 65, 66, 67, 68, 69, 70].

In this paper we consider the constraints that could be derived on HNLs after
an exploding PBH in close proximity to Earth is observed in a neutrino
telescope. We estimate sensitivities at IceCube [71, 72] as a benchmark
facility to observe high-energy neutrinos from a PBH burst. IceCube is an ice
Cherenkov detector located in the South Pole with a size of about 1 ${\rm
km^{3}}$. It has an impressive angular resolution of $\lesssim 1\degree$ for
high-energy muon-tracks with energies $E\sim{\cal O}({\rm TeV})$. Thus, we
expect IceCube to have the capability to determine the PBH position and
discriminate neutrino events from possible backgrounds [49]. We will assess
potential modifications to the neutrino spectra detected by IceCube if HNLs
are generated during the final stages of PBH evaporation and subsequently
decay into active neutrinos. While unresolved issues persist regarding black
hole evaporation physics, such as the information paradox [73, 74, 75] or the
thermal nature of the Hawking spectrum post-Page time [76, 77], and how
alterations would impact the mentioned constraints — see e. g. [78, 79, 80]
for an example of possible modifications in a specific prototype model — we
adopt a more neutral stance ahead, and assume the semi-classical
approximation’s validity until near the Planck scale.

The remainder of this work is organized as follows. We discuss the primary and
secondary SM contributions to the muon neutrino spectra expected at IceCube in
Section II. Section III is instead devoted to summarize the relevant secondary
contributions from HNL decays. We present in Sec. IV the details of the
statistical analysis performed to estimate the sensitivity at IceCube. We show
our results in Sec. V, in terms of sensitivities on the HNL parameter space.
Finally we conclude in Sec. VI. Throughout this work, we consider natural
units where $\hbar=c=k_{\rm B}=1$, and define the Planck mass to be
$M_{p}=1/\sqrt{G}$, with $G$ the gravitational constant.

##  II Standard Model contributions

Among different possibilities, large density fluctuations in the very early
Universe may have led to the creation of a PBH population that could persist
until the present day. Typically, the initial mass of PBHs is linked to the
particle horizon mass when they formed. Given the lack of clear constraints on
the Universe’s evolution before BBN, PBHs could have masses ranging from
around ${\cal O}(1)$ g to several orders of magnitude the Solar mass [28, 29].
Within this broad range of masses, a PBH could have formed with precisely the
right mass for its lifetime to match the age of the Universe, thus entering
the final explosive stages of its evolution at present. In the standard
evaporation scenario, considering only SM degrees of freedom emitted during
evaporation, a PBH with a mass of approximately $\sim 10^{15}$ g would be
evaporating today.

Therefore, by observing the absence of gamma-ray bursts in our Galaxy, we can
potentially limit the quantity of evaporating PBHs at present. If we assume a
constant universal number density, this translates into an upper limit on the
PBH number density $n_{\rm PBH}$ equal to [46]

$\displaystyle n_{\rm PBH}\lesssim
0.35\left(\frac{\beta^{\prime}}{10^{-29}}\right)\left(\frac{10^{15}~{}{\rm
g}}{M_{\rm in}}\right){\rm pc^{-3}},$ (1)

with $M_{\rm in}$ the initial PBH mass and $\beta^{\prime}$ the initial PBH
energy density at formation, normalized to the total energy density. Note that
gravitational clustering effects could relax this limit. This implies that in
a spherical volume with a radius of a parsec, we would anticipate
approximately $\sim 1.5$ PBHs evaporating presently to align with
observations. Therefore, it is conceivable that one of these PBHs may be in
close proximity to Earth, nearing the end of its lifespan, presenting an
opportunity for observation. This is the scenario we will presume henceforth.

Assuming a semi-classical picture, where matter fields, treated as quantum,
propagate in a classical background describing the gravitational field of a
black hole, Hawking determined the particle emission after the gravitational
collapse [13, 14]. The particle emission rate depends on the black hole
instantaneous properties, such as its mass $M$, angular momentum $J$ and
charge $Q$. Nevertheless, the charge and angular momentum111However, note that
if there exist a dark sector containing scalar degrees-of-freedom in large
quantities $N>100$, the angular momentum evaporates in a significantly
different rate, see e. g. Refs. [55, 56, 57, 51]. are expected to be depleted
in a faster rate than the mass. Thus, we can safely assume that PBH reaching
the final stages of their lifetime are Schwarzschild, described uniquely by
their masses $M$.

If we accept the semi-classical approach as valid until close to the Planck
scale, the PBH mass loss rate can be expressed as

$\displaystyle\frac{dM}{dt}$
$\displaystyle=-\varepsilon(M)\frac{M_{p}^{4}}{M^{2}}\,,$ (2)

where $\varepsilon(M)$ is the evaporation function that describes the degrees-
of-freedom (dof) that could be emitted for a given PBH mass [81, 82, 83, 84]

$\displaystyle\varepsilon(M)$
$\displaystyle=\sum_{i=\text{dofs}}\frac{g_{i}}{128\pi^{3}}\int_{0}^{\infty}\frac{x\,\Gamma_{s_{i}}(x)}{\exp(x)-(-1)^{2s_{i}}}\,dx\,.$
(3)

Here $g_{i}$ are the internal dofs associated to the particle species $i$
having spin $s_{i}$, and $\Gamma_{s_{i}}$ indicates the absorption
probability, also known as greybody factors, describing potential
backscattering arising from centrifugal and/or gravitational potentials
related to the curved spacetime around the black hole. Such quantities are
computed using a numerical approach as the one described in Refs. [34, 35,
85]. For convenience, we rely on the absorption probabilities from the code
BlackHawk [18, 86], which are valid for massless particles. We have verified
that using the correct probabilities for massive fermions from Refs. [87, 88]
yields similar results to those obtained from the tables provided in
BlackHawk. The variable $x$ is defined as the ratio $x=E/T$, where $E$ is the
particle’s energy and $T$ the Hawking temperature, related to the PBH mass via

$\displaystyle T=\frac{1}{8\pi GM}\sim 1~{}{\rm TeV}\left(\frac{10^{10}~{}{\rm
g}}{M}\right).$ (4)

From this we can infer that dofs having masses $m\lesssim 1~{}{\rm TeV}$ will
be emitted by a black hole with $M=10^{10}~{}{\rm g}$ without any Boltzmann
suppression, even those that do not interact with the SM. If the number of
dofs is constant during the whole PBH lifetime, the mass lose rate is exactly
solvable,

$\displaystyle M(t)=M_{\rm in}\left(1-\frac{t}{\tau_{d}}\right)^{1/3},$ (5)

with $\tau_{d}$ representing the PBH lifetime. The lifetime is given by

$\displaystyle\tau_{d}=\frac{M_{\rm in}^{3}}{3\varepsilon M_{p}^{4}}\sim
428~{}{\rm s}\left(\frac{4.07\times
10^{-3}}{\varepsilon}\right)\left(\frac{M_{\rm in}}{10^{10}~{}{\rm
g}}\right)^{3},$ (6)

where we used $\varepsilon=4.07\times 10^{-3}$ corresponding to the value of
the evaporation function for the SM dofs.

The instantaneous emission rate for a fermion (boson) $i$ from PBH evaporation
is given by a Fermi-Dirac (Bose-Einstein) thermal spectrum corrected by the
absorption cross section mentioned above,

$\displaystyle\frac{d^{2}N_{i}}{dEdt}(M)=\frac{g_{i}}{2\pi}\frac{\Gamma_{s_{i}}(E)}{\exp(8\pi
GME)-(-1)^{2s_{i}}}.$ (7)

In experimental terms, optical and neutrino telescopes are anticipated to
detect the final burst from PBHs within a finite time frame, measuring a burst
of photons or neutrinos over a specific duration. Therefore, integrating the
instantaneous spectrum over the assumed time interval is essential to
calculate the fluence for a given experiment. Since any particles emitted
prior would evade detection, we set the initial mass of an exploding PBH to
yield a lifetime matching the observation period. For instance, if the
observation time is 100 s, we adjust $M_{\rm in}=6.2\times 10^{9}$ g. The
time-integrated spectrum in a time interval $\tau$ is therefore given by

$\displaystyle\frac{dN_{i}}{dE}=\int_{0}^{\tau}dt\,\frac{d^{2}N_{i}}{dEdt}(M(t))\,.$
(8)

Once such time-integrated spectrum is determined for a given particle physics
model, it will be possible to estimate the observed events in a neutrino
telescope. However, it is important to note that observable neutrinos,
specifically active neutrinos, will be produced in two different forms: first,
directly from evaporation, and secondly from the decay of unstable SM
particles and HNL of interest. Next, we provide a detailed account of these
emissions.

###  II.1 Primary emission

Direct neutrino emission from evaporation has been considered since the
discovery of such a phenomenon. Moreover, neutrinos were considered as a
benchmark for analyzing the emission of massless fermions from a black hole,
due to the belief, at the time, that neutrinos were massless. Nevertheless,
following the discovery of neutrino oscillations and consequently neutrino
masses, the primary emission model has been revised in Ref. [89]. The main
consequence of being massive particles for the present case study is that
neutrinos are assumed to be emitted as mass eigenstates $\nu_{1,2,3}$, instead
of the flavour eigenstates $\nu_{e,\mu,\tau}$ considered in early studies.

To detect the primary neutrinos from evaporation, we would still need to
utilize weak interactions. The time-integrated spectrum of an observable state
with flavour $\alpha$ is then obtained by projecting the spectra of the mass
eigenstates into the flavour basis by employing the Pontecorvo-Maki-Nakagawa-
Sakata (PMNS) mixing matrix $U_{\alpha i}$,

$\displaystyle\left.\frac{dN_{\nu_{\alpha}}}{dE}\right|_{\rm
pri}=\sum_{i=1}^{3}|U_{\alpha i}|^{2}\frac{dN_{\nu_{i}}}{dE},$ (9)

with $dN_{\nu_{i}}/dE$ the time-integrated spectrum for the massive state
$\nu_{i}$ given by Eq. (8) for the fermionic case. The primary emission
reflects the thermal spectrum, corrected by the greybody factors, expected
from the evaporation. Therefore, such emission will be peaked close to a value
of $\sim 6T_{\rm PBH}^{\rm in}$, $T_{\rm PBH}^{\rm in}$ being the initial PBH
temperature. For the exploding PBHs that we are interested in, we find that
such a primary spectrum is peaked at large energies, $E\sim{\cal O}(10~{}{\rm
TeV})$, but, as we will see shortly, is subdominant to secondary emission.

###  II.2 Secondary emission

The largest part of the overall time-integrated spectrum of neutrinos does not
come from the direct black hole emission but from the weak decays of unstable
particles produced during the evaporation [81, 82]. If we only consider the SM
particle content for the moment, quark hadronization and hadron posterior
decays, together with heavy leptons and massive gauge bosons decays, will
produce a large amount of neutrinos since these dofs dominate the SM spectrum.
To estimate the time-integrated spectrum of these neutrinos, known as the
secondary component, we need to use tools that compute the spectra from
hadronization plus decay considering SM gauge interactions, that are
convoluted with the Hawking spectrum of parent particles $j$,

$\displaystyle\left.\frac{dN^{\nu_{\beta}}_{\rm PBH}}{dE}\right|_{\rm
sec}=\sum_{j}\int_{0}^{\infty}dE_{j}\frac{dN_{j}}{dE_{j}}\frac{dN(j\longrightarrow\nu_{\beta})}{dE}(E,E_{j}),$
(10)

where $dN(j\longrightarrow\nu_{\beta})/dE$ are the neutrino spectra from the
hadronization/decay of the parent particle $j$ and $dN_{j}/dE_{j}$ is its
primary time-integrated spectrum. Finally, the subscript “PBH” indicates that
the time-integrated spectrum is computed close to the PBH. As an example, we
have that quark hadronization will produce large amounts of $\pi^{\pm}$ that
subsequently decay generating neutrinos in their flavour eigenstates,
$\pi^{\pm}\to\mu^{\pm}+\accentset{(-)}{\nu}_{\mu}\to
e^{\pm}+\accentset{(-)}{\nu}_{e}+\accentset{(-)}{\nu}_{\mu}$, thus producing
three neutrinos per each pion produced. Therefore, we expect the secondary
neutrino spectrum to dominate over the primary, specially at low energies. We
use BlackHawk222We use BlackHawk v2.2 for the SM contributions and BlackHawk
v2.1 for the BSM primary contribution, since v2.2 lead to some issues in the
time-dependent spectrum when a HNL is added. [18, 86] to evaluate both the
primary and secondary spectra of neutrinos. In particular, we rely on
HDMSpectra [90], already integrated in BlackHawk, to compute the secondary
emission through the hadronization of primary particles, whose range of
validity extends up to the Planck scale.

As such secondary neutrinos are emitted as flavour eigenstates, we need to
consider neutrino oscillations to compute the correct spectrum. Given the
distances involved here, significantly larger than oscillation distances
$\sim{\cal O}({\rm km})$, neutrino oscillations undergo decoherence in their
route from the PBH to the Earth. The secondary contribution far from the PBH,
including decohered oscillations, is given by

$\left.\frac{dN_{\nu_{\alpha}}}{dE}\right|_{\rm
sec}=\sum_{\beta=e,\mu,\tau}\sum_{i=1,2,3}|U_{\alpha i}|^{2}|U_{\beta
i}|^{2}\left.\frac{dN^{\nu_{\beta}}_{\rm PBH}}{dE}\right|_{\rm sec}.$ (11)

Thus, events from a PBH burst would consist in the combined contribution of
both primary and secondary contributions defined above. However, if additional
BSM dofs exist which decay producing additional neutrinos, the prediction of
the secondary contribution would be significantly modified, and therefore the
observation of a PBH burst could provide additional information to constrain
the presence of such BSM states. We will consider next a simple, but well
motivated, scenario for BSM physics, and the implications in a PBH burst.

##  III Heavy Neutral Lepton contributions

###  III.1 Heavy Neutral Lepton Generalities

The experimental observation of neutrino oscillations has firmly established
the existence of neutrino masses, indicating the existence of BSM physics. To
accommodate neutrino masses, a minimal possibility consists in including $n$
right-handed fermions $N_{i}$, $i=\\{1,2,\ldots,n\\}$, singlets under the SM
gauge group. A simple model that can realize this setup is the type-I _seesaw_
mechanism [91, 92, 93, 94, 95]. Crucially, the SM symmetries allow for masses
$M^{ij}_{R}$ for such right-handed neutrinos (RHN) of the Majorana type. The
most general mass lagrangian is given by

$\displaystyle\mathscr{L}_{\rm RHN}^{m}=-Y_{\alpha
i}\overline{L}_{\alpha}\widetilde{H}N_{i}-\frac{1}{2}M^{ij}_{R}\overline{N_{i}^{c}}N_{j}+{\rm
h.c.},$ (12)

with $\overline{L}_{\alpha}$ the lepton doublets, $\widetilde{H}$ the
conjugate Higgs doublet, and $Y_{\alpha i}$ Yukawa parameters related to the
couplings between the aforementioned doublets and the RHNs. To address the
observed neutrino oscillation pattern [96, 97, 98], where two non-zero
quadratic mass differences have been measured, the minimal number of RHNs is
$n=2$. However, this does not forbid the existence of additional singlets.

After electroweak symmetry breaking, the lagrangian above can be rewritten as

$\displaystyle\mathscr{L}_{\rm RHN}^{m}=-\frac{1}{2}\overline{{\cal
N}_{L}^{c}}M_{\nu}{\cal N}_{L}+h.c.,$ (13)

where

$\displaystyle{\cal N}_{L}=\begin{pmatrix}\nu_{L}\\\
(N_{R})^{c}\end{pmatrix},\quad
M_{\nu}=\begin{pmatrix}\mathbf{0}_{3}&Yv/\sqrt{2}\\\
Y^{T}v/\sqrt{2}&M_{R}\end{pmatrix}\,,$ (14)

being $\nu_{L}=(\nu_{e},\nu_{\mu},\nu_{\tau})^{T}$ and
$N_{R}=\\{N_{1},\ldots,N_{n}\\}$ active and sterile neutrino states,
respectively, and $v$ the Higgs vacuum expectation value. After diagonalizing
the Lagrangian in Eq. (13) to form $M_{\nu}={\mathcal{U}}M_{\nu}^{\rm
diag}{\mathcal{U}}^{T}$, utilizing a unitary rotation matrix ${\mathcal{U}}$
with $(3+n)\times(3+n)$ entries, the resulting spectrum would encompass $3+n$
eigenstates denoted as ${\cal N}_{L}^{m}=\\{\nu_{1},\ldots,\nu_{3+n}\\}$.
Among these, 3 are considered “active”, representing weakly interacting
states, while the remaining $n$ are predominantly “sterile”, as their
interactions with other SM dofs occur solely through mixing parametrized via
${\mathcal{U}}$. The modified PMNS matrix now corresponds to a $3\times 3$
block of ${\mathcal{U}}$ and can show deviations from unitarity.

If we assume a large hierarchy between the RHN mass and electroweak scales,
i.e., $M_{R}\gg vY$, Majorana masses for the active neutrinos
$\\{\nu_{1},\nu_{2},\nu_{3}\\}$ can be explicitly written as333Note that while
being physical, related to CP-properties of active neutrinos, the negative
sign can be reabsorbed through a redefinition of fields.,

$\displaystyle m_{\nu}\approx-\frac{v^{2}}{2}YM_{R}^{-1}Y^{T}.$ (15)

In the scenario where the Yukawa matrix has entries of $|Y_{\alpha
i}|\sim{\cal O}(1)$, we would require $M_{R}\sim{\cal O}(10^{12})~{}{\rm
GeV}$, that is, the RHN masses are close to the grand-unification scale, to
have active neutrino masses $m_{\nu}\sim{\cal O}({\rm eV})$, consistent with
oscillation experiments. This case, corresponding to the well-known seesaw
mechanism, would then explain the observed smallness of neutrino masses. If
such large hierarchy does exist, then the mixing between active and RHN would
be extremely suppressed, ${\mathcal{U}}\approx\frac{v}{\sqrt{2}}YM_{R}^{-1}$.
Therefore, the RHN would have tiny mixings with the active sector, making a
direct experimental observation of such dofs quite challenging.

However, since the SM does not constrain either the scale of Majorana masses
or the values of the Yukawa parameters $Y_{\alpha i}$ in Eq. (12), it is
possible for $M_{R}$ to be at the electroweak scale and still explain neutrino
masses. For instance, this occurs in low-energy realizations of the seesaw
scenario, like the inverse [99] and the linear seesaw [100, 101, 102]. In such
a case, the seesaw formulae for the masses and mixing are not valid and need
to be extended. Thus, in general, the diagonalization of Eq. (13) would imply
the presence of light sterile neutrinos, which could exhibit substantial
mixing with the active states. This opens up the possibility of experimental
detection of these light sterile states, i. e., HNLs.

Various experimental searches of these states have been conducted and are
planned for future facilities, as comprehensively reviewed in Ref. [60]. In
this context, our focus lies on understanding the impact of such additional
states on the observed neutrino events arising from a PBH burst. To simplify
matters, we assume the presence of a single HNL state $\nu_{4}$ with a mass
$m_{4}$, which interacts with active neutrinos through the mixing $U_{\alpha
4}$, where $\alpha={e,\mu,\tau}$. In such a simplified phenomenological
scenario

$\displaystyle{\cal U}=\begin{pmatrix}U_{\alpha i}&U_{\alpha 4}\\\
U_{4i}&U_{44}\end{pmatrix}\,,$ (16)

while the flavor states $\nu_{\alpha}$ can be written as a superposition of
the mass eigenstates as $\nu_{\alpha}=\sum_{i=1}^{3}U_{\alpha
i}\nu_{i}+U_{\alpha 4}\nu_{4}$.

Before we explore the effects of an additional HNL on a PBH burst, it is
important to address the unitarity of the neutrino mixing matrix used to
calculate the active neutrino events in Sec. II. In a standard scenario, where
no sterile states are included, the mixing angles are fixed according to [96],
and the PMNS is a priori unitary. However, in a non-standard scenario, where
sterile neutrinos are introduced in the picture and mix with at least one
neutrino species, the unitarity of the $3\times 3$ block of ${\mathcal{U}}$ is
not guaranteed, and hence

$\displaystyle\sum_{i=1}^{3}|U_{\alpha i}|^{2}\neq 1\,,$ (17)

in general. Nevertheless, we have verified that taking the $3\times 3$
“active” block of ${\mathcal{U}}$ as unitary is a good approximation as long
as the mixing between active neutrinos and the additional state is smaller
than $|U_{\alpha 4}|^{2}\lesssim 10^{-3}$.

Should HNLs exist and be produced during PBH evaporation, they will impact its
mass evolution and give rise to additional secondary contributions. The rate
at which PBHs lose mass is primarily influenced by the number of emitted
degrees of freedom [53], with only minor dependence on the properties of the
emitted species, whether fermionic or bosonic, and their actual spin [34]. As
a consequence, if the SM particle content is extended only by one type of HNL,
the PBH mass evolution will not be significantly modified. On the other hand,
the possible emission of HNLs in multiple generations would increase the
evaporation rate accelerating the PBH explosion. We have checked that a large
($\gtrsim 10$) number of generations is required in order to observe sizeable
effects in the primary emission, in agreement with [53].

Since we will work under the assumption that only one generation of HNLs is
present, the primary neutrino spectra emitted in the PBH explosion are the
same as in the SM. However, the HNLs produced during the PBH evaporation will
generate an additional secondary flux of active neutrinos from their decays.
The morphology of such a secondary contribution, expected to induce visible
signals at IceCube, will depend on the HNL mass regime, on the available decay
channels, and which active neutrino flavour the sterile neutrino mixes with.

At the scope of our analysis we assume an observation time of 100 s in a
neutrino telescope like IceCube, corresponding to the explosion of a PBH with
mass $M\sim 6\times 10^{9}$ g. This will constrain the maximum HNL mass that
can be produced in the evaporation process to $m_{4}\lesssim 2$ TeV, a value
that is dictated by the temperature of the PBH, as above that threshold the
HNL emission is heavily suppressed. HNLs lighter than this value can be
produced and give rise to an observable flux of daughter active neutrinos. For
the sake of illustration, in the following we will consider two separate HNL
regimes: a light-mass ([0.1-1] GeV) and a heavy-mass ([0.5-2] TeV) regime.

###  III.2 Light-mass regime: [0.1-1] GeV

HNL decay modes have been studied extensively in the literature, see e.g.
Refs. [103, 58, 104, 105, 106, 107, 108, 109, 60, 110]. HNLs with masses
$\mathcal{O}(100)$ MeV mainly decay into final states that include one lepton
and one meson or three leptons. Below $1$ GeV, the following two channels with
active neutrinos in the final state dominate over the others

$\displaystyle\nu_{4}$ $\displaystyle\to$
$\displaystyle\nu_{\alpha}\pi^{0}\,,$ (18) $\displaystyle\nu_{4}$
$\displaystyle\to$
$\displaystyle\nu_{\alpha}\nu_{\ell}\bar{\nu}_{\ell}\,~{}~{}(\ell=e,\mu,\tau)\,,$

where $\alpha$ indicates the neutrino flavor accessible depending on which
HNL-active neutrino mixing $U_{\alpha 4}$ is activated. In what follows, we
assume neutrinos to be Majorana, so that the charge-conjugated process will be
also taken into account for each HNL decay process that includes neutrinos in
the final states. Moreover, we will consider three different scenarios for the
mixings: a first case where only the mixing with the electron neutrino is
activated, $|U_{e4}|^{2}\neq 0,|U_{\mu 4}|^{2}=0,|U_{\tau 4}|^{2}=0$ (dubbed
1:0:0); a second case where only mixing with muon neutrinos is allowed,
$|U_{e4}|^{2}=0,|U_{\mu 4}|^{2}\neq 0,|U_{\tau 4}|^{2}=0$ (case 0:1:0); and
case 0:0:1, $|U_{e4}|^{2}=0,|U_{\mu 4}|^{2}=0,|U_{\tau 4}|^{2}\neq 0$. While
in the case 0:1:0 the muon neutrino contribution is directly produced, in the
other two cases $\nu_{\mu}$ appears due to oscillations of the other flavor
components during their propagation from the PBH to the Earth.

Given these considerations, we estimate the time-integrated secondary spectrum
of neutrinos with flavor $\alpha$, produced through the decay of HNLs, as
[111, 108, 112, 110]

$\frac{dN_{\alpha}}{dE}=\mathcal{B}_{a}\int d\cos{\theta}\int_{E_{s,\rm
min}}^{E_{s,\rm
max}}dE_{s}\frac{1}{\gamma_{s}\left(1+\beta_{s}\cos{\theta}\right)}\frac{dN_{s}}{dE_{s}}\mathcal{F}_{\alpha}\left[\frac{E}{\gamma_{s}\left(1+\beta_{s}\cos{\theta}\right)},\cos{\theta}\right]\,,$
(19)

where $E$, $E_{s}$ are the energies of $\nu_{\alpha}$ and the HNL in the
laboratory frame, respectively, and $\theta$ is the angle formed between the
direction of $\nu_{\alpha}$ in the HNL rest frame and its velocity in the
laboratory frame [108]. Moreover, in Eq. (19), $\gamma_{s}=E_{s}/m_{4}$ and
$\beta_{s}=p_{s}/E_{s}$, $m_{4}$ are the HNL mass and $p_{s}$ its momentum,
while $\mathcal{B}_{\alpha}$ = $\Gamma_{\alpha}/\Gamma_{\mathrm{tot}}$
indicates the branching ratio of the decay process under consideration, with
$\Gamma_{\mathrm{tot}}$ and $\Gamma_{\alpha}$ being respectively the total and
partial decay widths. $\frac{dN_{s}}{dE_{s}}$, instead, denotes the total
primary spectrum of HNLs directly emitted by the PBH through Hawking
evaporation.

The integration limits $E_{s,\rm min}$ and $E_{s,\rm max}$ in Eq. (19) refer
to the minimal and maximal energy that the HNL can have in the laboratory
frame, respectively. In principle, these values can vary in the range
$[10^{-5}-10^{5}]$ $T^{\rm in}_{\mathrm{PBH}}$ [18, 86], as for a
Schwarzschild black hole the Hawking radiation is mainly peaked around its
temperature. However, they will eventually depend on the specific features of
the decay channel under consideration. In particular, $E_{s,\rm min}$ is the
minimal energy the HNL should reach in order to produce an active neutrino
$\nu_{\alpha}$ with energy $E$, whereas $E_{s,\rm max}$ is the maximal energy
value above which $\frac{dN_{s}}{dE_{s}}$ becomes negligible. Finally,
$\mathcal{F}_{\alpha}$ in Eq. (19) represents the angular and energetic
distribution $\mathcal{F}_{\alpha}\equiv dN_{\alpha}/dE^{\prime}d\rm
cos\theta$ that the active neutrino has in the HNL rest frame, and depends on
the specific decay channel under consideration [108, 110, 112].

Note that the energies of the active neutrino in the laboratory frame $E$ and
in the HNL rest frame $E^{\prime}$, the sterile neutrino energy $E_{s}$ and
$\theta$ are not independent quantities, but they are all related through the
following boost relation

$E=\gamma_{s}E^{\prime}\left(1+\beta_{s}\cos{\theta}\right).$ (20)

Therefore, Eq. (19) can be conveniently rewritten in terms of the energy
$E^{\prime}$ that the active neutrino $\nu_{\alpha}$ has in the HNL rest frame
by simply performing the following change of variables [108]

$d\cos{\theta}=-\frac{1}{\beta_{s}\gamma_{s}}\frac{E}{E^{{}^{\prime}2}}dE^{\prime}\,.$
(21)

By substituting Eq. (21) into Eq. (19) we obtain

$\frac{dN_{\alpha}}{dE}=\mathcal{B}_{\alpha}\int_{E^{\prime}_{\mathrm{min}}}^{E^{\prime}_{\mathrm{max}}}dE^{\prime}\int_{E_{s,\rm
min}}^{E_{s,\rm
max}}dE_{s}\frac{1}{\beta_{s}\gamma_{s}}\frac{1}{E^{\prime}}\frac{dN_{s}}{dE_{s}}\mathcal{F}_{\alpha}\left(E^{\prime},\cos\theta\right)\,.$
(22)

Taking into account that $\gamma_{s}\beta_{s}=p_{s}/m_{4}$, Eq. (19) can
finally be written as

$\frac{dN_{\alpha}}{dE}=\mathcal{B}_{\alpha}m_{4}\int_{E^{\prime}_{\mathrm{min}}}^{E^{\prime}_{\mathrm{max}}}dE^{\prime}\int_{E_{s,\rm
min}}^{E_{s,\rm
max}}dE_{s}\frac{1}{p_{s}}\frac{1}{E^{\prime}}\frac{dN_{s}}{dE_{s}}\mathcal{F}_{\alpha}\left(E^{\prime},\cos\theta\right)\,.$
(23)

The minimal and maximal energy $E^{\prime}_{\mathrm{min/max}}$ that an active
neutrino can obtain in the HNL rest frame depend on the kinematic features of
the decay under consideration, and more precisely on the shape of
$\mathcal{F}_{\alpha}$ and on the energetic limits imposed by the boost in Eq.
(20). Indeed, regardless of which distribution $\mathcal{F}_{\alpha}$ the
active neutrino may follow, the energy of the latter in the HNL frame can
never exceed $E/\gamma_{s}\left(1-\beta_{s}\right)$ nor be smaller than
$E/\gamma_{s}\left(1+\beta_{s}\right)$.

As previously mentioned, in the mass range $[0.1-1]$ GeV the two-body decay
$\nu_{4}\rightarrow\nu_{\alpha}\pi^{0}$ and the three-body decay
$\nu_{4}\rightarrow\nu_{\alpha}\nu_{\ell}\bar{\nu}_{\ell}$ in Eq. (18)
dominate over all the other decay modes. Given the different number of final
particles and, consequently, the different kinematic features, the two
channels will be separately discussed in Subsections III.2.1 and III.2.2.

####  III.2.1 Two-body decays

In the two-body decay process $\nu_{4}\rightarrow\nu_{\alpha}\pi^{0}$ the
daughter active neutrino $\nu_{\alpha}$ has in the HNL rest frame a Dirac
delta-like distribution in energy and angle [110, 108, 112, 113]

$\mathcal{F}_{\alpha}^{\mathrm{2b}}=\left.\frac{dN_{\alpha}}{dE^{\prime}d\rm
cos\theta}\right|^{\mathrm{2b}}_{\alpha}=\frac{1}{2}\delta\left(E^{\prime}-\frac{m_{4}^{2}-m_{\pi^{0}}^{2}}{2m_{4}}\right)\,.$
(24)

Therefore, the active neutrino $\nu_{\alpha}$ produced by the HNL together in
a two-body decay with a $\pi_{0}$ will have the following energy spectrum
[108, 113, 60, 112]

$\left.\frac{dN_{\alpha}}{dE}\right|^{\rm
2b}=\frac{\mathcal{B}_{\alpha}m_{4}^{2}}{m_{4}^{2}-m_{\pi^{0}}^{2}}\int_{E_{s,\rm
min}}^{E_{s,\rm max}}dE_{s}\frac{1}{p_{s}}\frac{dN_{s}}{dE_{s}}\,.$ (25)

The decay rate for this two-body channel is [109]

$\Gamma_{\alpha}\left(\nu_{4}\rightarrow\nu_{\alpha}\pi^{0}\right)=2\frac{G_{F}^{2}m_{4}^{3}}{32\pi}f_{\pi}^{2}|U_{\alpha
4}|^{2}\left[1-\left(\frac{m_{\pi_{0}}}{m_{4}}\right)^{2}\right]^{2}\,,$ (26)

where the factor 2 takes into account that for Majorana neutrinos the charged-
conjugated channel $\nu_{4}\rightarrow\bar{\nu}_{\alpha}\pi^{0}$ is open as
well [109, 108, 114], and yields the same value as the main channel:
$\Gamma_{\alpha}\left(\nu_{4}\rightarrow\nu_{\alpha}\pi^{0}\right)=\Gamma_{\alpha}\left(\nu_{4}\rightarrow\bar{\nu}_{\alpha}\pi^{0}\right)$.
In Eq. (26), we assume the meson decay constant $f_{\pi}=130$ MeV [115, 109].
The branching ratio $\mathcal{B}_{\alpha}$ in Eq. (25) has been computed
following [109] considering all the relevant decay processes between 0.1 and
1.0 GeV.

Finally, in a two-body decay the minimal energy of the HNL in the laboratory
frame is obtained by maximizing the boost in Eq. (20) ($\cos{\theta}$ = 1),
giving [108]

$\displaystyle E_{s,\rm
min}=\frac{m_{4}^{2}-m_{\pi_{0}}^{2}}{4E}+\frac{Em_{4}^{2}}{m_{4}^{2}-m_{\pi_{0}}^{2}}\,.$
(27)

The maximal energy of the HNL, instead, has been set to $E_{s,\rm
max}=10^{5}~{}T^{\mathrm{i}n}_{\mathrm{PBH}}$.

####  III.2.2 Three-body decays

In the case of the three body decay
$\nu_{4}\rightarrow\nu_{\alpha}\nu_{\ell}\bar{\nu}_{\ell}$, the produced
active neutrino and antineutrino follow two different angular and energetic
distributions in the HNL rest frame [108, 112, 110]

$\displaystyle\mathcal{F}_{\mathrm{I},\rm\alpha}^{\mathrm{3b}}\left(E^{\prime}\right)=\left.\frac{dN_{\alpha}}{dE^{\prime}d\rm
cos\theta}\right|^{\mathrm{3b}}_{\mathrm{I}}=\frac{1}{2}16\frac{E^{\prime
2}}{m_{4}^{3}}\left(3-4\frac{E^{{}^{\prime}}}{m_{4}}\right)\,,$ (28)
$\displaystyle\mathcal{F}^{\mathrm{3b}}_{\mathrm{II},\rm\alpha}\left(E^{\prime}\right)=\left.\frac{dN_{\alpha}}{dE^{\prime}d\rm
cos\theta}\right|^{\mathrm{3b}}_{\mathrm{II}}=\frac{1}{2}96\frac{E^{\prime
2}}{m_{4}^{3}}\left(1-2\frac{E^{\prime}}{m_{4}}\right)\,.$ (29)

Notice that $\mathcal{F}_{\mathrm{I},\rm\alpha}^{\mathrm{3b}}$ and
$\mathcal{F}^{\mathrm{3b}}_{\mathrm{II},\rm\alpha}$ are double differential
distributions in energy $E^{\prime}$ and angle $\cos{\theta}$ [or
equivalently, in $E$ and $\cos{\theta}$, see Eq. (20)]. For this reason, they
require to be normalized to 1 when integrated over both energy and angle.
However, we assume that the decay products are isotropically distributed [111,
108]. Therefore, when integrated only over the energy $E^{\prime}$,
$\mathcal{F}_{\mathrm{I},\rm\alpha}^{\mathrm{3b}}$ and
$\mathcal{F}^{\mathrm{3b}}_{\mathrm{II},\rm\alpha}$ in Eqs. (28) and (29) must
be normalized to 1/2 in the integration range $[0,m_{4}/2]$ [110], which is
the only kinematically allowed energy range for the active neutrino in the HNL
rest frame.

The active neutrinos $\nu_{\alpha}$ and $\nu_{\ell}$ produced in the decay
process $\nu_{4}\rightarrow\nu_{\alpha}\nu_{\ell}\bar{\nu}_{\ell}$ follow
$\mathcal{F}_{\mathrm{I},\rm\alpha}^{\mathrm{3b}}$, while
$\mathcal{F}^{\mathrm{3b}}_{\mathrm{II},\rm\alpha}$ describes the spectrum of
the antineutrino $\bar{\nu}_{\ell}$. The opposite applies to the charge-
conjugated process
$\nu_{4}\rightarrow\bar{\nu}_{\alpha}\bar{\nu}_{\ell}\nu_{\ell}$, where
$\mathcal{F}_{\mathrm{I},\rm\alpha}^{\mathrm{3b}}$ traces the spectrum of the
antineutrinos $\bar{\nu}_{\alpha},\bar{\nu}_{\ell}$ whereas the neutrino
$\nu_{\ell}$ follows the distribution
$\mathcal{F}^{\mathrm{3b}}_{\mathrm{II},\rm\alpha}$ [108, 110].

As previously considered for the two-body decay of Subsec. III.2.1, we will
focus on scenarios in which only one neutrino species $\alpha$ mixes with the
HNL. Under this assumption and given the above considerations, the total
angular and energetic distribution for the two neutrinos with flavor $\ell$ —
that do not mix with the HNL — is given by [108]

$\mathcal{F}_{\ell}^{\mathrm{3b}}\left(E^{\prime}\right)=\frac{1}{8}\left(\mathcal{F}^{\mathrm{3b}}_{\mathrm{I}}+\mathcal{F}^{\mathrm{3b}}_{\mathrm{II}}\right)\,.$
(30)

The active neutrino with flavor $\alpha$ that mixes with the HNL will instead
have the following distribution function, in the HNL rest frame

$\mathcal{F}_{\alpha}^{\mathrm{3b}}\left(E^{\prime}\right)=\frac{1}{4}\left(3\mathcal{F}^{\mathrm{3b}}_{\mathrm{I}}+\mathcal{F}^{\mathrm{3b}}_{\mathrm{II}}\right)\,.$
(31)

Thus, the active neutrino time-integrated spectrum for the three-body decay is
given by

$\left.\frac{dN_{\alpha(\ell)}}{dE}\right|^{\rm 3b}\
=\mathcal{B}_{\alpha}m_{4}\int_{E_{s},\rm min}^{E_{s},\rm
max}dE_{s}\frac{1}{p_{s}}\frac{dN_{s}}{dE_{s}}\int_{E^{\prime}_{\rm
min}}^{E^{\prime}_{\rm
max}}dE^{\prime}\frac{1}{E^{\prime}}\mathcal{F}_{\alpha(\ell)}^{\mathrm{3b}}\left(E^{\prime}\right)\,,$
(32)

as it can be derived from Eq. (23), for an active neutrino produced in the
three-body decay that mixes with the HNL ($\nu_{\alpha}$) or not
($\nu_{\ell}$).

The minimal and maximal HNL energies in the laboratory frame are set to the
standard values of $E_{s,\rm min}$ = $10^{-5}$ $T^{\rm in}_{\mathrm{PBH}}$ and
$E_{s,\rm max}$ = $10^{5}$ $T^{\rm in}_{\mathrm{PBH}}$, differently from the
two-body case.

The minimal energy of the active neutrino $E^{\prime}_{\rm min}$ in the HNL
rest frame is obtained by inserting the maximal cosine of the angle $\theta$,
i.e. $\cos{\theta}$ = 1, into Eq. (20)

$E^{\prime}_{\rm min}=\frac{E}{\gamma_{s}\left(1+\beta_{s}\right)}\,.$ (33)

On the other hand, the maximal energy of the active neutrino $E^{\prime}_{\rm
max}$ is set to $m_{4}/2$ as it constitutes the upper limit of the
kinematically allowed energy range for $\nu_{\alpha}$ in the HNL rest frame
[110]. Similarly to the two-body decay, the branching ratio of the process
$\mathcal{B}_{\alpha}$ has been computed following Ref. [109].

The partial decay width for the three-body channel reads

$\Gamma_{\alpha}\left(\nu_{4}\rightarrow\nu_{\alpha}\sum_{\ell}\nu_{\ell}\bar{\nu}_{\ell}\right)=\sum_{\ell}\left[\Gamma\left(\nu_{4}\rightarrow\nu_{\alpha}\nu_{\ell}\bar{\nu}_{\ell}\right)+\left(\nu_{4}\rightarrow\bar{\nu}_{\alpha}\nu_{\ell}\bar{\nu}_{\ell}\right)\right]=2\frac{G_{F}^{2}m_{4}^{5}}{64\pi^{3}}|U_{\alpha
4}|^{2}\,.$ (34)

###  III.3 Heavy-mass regime: [0.5-2] TeV

Above the electroweak scale, HNLs mainly decay into bosonic states [58]. In
this scenario, we assume for simplicity that only the muonic neutrino mixes
with the HNL, i.e., we consider only the case 0:1:0. Under this assumption,
the three dominant channels that contribute to the muonic neutrino production
in the [0.5-2] TeV mass range are

$\displaystyle\nu_{4}$ $\displaystyle\to$ $\displaystyle W^{\pm}\mu^{\mp}\,,$
(35) $\displaystyle\nu_{4}$ $\displaystyle\to$ $\displaystyle
Z^{0}\nu_{\mu}\,,$ $\displaystyle\nu_{4}$ $\displaystyle\to$ $\displaystyle
H^{0}\nu_{\mu}\,.$

The produced bosons and muons will subsequently hadronize and/or decay
producing additional fluxes of visible $\nu_{\mu}$, similarly to the secondary
contributions already present in the SM and previously described in Sec. II.2.
In the following, we will dub the intermediate states of the previous
processes, $W^{\pm}$, $\mu^{\pm}$, $Z^{0}$ and $H^{0}$, as “i.s.”.

Similarly to the light-mass regime [see Eq. (23)], the secondary spectrum of
$\nu_{\mu}$ from heavy-mass HNLs can be obtained by a convolution of the HNL
primary spectrum, $\frac{dN_{s}}{dE_{s}}$, with the spectrum of $\nu_{\mu}$
produced at a fixed parent energy: $\frac{dN}{dE^{\prime}}(\nu_{4}\to
i.s.\to\nu_{\mu})$. The latter corresponds to the energetic and angular
distribution that the muon neutrino can have in the HNL rest frame, and will
be evaluated using numerical tools.

To compute such energetic and angular distributions, we use PPPC4DMID [116], a
code that relies on high-statistics simulations and the PYTHIA [117] event
generator. As mentioned before, given the observation time frame under
consideration, the production of HNLs with masses above 1 TeV is heavily
suppressed. Although PPPC4DMID aims at computing the secondary emission from
DM annihilation/decay, it can be used for PBH evaporation following few steps.

We compute the muonic neutrino spectra from each bosonic or leptonic
intermediate state separately, i.e. $\frac{dN}{dE^{\prime}}(\nu_{4}\to
W/Z/H\to\nu_{\mu})$ and
$\frac{dN}{dE^{\prime}}(\nu_{4}\to\mu/\nu_{\mu}\to\nu_{\mu})$ from DM decay
with a centre-of mass energy equal to the energy of the intermediate states in
the HNL rest frame. Afterwards, we combine the resulting neutrino spectra for
each HNL decay channel to obtain the overall neutrino spectrum.

For the sake of comparison and to test our method, we computed the total muon
neutrino spectrum $\frac{dN}{dE^{\prime}}(\nu_{4}\to i.s.\to\nu_{\mu})$ for
$m_{4}=1$ TeV using another tool, HDMSpectra [90] and checked that the two
methods provide same results444Contrarily to PPPC4DMID, HDMSpectra’s range of
validity is above the scale of electroweak symmetry breaking, and in
particular $m_{\mathrm{DM}}\gtrsim$ TeV. Differences in the spectra among the
two codes due to the implementation of electroweak corrections turn out not to
affect the overall spectra of our interest..

The resulting time-integrated spectrum in the laboratory frame of active
muonic neutrinos from the decay of HNLs with masses above the electroweak
scale is

$\frac{dN_{\nu_{\mu}}}{dE}=\sum_{i.s.}\mathcal{B}(\nu_{4}\rightarrow
i.s.)\;m_{4}\int_{E_{s},\rm min}^{E_{s},\rm
max}dE_{s}\frac{1}{p_{s}}\frac{dN_{s}}{dE_{s}}\int_{E^{\prime}_{\rm
min}}^{E^{\prime}_{\rm
max}}dE^{\prime}\frac{1}{E^{\prime}}\frac{dN}{dE^{\prime}}\left(\nu_{4}\rightarrow
i.s.\rightarrow\nu_{\mu}\right)\,,$ (36)

where $\mathcal{B}(\nu_{4}\rightarrow i.s.)$ is the branching ratio of the
process $\nu_{4}\to i.s.$ computed taking into account all the heavy HNL decay
modes of Eq. (35) in the total decay width.

The partial decay widths for each of these processes have been taken from Ref.
[58] and read

$\displaystyle\Gamma\left(\nu_{4}\to\mu W_{L}\right)$
$\displaystyle=2\frac{g^{2}}{64\pi M^{2}_{W}}|U_{\mu
4}|^{2}m_{4}^{3}\left[1-\left(\frac{M_{W}}{m_{4}}\right)^{2}\right]^{2}\,,$
$\displaystyle\Gamma\left(\nu_{4}\to\mu W_{T}\right)$
$\displaystyle=2\frac{g^{2}}{32\pi}|U_{\mu
4}|^{2}m_{4}\left[1-\left(\frac{M_{W}}{m_{4}}\right)^{2}\right]^{2}\,,$
$\displaystyle\Gamma\left(\nu_{4}\to\nu_{\mu}Z^{0}_{L}\right)$
$\displaystyle=\frac{g^{2}}{64\pi M^{2}_{Z}}|U_{\mu
4}|^{2}m_{4}^{3}\left[1-\left(\frac{M_{Z}}{m_{4}}\right)^{2}\right]^{2}\,,$
$\displaystyle\Gamma\left(\nu_{4}\to\nu_{\mu}Z^{0}_{T}\right)$
$\displaystyle=\frac{g^{2}}{32\pi\cos^{2}\theta_{W}}|U_{\mu
4}|^{2}m_{4}\left[1-\left(\frac{M_{Z}}{m_{4}}\right)^{2}\right]^{2}\,,$
$\displaystyle\Gamma\left(\nu_{4}\to\nu_{\mu}H^{0}\right)$
$\displaystyle=\frac{g^{2}}{64\pi M^{2}_{H}}|U_{\mu
4}|^{2}m_{4}^{3}\left[1-\left(\frac{M_{H}}{m_{4}}\right)^{2}\right]^{2}\,,$
(37)

where $M_{W,Z}$ are respectively the mass of W and Z bosons, $\cos\theta_{W}$
is the cosine of the Weinberg angle, and $g$ is the weak coupling. The Z and W
unpolarized decay widths are obtained by summing the decay widths into the
transverse $X_{T}$ and longitudinal $X_{L}$ components, $X=\\{W,Z\\}$, of Eq.
(37) [116]. In the decay widths into a W boson, the factor 2 accounts for the
two possible charged channels $\nu_{4}\to\mu^{\pm}W^{\pm}$.

###  III.4 Summary of neutrino spectra from HNL decays

The overall muon neutrino time-integrated spectrum expected from an
evaporating PBH will be given by the sum of all SM contributions
$\left.\frac{dN_{\nu_{\mu}}}{dE}\right|_{\mathrm{SM}}$ and the secondary
contributions from HNL decay modes
$\left.\frac{dN_{\nu_{\mu}}}{dE}\right|_{\mathrm{HNL}}$ described above

$\left.\frac{dN_{\nu_{\mu}}}{dE}\right|_{\mathrm{tot}}=\left.\frac{dN_{\nu_{\mu}}}{dE}\right|_{\mathrm{SM}}+\left.\frac{dN_{\nu_{\mu}}}{dE}\right|_{\mathrm{HNL}}\,.$
(38)

More specifically, in the case of HNLs with masses in the range $[0.1-1]$ GeV
the muon-neutrino time-integrated spectrum from HNL decays will receive
contributions from both the two- and three-body decays

$\left.\frac{dN_{\nu_{\mu}}}{dE}\right|_{\mathrm{HNL}}=\left.\frac{dN_{\nu_{\mu}}}{dE}\right|^{\mathrm{2b}}+\left.\frac{dN_{\nu_{\mu}}}{dE}\right|^{\mathrm{3b}}\,,$
(39)

whereas the spectrum of muon neutrinos from heavy HNLs with masses in the
$[0.5-2]$ TeV range is directly given in Eq. (36).

Figure 1: Total time-integrated spectrum of muon neutrinos expected at the
Earth from the evaporation of a PBH multiplied by squared-energy, for an
observation time of $\tau=100$ s. In each panel, we show the SM-only
contribution (black, dashed) and the SM + HNL spectrum for different HNL
benchmark masses (solid, color) and mixings.

We show in Fig. 1 the total muon neutrino time-integrated spectrum multiplied
by energy-squared including both SM and HNL contributions as function of the
muon neutrino energy. The upper left (1:0:0), upper right (0:1:0) and lower
left (0:0:1) panels refer to the expected time-integrated spectrum in the
presence of low-mass HNLs, for three different mixing scenarios. Namely, we
fix two HNL benchmark masses $m_{4}=0.2$ (cyan solid) and $m_{4}=1$ (pink
solid) GeV. The black, dashed curve indicates the SM-only contribution. Let us
recall also that we are considering an observation time of $\tau=100$ s. In
all panels, the contribution from the primary neutrinos is observed as a bump
at energies of $E\sim{10^{4}~{}{\rm GeV}}$.

We can see that, regardless of the mixing mode, the light HNLs induce a
significant increase in the expected muon neutrino time-integrated spectrum
between $10^{6}$ and $10^{7}$ GeV with respect to the SM contribution. For
$m_{4}=0.2$ GeV, the HNLs mainly contribute to the muon-neutrino time-
integrated spectrum through the three-body decay given in Eq. (32).

Higher HNLs masses instead give rise to additional two-body decay
contributions that more sizably affect the expected time-integrated spectrum.
However, such increase is non-linear, as at higher HNLs masses other decay
channels without neutrinos as final states open and the branching ratios of
the two and three-body decay are dumped. Eventually, the expected spectrum for
$m_{4}=1$ GeV shows mild differences with respect to the $m_{4}=0.2$ GeV case.

Finally, the lower-right panel of Fig. 1 shows the muon-neutrino time-
integrated spectrum with contributions from HNLs in the mass range $[0.5-2]$
TeV. In particular, we fix two HNL benchmark masses $m_{4}=500$ (blue solid)
and $m_{4}=2000$ (orange solid) GeV and we consider the mixing scenario 0:1:0.
Similarly to the light-mass regime, HNLs decays mainly contribute at high
energies, between $10^{6}$ and $10^{7}$ GeV, and the resulting spectra are
similar for both benchmark masses.

##  IV Statistical Analysis

Neutrino telescopes like IceCube [71, 72] may be able to observe exploding
PBHs as transient point sources. We are interested in evaluating the
signatures at IceCube of muon neutrinos produced from the last stages of a PBH
evaporation, more specifically during an observation time of 100 s. Muon-
tracks events at IceCube benefit from a high-quality angular resolution [50]
that allows to pinpoint with precision the declination angle, i.e. the angle
with respect to the ground, of the incoming astrophysical muon neutrinos. We
evaluate the fluence of $\nu_{\mu}$ from the last instants of an evaporating
PBH located at a distance $d_{\mathrm{PBH}}$ from Earth as

$F_{\nu_{\mu}}(E)=\frac{1}{4\pi
d_{\mathrm{PBH}}^{2}}\left.\frac{dN_{\nu_{\mu}}}{dE}\right|_{\mathrm{tot}}\,,$
(40)

where $\left.\frac{dN_{\nu_{\mu}}}{dE}\right|_{\mathrm{tot}}$ takes into
account all HNL and SM contributions as defined in Eq. (38).

The total number of muon-neutrino events expected at IceCube, for a given
declination angle $\delta$, is obtained as

$N_{\nu_{\mu}}(\delta)=\int_{E_{\mathrm{min}}}^{E_{\mathrm{max}}}F_{\nu_{\mu}}(E)\mathcal{A}_{\mathrm{eff}}(E,\delta)\,,$
(41)

where the integration limits are set to ${E_{\mathrm{min}}}=100$ GeV, that is
the energy threshold of IceCube, and ${E_{\mathrm{max}}}=10^{8}$ GeV,
corresponding to the energy above which the total muon-neutrino fluence
becomes negligible. IceCube’s effective area,
$\mathcal{A}_{\mathrm{eff}}(E,\delta)$, indicates the efficiency to observe an
astrophysical neutrino flux as a function of energy and declination555The
declination angle $\delta$ is related to the azimuthal angle $\zeta$ through
$\cos\zeta=\cos\left(\delta+\pi/2\right)$. . We use the publicly available
effective area given in Fig. 1 in the supplemental material of [72]666Notice
that the effective area strongly depends on the analysis cuts and detector
effects. We rely on the effective area for the IC86 2012-2018 event selection,
corresponding to an analysis sensitivity for a point-like neutrino source that
most resembles our case. and we fix the declination angle to be in the range
$[30\degree<\delta<90\degree]$. This choice corresponds to muon neutrinos from
a portion of the Northern sky and allows us to neglect the atmospheric muon
background with more reliability, since almost all atmospheric muons from the
Northern hemisphere are filtered out by the Earth.

Figure 2: Expected number of muon-neutrino events at IceCube as a function of
the HNL mass, from the last stages (100 s) of an evaporating PBH located at a
distance $d_{\mathrm{PBH}}=10^{-4}$ pc from Earth and at a declination angle
$[30\degree<\delta<90\degree]$. The black dashed curve corresponds to the SM-
only case.

We present in Fig. 2 the number of muon neutrino events expected at IceCube
from an evaporating PBH located at $d_{\mathrm{PBH}}=10^{-4}$ pc from Earth,
as a function of the HNL mass. As before, the PBH burst is assumed to last 100
s. The different colored curves represent the number of muon-tracks events
from HNL in the $[0.1,1]$ GeV range for three mixing scenarios, 1:0:0 (cyan),
0:1:0 (red), and 0:0:1 (green). The black dashed curve indicates the SM-only
contribution. The behavior of the curves with HNLs is mainly determined by the
oscillations among active states and the interplay between the two- and three-
body decay contributions. The dip visible around $m_{4}\sim 0.15$ GeV marks
the transition from the regime where the three-body decay
$\nu_{4}\to\nu_{\alpha}\nu_{\ell}\bar{\nu}_{\ell}$ dominates to the regime
where $\nu_{4}\to\nu_{\alpha}\pi^{0}$ is instead more important. A similar
dependence of $N_{\nu_{\mu}}$ on the HNL mass can be seen in muon-tracks
events where the muon neutrino does not directly mix with the HNL, i.e. in the
1:0:0 and 0:0:1 mixing modes. However, since $\theta_{23}$ is close to
maximal, the number of events for the 0:0:1 case is larger than for 1:0:0.

Concerning HNLs with masses above the electroweak scale, we consider only the
mixing 0:1:0 for illustration and restrict our analysis to HNLs with masses
above 500 GeV, to ensure the dominance of the bosonic channels in Eq. (35). In
this scenario, the number of expected muon-tracks events in the presence of
HNLs hardly depends on the HNL mass, since the mass dependence in the
different channels of Eq. (37) is approximately the same and $m_{4}$
eventually simplifies in $\mathcal{B}_{\alpha}$. We predict
$N_{\nu_{\mu}}^{\mathrm{SM+HNL}}\sim 1360$, independently of $m_{4}$, in the
mass range $[0.5,2]$ TeV.

We have also quantified how much using different ranges of declination angles
affect the number of muon neutrinos expected at IceCube. In particular,
declination angles between $-5\degree$ and $30\degree$, i.e. covering the
remaining portion of the Northern sky, would imply $8\%$ more events. On the
other hand, neutrino fluxes incoming from the Southern sky, i.e.,
$[-90\degree<\delta<-5\degree]$, are almost completely shadowed by the
background of atmospheric muons leading to $N_{\nu_{\mu}}<100$.

Next, to estimate the sensitivity at IceCube we rely on the following
$\chi^{2}$ test-statistics

$\chi^{2}=\left(\dfrac{N^{\mathrm{SM+HNL}}_{\nu_{\mu}}-N^{\mathrm{SM}}_{\nu_{\mu}}}{\sigma_{\nu_{\mu}}}\right)^{2}+\left(\dfrac{\alpha}{\sigma_{\alpha}}\right)^{2}\,,$
(42)

where $N^{\mathrm{SM+HNL}}_{\nu_{\mu}}$ is the total number of $\nu_{\mu}$
events expected from the SM (primary + secondary) and HNL contributions,
whereas $N^{\mathrm{SM}}_{\nu_{\mu}}$ indicates the SM-only contribution. The
main source of background could be high-energy atmospheric neutrino events,
however the short duration of the PBH evaporation and the choice of
declination angle make this observation essentially background-free [50, 46].
Indeed, atmospheric neutrino events are estimated to be $\mathcal{O}(10^{-4})$
for an observation time of 100 s [50, 118, 119]. Moreover, we have checked
that the inclusion of a nuisance parameter $\alpha$ related to the measurement
of $d_{\mathrm{PBH}}$ with an uncertainty of $\sigma_{\alpha}=10\%$ does not
significantly affect the estimated sensitivity at IceCube. For this reason, in
the rest of our analysis we do not include any nuisance parameter. Given these
premises, we fix $\sigma_{\nu_{\mu}}$ to be dominated by statistical errors,
$\sigma_{\nu_{\mu}}=N^{\mathrm{SM}}_{\nu_{\mu}}$.

It should be noted that the estimated sensitivity at IceCube highly depends on
the chosen declination angle. The $\chi^{2}$ at $d_{\rm PBH}$ = $10^{-4}$ is
increased by 10% when the declination angle is set between $-5\degree$ and
$30\degree$ and drops below $\chi^{2}<0.2$ for muon-tracks events from the
Southern sky.

##  V Results

We compute the sensitivity at IceCube using the $\chi^{2}$ function defined in
Eq. (42) and varying the following free parameters: $m_{4}$ and
$d_{\mathrm{PBH}}$. While relevant to obtain the partial decay widths of the
HNL, $|U_{\alpha 4}|^{2}$ does not affect the estimation of the muon-track
events. Indeed, when only one active neutrino species mixes with the HNL, as
in the mixing scenarios considered in this work, only the branching ratio
enters in the computation of $N_{\nu_{\mu}}$ and the dependence on $|U_{\alpha
4}|^{2}$ cancels out. However, $|U_{\alpha 4}|^{2}$ enters as a relevant
ingredient in the computation of the HNL lifetime. Therefore, if this mixing
is small, the decay length may exceed $d_{\rm PBH}$, causing the HNL to decay
after traversing the Earth. This dependence on the mixing will put a lower
bound to our sensitivity as shown below.

Figure 3: IceCube sensitivity to HNLs from a PBH burst lasting 100 s, in terms
of the HNL mass and the PBH distance to Earth, for different scenarios for the
HNL-active neutrino mixing. See text for more details.

We show in Fig. 3 the $1,2,3,5\sigma$ sensitivities (lighter to darker shades)
in terms of $m_{4}$ and $d_{\rm PBH}$. Each panel corresponds to a different
mixing scenario, as indicated in its title, and we assume the mixing to be
large enough for the HNL decay to occur before reaching the Earth. As
anticipated, the $\chi^{2}$ does not depend on $|U_{\alpha 4}|^{2}$. In the
case of a HNL with mass in the $[0.1,1]$ GeV range, the mixing scenario 1:0:0
leads to the smallest number of muon neutrino events, as was evident already
from Fig. 2. This translates into sensitivities shifted towards smaller
distances, as in order to obtain a sizable signal the PBH must be located
closer to Earth. In the other two scenarios, 0:1:0 and 0:0:1, we see that if
the PBH explosion occurs at a distance $d_{\mathrm{PBH}}\sim 10^{-4}$ pc777For
the sake of comparison, this distance is of the same order as the Uranus-Sun
distance., IceCube could probe basically all the mass range $m_{4}=[0.1,1]$
GeV at $1\sigma$ CL. HNLs with masses above the electroweak scale (lower-right
panel in Fig. 3) would produce instead larger fluxes of muon neutrinos, and
would be visible even if occurring at a larger distance.

In order to compare with previous bounds, we depict in Fig. 4 our 90% CL
sensitivities in the plane ($|U_{\alpha 4}|^{2}$, $m_{4}$) as colored shaded
regions. We show a different mixing scenario in each panel, and we fix some
benchmark values for the PBH distance $d_{\mathrm{PBH}}$, as indicated in the
legends. Although not directly dependent on $|U_{\alpha 4}|^{2}$ as already
explained, our exclusion regions show a diagonal cut. As mentioned before, we
impose a lower limit on the HNL lifetime, or, equivalently, an upper limit on
the HNL decay length that cannot exceed $d_{\mathrm{PBH}}$. This constraint
excludes mixings below a certain value of $|U_{\alpha 4}|^{2}$, depending on
the HNL mass. Moreover, in each panel we superimpose existing bounds (see e.g.
Ref. [60] for a review), indicated as a dark-blue shaded area. In particular,
for scenario 1:0:0 bounds from NA62 [120], T2K [121], PiENU [122], BEBC [123]
and PS191 [124] are relevant; for 0:1:0, T2K [121], MicroBooNE [125], NuTeV
[126], E949 [127]. Finally, for 0:0:1 the most important bounds in the
considered mass range are given by T2K [121], CHARM [128], and constraints
from IceCube [129] looking for low-energy “double-bang” events.

Our results show that the observation of muon-track events at IceCube from the
explosion of a PBH located at a distance $d_{\mathrm{PBH}}=10^{-4}$ pc would
allow to explore the HNL scenario at masses $0.3~{}(0.1)\lesssim m_{4}\lesssim
0.6~{}(1)$ GeV, in the case where the HNL mixes only with the muon (tau)
flavor, testing values of mixing below those currently excluded by other
searches. In the case of mixing only with the electron flavor, the expected
flux of muon neutrinos would be smaller and therefore would require the PBH
explosion to occur as close as $d_{\mathrm{PBH}}=2\times 10^{-5}$ pc to Earth
to be able to probe part of the parameter space not yet excluded.

Moving to $m_{4}>500$ GeV instead, the expected fluxes of muon neutrinos are
larger and the PBH would be visible already at a distance of
$d_{\mathrm{PBH}}=2.5\times 10^{-4}$ pc. In this case, the inferred bounds
would be much more stringent as they would cover all the parameter space down
to very tiny values of $|U_{\alpha 4}|^{2}$. Indeed, such heavy HNLs would be
very short-lived and decay right after their production from PBH evaporation.
As a consequence, the HNL decay lifetime bound would only apply to mixings
$|U_{\alpha 4}|^{2}\lesssim 10^{-29}$, as visible in the lower-right panel of
Fig. 4.

Figure 4: Expected IceCube sensitivity (at $90\%$ CL) to the mixing matrix
elements $|U_{\alpha 4}|^{2}$ as a function of the HNL mass, for a PBH burst
of $\tau=100$ s. Each panel correspond to a different mixing scenario, in
which the HNL only mixes with one of the active neutrinos with flavor
$\alpha$, while the other two mixings are set to zero. Different colors
correspond to different PBH distances to Earth. Dark blue shaded areas denote
existing constraints [130, 60].

As a final comment, let us recall that all sensitivities have been obtained
assuming a burst of $\tau=100$ s, corresponding to a PBH with mass $M\sim
6\times 10^{9}$ g at the beginning of the explosion. Should we consider
smaller observation time frames, and hence smaller $M$, our results would
change as a $90\%$ CL detection would require the explosion to occur closer to
Earth. To quantify this effect, for an observation time of $\tau=10$ s, i.e.,
$M\sim 2.9\times 10^{9}$, depending on the mixing scenario, the PBH should be
$\sim$ 30-50% closer to obtain the same $90\%$ CL sensitivities.

##  VI Conclusions

Several theories of the early Universe predict the presence of a population of
black holes of primordial origin, which in turn may be linked to the DM
mystery. As theorized by Hawking, these objects are expected to evaporate,
emitting radiation with thermal spectra. The characteristics of such an
emission are intimately related to the properties of emitted particles (mass,
charge, spin and number of degrees of freedom) and therefore provide valuable
opportunities to probe the existence of new physics. The final moments of this
evaporation phase are particularly relevant, as the emission increases while
the PBH loses mass and raises its temperature, leading to the emission of
heavier degrees of freedom.

In this work, we investigated a scenario in which one light PBH explodes, i.e.
ends its evaporation phase close to Earth thus giving rise to a “burst” of
visible particles. We were particularly interested in a BSM scenario in which
the SM is extended with one HNL that mixes with active neutrinos. Such a new
sterile state will be emitted in the evaporation process and consequently
decay into active neutrinos, sizably contributing to the expected flux of
observable neutrinos at Earth.

We estimated the spectra of muon neutrinos expected at a neutrino telescope
like IceCube, taking into account all primary and secondary contributions from
the PBH explosion, including the HNL ones. Under the assumption that the PBH
explosion occurs sufficiently close to Earth to be visible at IceCube and
lasts 100 s, we inferred sensitivities on the relevant HNL parameter space, in
terms of the HNL mass and mixing with the active flavors. Depending on how
close the explosion occurs, IceCube could be able to set stringent
constraints, even improving existing ones from other facilities, in some
ranges of the HNL mass and relevant mixing.

Certainly, our estimates are strongly dependent on the assumptions made for
the analysis, namely current IceCube’s effective area, its angular resolution
and the capacity of reducing backgrounds. Future upgrades are anticipated to
improve the effective area by a factor $\sim 5$ [131], thus implying a
substantial enhancement in the detection of high-energy muon tracks. The
development of new analysis and detection techniques will further enhance the
possibility of observing such events, and eventually improve the sensitivities
to BSM physics, translating into observable signals even if the explosion
occurs at farther distances. In the near future, other neutrino telescopes
like KM3NeT [132], P-ONE [133], and Baikal-GVD [134] will also relevantly
contribute to the search for PBH explosions in the proximity of Earth, being
sensitive to different directions in the sky and thus allowing to provide a
more complete picture.

Finally, let us comment that the rapid evaporation of these objects in their
final instants will also give rise to high-energy photon bursts definitely
within the reach of current and forthcoming gamma-ray experiments, in
particular Cherenkov telescope arrays [37]. A clear identification of the PBH
origin of such a burst will call for a multimessenger approach, thus making
observations with both neutrinos and photons essential at the scope of a
precise analysis. On the other hand, in the specific BSM scenario considered
here, we do not expect sizable contributions into photons from HNL decays
compared to the SM one. For this reason, we expect neutrino telescopes to be
the preferred facilities to infer bounds on the HNL scenario from PBH
explosions. In any case, we leave more investigations on the photon emission
in the HNL scenario for future work.

## Acknowledgments

We are grateful to A. Arbey for support with BlackHawk. Y.F.P.G. would like to
thank the warm hospitality of the Astroparticles and High Energy Physics group
of IFIC during the visit that marked the inception of this work. V.D.R. and
A.T. acknowledge financial support by the CIDEXG/2022/20 grant (project
“D’AMAGAT”) funded by Generalitat Valenciana and by the Spanish grants
CNS2023-144124 (MCIN/AEI/10.13039/501100011033 and “Next Generation EU”/PRTR)
and PID2020-113775GB-I00 (MCIN/AEI/10.13039/501100011033). Y.F.P.G. has been
funded by the UK Science and Technology Facilities Council (STFC) under grant
ST/T001011/1. This project has received funding/support from the European
Union’s Horizon 2020 research and innovation programme under the Marie
Skłodowska-Curie grant agreement No 860881-HIDDeN.

## References

  * [1] Planck Collaboration, N. Aghanim et al. Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO].
  * [2] M. Schumann J. Phys. G 46 no. 10, (2019) 103003, arXiv:1903.03026 [astro-ph.CO].
  * [3] J. Billard et al. arXiv:2104.07634 [hep-ex].
  * [4] B. Carr and F. Kuhnel Ann. Rev. Nucl. Part. Sci. 70 (2020) 355–394, arXiv:2006.02838 [astro-ph.CO].
  * [5] S. Bird et al. Phys. Dark Univ. 41 (2023) 101231, arXiv:2203.08967 [hep-ph].
  * [6] A. M. Green  2, 2024.  arXiv:2402.15211 [astro-ph.CO].
  * [7] Y. B. . N. Zel’dovich, I. D. Soviet Astron. AJ (Engl. Transl. ), 10 (1967) 602.
  * [8] S. Hawking Mon. Not. Roy. Astron. Soc. 152 (1971) 75.
  * [9] B. J. Carr and S. W. Hawking Mon. Not. Roy. Astron. Soc. 168 (1974) 399–415.
  * [10] G. F. Chapline Nature 253 no. 5489, (1975) 251–252.
  * [11] LIGO Scientific, Virgo Collaboration, B. P. Abbott et al. Phys. Rev. Lett. 116 no. 6, (2016) 061102, arXiv:1602.03837 [gr-qc].
  * [12] LIGO Scientific, Virgo Collaboration, R. Abbott et al. Astrophys. J. Lett. 913 no. 1, (2021) L7, arXiv:2010.14533 [astro-ph.HE].
  * [13] S. W. Hawking Nature 248 (1974) 30–31.
  * [14] S. W. Hawking Commun. Math. Phys. 43 (1975) 199–220. [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * [15] B. J. Carr Astrophys. J.  206 (May, 1976) 8–25.
  * [16] R. Lehoucq, M. Casse, J. M. Casandjian, and I. Grenier Astron. Astrophys. 502 (2009) 37, arXiv:0906.1648 [astro-ph.HE].
  * [17] E. L. Wright Astrophys. J. 459 (1996) 487, arXiv:astro-ph/9509074.
  * [18] A. Arbey and J. Auffinger Eur. Phys. J. C 79 no. 8, (2019) 693, arXiv:1905.04268 [gr-qc].
  * [19] G. Ballesteros, J. Coronado-Blázquez, and D. Gaggero Phys. Lett. B 808 (2020) 135624, arXiv:1906.10113 [astro-ph.CO].
  * [20] R. Laha, J. B. Muñoz, and T. R. Slatyer Phys. Rev. D 101 no. 12, (2020) 123514, arXiv:2004.00627 [astro-ph.CO].
  * [21] M. Boudaud and M. Cirelli Physical Review Letters 122 no. 4, (Jan, 2019) . http://dx.doi.org/10.1103/PhysRevLett.122.041104.
  * [22] B. Dasgupta, R. Laha, and A. Ray Phys. Rev. Lett. 125 no. 10, (2020) 101101, arXiv:1912.01014 [hep-ph].
  * [23] W. DeRocco and P. W. Graham Phys. Rev. Lett. 123 no. 25, (2019) 251102, arXiv:1906.07740 [astro-ph.CO].
  * [24] R. Laha Phys. Rev. Lett. 123 no. 25, (2019) 251101, arXiv:1906.09994 [astro-ph.HE].
  * [25] S. Wang, D.-M. Xia, X. Zhang, S. Zhou, and Z. Chang Phys. Rev. D 103 no. 4, (2021) 043010, arXiv:2010.16053 [hep-ph].
  * [26] V. De Romeri, P. Martínez-Miravé, and M. Tórtola JCAP 10 (2021) 051, arXiv:2106.05013 [hep-ph].
  * [27] N. Bernal, V. Muñoz Albornoz, S. Palomares-Ruiz, and P. Villanueva-Domingo JCAP 10 (2022) 068, arXiv:2203.14979 [hep-ph].
  * [28] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama Phys. Rev. D 81 (2010) 104019, arXiv:0912.5297 [astro-ph.CO].
  * [29] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama arXiv:2002.12778 [astro-ph.CO].
  * [30] C. Keith, D. Hooper, N. Blinov, and S. D. McDermott arXiv:2006.03608 [astro-ph.CO].
  * [31] P. He and L.-Z. Fang Astrophys. J. Lett. 568 (2002) L1–L4, arXiv:astro-ph/0202218.
  * [32] K. J. Mack and D. H. Wesley arXiv:0805.1531 [astro-ph].
  * [33] J. Auffinger Prog. Part. Nucl. Phys. 131 (2023) 104040, arXiv:2206.02672 [astro-ph.CO].
  * [34] D. N. Page Phys. Rev. D 13 (1976) 198–206.
  * [35] D. N. Page Phys. Rev. D 14 (1976) 3260–3273.
  * [36] J. H. MacGibbon, B. J. Carr, and D. N. Page Phys. Rev. D 78 (2008) 064043, arXiv:0709.2380 [astro-ph].
  * [37] X. Boluna, S. Profumo, J. Blé, and D. Hennings arXiv:2307.06467 [astro-ph.HE].
  * [38] H.E.S.S. Collaboration, J.-F. Glicenstein, A. Barnacka, M. Vivier, and T. Herr in 33rd International Cosmic Ray Conference, p. 0930.  7, 2013.  arXiv:1307.4898 [astro-ph.HE].
  * [39] T. Tavernier, J.-F. Glicenstein, and F. Brun PoS ICRC2019 (2020) 804, arXiv:1909.01620 [astro-ph.HE].
  * [40] A. A. Abdo et al. Astropart. Phys. 64 (2015) 4–12, arXiv:1407.1686 [astro-ph.HE].
  * [41] VERITAS Collaboration, S. Archambault PoS ICRC2017 (2018) 691, arXiv:1709.00307 [astro-ph.HE].
  * [42] Fermi-LAT Collaboration, M. Ackermann et al. Astrophys. J. 857 no. 1, (2018) 49, arXiv:1802.00100 [astro-ph.HE].
  * [43] HAWC Collaboration, A. U. Abeysekara et al. arXiv:1310.0073 [astro-ph.HE].
  * [44] HAWC Collaboration, A. Albert et al. JCAP 04 (2020) 026, arXiv:1911.04356 [astro-ph.HE].
  * [45] T. N. Ukwatta, J. H. MacGibbon, W. C. Parke, K. S. Dhuga, S. Rhodes, A. Eskandarian, N. Gehrels, L. Maximon, and D. C. Morris in 12th Marcel Grossmann Meeting on General Relativity, pp. 1588–1590.  3, 2010.  arXiv:1003.4515 [astro-ph.HE].
  * [46] Y. F. Perez-Gonzalez Phys. Rev. D 108 no. 8, (2023) 083014, arXiv:2307.14408 [astro-ph.HE].
  * [47] H.E.S.S. Collaboration, F. Aharonian et al. JCAP 04 (2023) 040, arXiv:2303.12855 [astro-ph.HE].
  * [48] F. Halzen, B. Keszthelyi, and E. Zas Phys. Rev. D 52 (1995) 3239–3247, arXiv:hep-ph/9502268.
  * [49] IceCube Collaboration, P. Dave and I. Taboada PoS ICRC2019 (2021) 863, arXiv:1908.05403 [astro-ph.HE].
  * [50] A. Capanema, A. Esmaeili, and A. Esmaili JCAP 12 no. 12, (2021) 051, arXiv:2110.05637 [hep-ph].
  * [51] M. Calzà and J. a. G. Rosa arXiv:2312.09261 [hep-ph].
  * [52] T. N. Ukwatta, D. R. Stump, J. T. Linnemann, J. H. MacGibbon, S. S. Marinelli, T. Yapici, and K. Tollefson Astropart. Phys. 80 (2016) 90–114, arXiv:1510.04372 [astro-ph.HE].
  * [53] M. J. Baker and A. Thamm SciPost Phys. 12 no. 5, (2022) 150, arXiv:2105.10506 [hep-ph].
  * [54] M. J. Baker and A. Thamm JHEP 01 (2023) 063, arXiv:2210.02805 [hep-ph].
  * [55] M. Calzà, J. March-Russell, and J. a. G. Rosa arXiv:2110.13602 [astro-ph.CO].
  * [56] M. Calzà and J. a. G. Rosa JHEP 12 (2022) 090, arXiv:2210.06500 [gr-qc].
  * [57] M. Calzà and J. a. G. Rosa arXiv:2311.12930 [gr-qc].
  * [58] A. Atre, T. Han, S. Pascoli, and B. Zhang JHEP 05 (2009) 030, arXiv:0901.3589 [hep-ph].
  * [59] M. Drewes and B. Garbrecht Nucl. Phys. B 921 (2017) 250–315, arXiv:1502.00477 [hep-ph].
  * [60] A. M. Abdullahi et al. J. Phys. G 50 no. 2, (2023) 020501, arXiv:2203.08039 [hep-ph].
  * [61] E. V. Bugaev, M. G. Elbakidze, and K. V. Konishchev Phys. Atom. Nucl. 66 (2003) 476–480, arXiv:astro-ph/0110660.
  * [62] T. Fujita, M. Kawasaki, K. Harigaya, and R. Matsuda Phys. Rev. D 89 no. 10, (2014) 103501, arXiv:1401.1909 [astro-ph.CO].
  * [63] L. Morrison, S. Profumo, and Y. Yu JCAP 05 (2019) 005, arXiv:1812.10606 [astro-ph.CO].
  * [64] A. Ambrosone, R. Calabrese, D. F. G. Fiorillo, G. Miele, and S. Morisi Phys. Rev. D 105 no. 4, (2022) 045001, arXiv:2106.11980 [hep-ph].
  * [65] D. Hooper and G. Krnjaic Phys. Rev. D 103 no. 4, (2021) 043504, arXiv:2010.01134 [hep-ph].
  * [66] Y. F. Perez-Gonzalez and J. Turner Phys. Rev. D 104 no. 10, (2021) 103021, arXiv:2010.03565 [hep-ph].
  * [67] N. Bernal, C. S. Fong, Y. F. Perez-Gonzalez, and J. Turner Phys. Rev. D 106 no. 3, (2022) 035019, arXiv:2203.08823 [hep-ph].
  * [68] R. Calabrese, M. Chianese, J. Gunn, G. Miele, S. Morisi, and N. Saviano arXiv:2311.13276 [hep-ph].
  * [69] R. Calabrese, M. Chianese, J. Gunn, G. Miele, S. Morisi, and N. Saviano Phys. Rev. D 107 no. 12, (2023) 123537, arXiv:2305.13369 [hep-ph].
  * [70] A. Ghoshal, Y. F. Perez-Gonzalez, and J. Turner arXiv:2312.06768 [hep-ph].
  * [71] IceCube Collaboration, M. G. Aartsen et al. Phys. Rev. D 99 no. 3, (2019) 032004, arXiv:1808.07629 [hep-ex].
  * [72] IceCube Collaboration, M. G. Aartsen et al. Phys. Rev. Lett. 124 no. 5, (2020) 051103, arXiv:1910.08488 [astro-ph.HE].
  * [73] S. W. Hawking Phys. Rev. D 14 (1976) 2460–2473.
  * [74] A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian, and A. Tajdini Rev. Mod. Phys. 93 no. 3, (2021) 035002, arXiv:2006.06872 [hep-th].
  * [75] L. Buoninfante, F. Di Filippo, and S. Mukohyama JHEP 10 (2021) 081, arXiv:2107.05662 [hep-th].
  * [76] D. N. Page Phys. Rev. Lett. 71 (1993) 3743–3746, arXiv:hep-th/9306083.
  * [77] D. N. Page JCAP 09 (2013) 028, arXiv:1301.4995 [hep-th].
  * [78] A. Alexandre, G. Dvali, and E. Koutsangelas arXiv:2402.14069 [hep-ph].
  * [79] V. Thoss, A. Burkert, and K. Kohri arXiv:2402.17823 [astro-ph.CO].
  * [80] S. Balaji, G. Domènech, G. Franciolini, A. Ganz, and J. Tränkle arXiv:2403.14309 [gr-qc].
  * [81] J. H. MacGibbon and B. R. Webber Phys. Rev. D 41 (1990) 3052–3079.
  * [82] J. H. MacGibbon Phys. Rev. D 44 (1991) 376–392.
  * [83] A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner Phys. Rev. D 105 no. 1, (2022) 015022, arXiv:2107.00013 [hep-ph].
  * [84] A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner Phys. Rev. D 108 no. 1, (2023) 015005, arXiv:2212.03878 [hep-ph].
  * [85] D. N. Page Phys. Rev. D 16 (1977) 2402–2411.
  * [86] A. Arbey and J. Auffinger Eur. Phys. J. C 81 (2021) 910, arXiv:2108.02737 [gr-qc].
  * [87] C. Doran, A. Lasenby, S. Dolan, and I. Hinder Phys. Rev. D 71 (2005) 124020, arXiv:gr-qc/0503019.
  * [88] S. Dolan, C. Doran, and A. Lasenby Phys. Rev. D 74 (2006) 064005, arXiv:gr-qc/0605031.
  * [89] C. Lunardini and Y. F. Perez-Gonzalez JCAP 08 (2020) 014, arXiv:1910.07864 [hep-ph].
  * [90] C. W. Bauer, N. L. Rodd, and B. R. Webber JHEP 06 (2021) 121, arXiv:2007.15001 [hep-ph].
  * [91] P. Minkowski Phys. Lett. 67B (1977) 421–428.
  * [92] T. Yanagida Conf. Proc. C 7902131 (1979) 95–99.
  * [93] M. Gell-Mann, P. Ramond, and R. Slansky Conf. Proc. C790927 (1979) 315–321, arXiv:1306.4669 [hep-th].
  * [94] R. N. Mohapatra and G. Senjanovic Phys. Rev. Lett. 44 (1980) 912.
  * [95] J. Schechter and J. W. F. Valle Phys.Rev.D 22 (1980) 2227.
  * [96] P. F. de Salas, D. V. Forero, S. Gariazzo, P. Martínez-Miravé, O. Mena, C. A. Ternes, M. Tórtola, and J. W. F. Valle JHEP 02 (2021) 071, arXiv:2006.11237 [hep-ph].
  * [97] F. Capozzi, E. Di Valentino, E. Lisi, A. Marrone, A. Melchiorri, and A. Palazzo Phys. Rev. D 104 no. 8, (2021) 083031, arXiv:2107.00532 [hep-ph].
  * [98] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou JHEP 09 (2020) 178, arXiv:2007.14792 [hep-ph].
  * [99] R. N. Mohapatra and J. W. F. Valle Phys. Rev. D 34 (1986) 1642.
  * [100] E. K. Akhmedov, M. Lindner, E. Schnapka, and J. W. F. Valle Phys. Lett. B 368 (1996) 270–280, arXiv:hep-ph/9507275.
  * [101] E. K. Akhmedov, M. Lindner, E. Schnapka, and J. W. F. Valle Phys. Rev. D 53 (1996) 2752–2780, arXiv:hep-ph/9509255.
  * [102] M. Malinsky, J. C. Romao, and J. W. F. Valle Phys. Rev. Lett. 95 (2005) 161801, arXiv:hep-ph/0506296.
  * [103] R. E. Shrock Phys. Rev. D 24 (1981) 1232.
  * [104] J. C. Helo, S. Kovalenko, and I. Schmidt Nucl. Phys. B 853 (2011) 80–104, arXiv:1005.1607 [hep-ph].
  * [105] K. Bondarenko, A. Boyarsky, D. Gorbunov, and O. Ruchayskiy JHEP 11 (2018) 032, arXiv:1805.08567 [hep-ph].
  * [106] D. A. Bryman and R. Shrock Phys. Rev. D 100 no. 5, (2019) 053006, arXiv:1904.06787 [hep-ph].
  * [107] D. A. Bryman and R. Shrock Phys. Rev. D 100 (2019) 073011, arXiv:1909.11198 [hep-ph].
  * [108] L. Mastrototaro, A. Mirizzi, P. D. Serpico, and A. Esmaili JCAP 01 (2020) 010, arXiv:1910.10249 [hep-ph].
  * [109] P. Coloma, E. Fernández-Martínez, M. González-López, J. Hernández-García, and Z. Pavlovic Eur. Phys. J. C 81 no. 1, (2021) 78, arXiv:2007.03701 [hep-ph].
  * [110] K. Akita, S. H. Im, M. Masud, and S. Yun arXiv:2312.13627 [hep-ph].
  * [111] L. Oberauer, C. Hagner, G. Raffelt, and E. Rieger Astropart. Phys. 1 (1993) 377–386.
  * [112] V. Syvolap arXiv:2301.07052 [hep-ph].
  * [113] L. Mastrototaro, P. D. Serpico, A. Mirizzi, and N. Saviano Phys. Rev. D 104 no. 1, (2021) 016026, arXiv:2104.11752 [hep-ph].
  * [114] IceCube Collaboration, L. Fischer PoS ICHEP2022 (11, 2022) 190.
  * [115] Particle Data Group Collaboration, R. L. Workman et al. PTEP 2022 (2022) 083C01.
  * [116] M. Cirelli, G. Corcella, A. Hektor, G. Hutsi, M. Kadastik, P. Panci, M. Raidal, F. Sala, and A. Strumia JCAP 03 (2011) 051, arXiv:1012.4515 [hep-ph]. [Erratum: JCAP 10, E01 (2012)].
  * [117] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands Comput. Phys. Commun. 191 (2015) 159–177, arXiv:1410.3012 [hep-ph].
  * [118] A. Esmaili, A. Ibarra, and O. L. G. Peres JCAP 11 (2012) 034, arXiv:1205.5281 [hep-ph].
  * [119] M. Honda, M. Sajjad Athar, T. Kajita, K. Kasahara, and S. Midorikawa Phys. Rev. D 92 no. 2, (2015) 023004, arXiv:1502.03916 [astro-ph.HE].
  * [120] NA62 Collaboration, E. Cortina Gil et al. Phys. Lett. B 807 (2020) 135599, arXiv:2005.09575 [hep-ex].
  * [121] T2K Collaboration, K. Abe et al. Phys. Rev. D 100 no. 5, (2019) 052006, arXiv:1902.07598 [hep-ex].
  * [122] PIENU Collaboration, A. Aguilar-Arevalo et al. Phys. Rev. D 97 no. 7, (2018) 072012, arXiv:1712.03275 [hep-ex].
  * [123] R. Barouki, G. Marocco, and S. Sarkar SciPost Phys. 13 (2022) 118, arXiv:2208.00416 [hep-ph].
  * [124] G. Bernardi et al. Phys. Lett. B 203 (1988) 332–334.
  * [125] MicroBooNE Collaboration, P. Abratenko et al. Phys. Rev. D 101 no. 5, (2020) 052001, arXiv:1911.10545 [hep-ex].
  * [126] NuTeV, E815 Collaboration, A. Vaitaitis et al. Phys. Rev. Lett. 83 (1999) 4943–4946, arXiv:hep-ex/9908011.
  * [127] E949 Collaboration, A. V. Artamonov et al. Phys. Rev. D 91 no. 5, (2015) 052001, arXiv:1411.3963 [hep-ex]. [Erratum: Phys.Rev.D 91, 059903 (2015)].
  * [128] CHARM II Collaboration, P. Vilain et al. Phys. Lett. B 343 (1995) 453–458.
  * [129] P. Coloma, P. A. N. Machado, I. Martinez-Soler, and I. M. Shoemaker Phys. Rev. Lett. 119 no. 20, (2017) 201804, arXiv:1707.08573 [hep-ph].
  * [130] P. D. Bolton, F. F. Deppisch, and P. S. Bhupal Dev JHEP 03 (2020) 170, arXiv:1912.03058 [hep-ph].
  * [131] IceCube-Gen2 Collaboration, M. G. Aartsen et al. J. Phys. G 48 no. 6, (2021) 060501, arXiv:2008.04323 [astro-ph.HE].
  * [132] KM3Net Collaboration, S. Adrian-Martinez et al. J. Phys. G 43 no. 8, (2016) 084001, arXiv:1601.07459 [astro-ph.IM].
  * [133] P-ONE Collaboration, M. Agostini et al. Nature Astron. 4 no. 10, (2020) 913–915, arXiv:2005.09493 [astro-ph.HE].
  * [134] BAIKAL Collaboration, I. A. Belolaptikov et al. Astropart. Phys. 7 (1997) 263–282.