# Identifying Spin Properties of Evaporating Black Holes through
Asymmetric Neutrino and Photon Emission

Yuber F. Perez-Gonzalez<EMAIL_ADDRESS>Institute for
Particle Physics Phenomenology, Durham University, South Road DH1 3LE, Durham,
United Kingdom

###### Abstract

Kerr black holes radiate neutrinos in an asymmetric pattern, preferentially in
the lower hemisphere relative to the black hole’s rotation axis, while
antineutrinos are predominantly produced in the upper hemisphere. Leveraging
this asymmetric emission, we explore the potential of high-energy,
$E_{\nu}\gtrsim 1$ TeV, neutrino and antineutrino detection to reveal crucial
characteristics of an evaporating primordial black hole at the time of its
burst when observed near Earth. We improve upon previous calculations by
carefully accounting for the non-isotropic particle emission, as Earth
occupies a privileged angle relative to the black hole’s rotation axis.
Additionally, we investigate the angular dependence of primary and secondary
photon spectra and assess the evaporating black hole’s time evolution during
the final explosive stages of its lifetime. Since photon events outnumber
neutrinos by about three orders of magnitude, we find that a neutrino
measurement can aid in identifying the initial angular momentum and the black
hole hemisphere facing Earth only for evaporating black holes within our solar
system, at distances $\lesssim 10^{-4}$ pc, and observed during the final 100
s of their lifetime. Codes used in this work will be publicly available in .

††preprint: IPPP/23/36

##  I Introduction

The Early Universe’s high density provides an ideal environment for the
formation of black holes (BH) with masses significantly smaller than those
observed in Gravitational Waves or direct astrophysical measurements [1, 2,
3]. These _primordial_ black holes (PBHs) could profoundly impact the
evolution of our Universe [4, 5, 6]. Given their potentially tiny masses, as
small as the Planck mass, it becomes imperative to take into account quantum
effects in the evolution of PBHs. Using a semi-classical approximation,
Hawking demonstrated that BHs evaporate by emitting a flux of particles with a
thermal spectrum [7, 8]. As a consequence of energy conservation arguments,
the BH mass decreases at a rate $\dot{M}\sim M^{-2}$, initially leading to a
slow evaporation. However, particle emission accelerates, triggering a runaway
effect that culminates in an explosive stage. Such evaporation process,
assuming the Standard Model’s (SM) degrees-of-freedom (dofs), enables to
estimate the initial mass of a PBH whose lifetime matches the age of the
Universe, yielding a value of $\sim 10^{15}$ g. As a result, PBHs with masses
smaller than this value would have already evaporated, establishing
constraints to be derived on the initial PBH abundance if they evaporated
during the Big-Bang Nucleosynthesis [9, 10, 11], reionization [12, 13] or the
formation of the CMB [9, 10]. Furthermore, due to the universality of Hawking
evaporation, PBHs could have produced dofs that do not interact with the SM
sector, potentially contributing to the observed Dark Matter (DM) [14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32], to the
relativistic dofs, parameterized by $\Delta N_{\rm eff}$ [33, 16, 34, 17, 35,
27, 33], or modify the generation of the matter-antimatter asymmetry [36, 37,
38, 18, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50]. Conversely, if their
initial mass exceeds $10^{15}$ g, these PBHs could potentially contribute to a
portion of the observed DM [51, 10, 52, 53, 54]. For a comprehensive review,
see Ref. [55].

There remains the possibility that some of the PBHs would have an initial mass
such that they would be evaporating today. These evaporating primordial black
holes (EPBHs) would emit a large flux of particles during the final explosive
stages of their lifetimes, presenting a potential opportunity for observation
in different facilities. Thus, experiments designed to search for gamma-ray
bursts, such as H.E.S.S [56, 57], Milagro [58], VERITAS [59], HAWC [60, 61],
and Fermi-LAT [62], have already placed constraints on the existence of PBHs
in our local region. However, it is worth noting that other particles, such as
neutrinos [63, 64, 65], could also be measurable if an EPBH happens to occur
close enough to Earth. Undoubtedly, observing PBH evaporation in its final
stages would not only be a triumph for theoretical physics but also hold the
potential to unlock a wealth of new physics. For instance, this unique
phenomenon could provide invaluable insights into the complete spectrum of
particles, including those beyond-the-Standard-Model (BSM), that exist nature
[66, 67, 68, 69].

Another avenue to explore BSM physics in EPBHs involves determining whether
these objects had a substantial angular momentum at the onset of their final
bursts. In a scenario that includes only the SM set of particles, even a
close-to-maximally rotating PBH should have lost most of its angular momentum
in the early stages of the Hawking evaporation [70, 71]. However, the presence
of an abundant number of scalar dofs, as expected in string axiverse models
[72], could significantly modify this behaviour. In such cases, the BH’s
angular momentum might not deplete as rapidly as in the standard scenario
since scalar particles only reduce the BH mass but not its angular momentum
[73, 74]. This leads to PBHs approaching their final stages with substantial
angular momenta [75, 76]. Thus, investigating the angular momentum properties
of EPBHs offers an intriguing window into potential BSM physics.

Fig. 1:  An illustration of the Earth-Evaporating Kerr Black Hole system
considered in this work. Earth is placed at an angle $\theta$ with respect to
the rotation axis and at a distance $d_{L}$. Earthrise picture reproduced from
[77].

Previous studies have examined the impact of non-zero angular momentum on the
photon and neutrino emissions from an EPBH [65, 76]. These works used the
BlackHawk code [78, 79], which calculates both the primary particle spectrum
directly emitted during evaporation and the secondary spectrum originating
from the decay of unstable particles produced by the black hole. However, it
is important to note that BlackHawk computes the angle-integrated spectrum,
which is not suitable for a Kerr EPBH since rotating EPBHs exhibit a non-
isotropic particle emission due to their axisymmetric nature. Consequently, if
an EPBH were to burst close to Earth, observatories would measure the particle
flux at a unique angle with respect to the rotation axis, leading to distinct
observational characteristics, see Fig. 1 for an artistic depiction of the
Earth–EPBH system considered in this work. Additionally, in the specific case
of neutrinos, it has been observed that due to the parity violation in weak
interactions, their emission from an EPBH follows an asymmetric pattern:
neutrinos are preferentially emitted in the EPBH’s “southern” hemisphere,
while antineutrinos are predominantly produced in the “northern” hemisphere
[80, 81, 82, 83, 84]. This well-established behavior opens up new avenues for
exploring the angular momentum of an EPBH since neutrino emission is
intrinsically linked to the black hole’s spin and varies with the polar angle
relative to its rotation axis.

In this work, we conduct a comprehensive analysis of the angular distribution
of neutrino and photon emissions, both primary and secondary, from Kerr EPBHs.
Additionally, we investigate the potential of gamma-ray and neutrino
observatories to discern crucial EPBH properties, such as its angular momentum
and the hemisphere facing Earth. Despite the photon events outnumbering
neutrino events by approximately three orders of magnitude, the measurement of
neutrinos is instrumental in determining the EPBH hemisphere facing Earth,
particularly in the scenario where such an object is in close proximity to our
planet.

The structure of this paper is as follows. In Sec. II, we provide a detailed
description of the neutrino asymmetric emission and the photon angular
distribution. Furthermore, within the same section, we first outline a
numerical procedure to compute the angular dependence of the secondary
particle emission. Moving forward, in Sec. III, we delve into the time
evolution of a Kerr EPBH during its final burst, critically assessing its
impact on the neutrino emission asymmetry. Sec. IV is exclusively dedicated to
exploring the detection prospects, with specific emphasis on the IceCube [85,
86] and HAWC [61] observatories. Moreover, we thoroughly examine the
feasibility of determining the initial characteristics of EPBHs through the
combination of neutrino and photon measurements. In Sec. V, we draw our
conclusions. We have included two appendices: App. A provides a short review
on the Dirac equation in the Kerr spacetime, and App. B contains the effective
potentials for scalar and vectors used to obtain absorption probabilities. 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, throughout this
manuscript.

##  II Angular distribution of Particle Emission from Kerr Black Holes

Let us consider the Hawking radiation emitted from a Kerr BH with an
instantaneous mass $M$ and a dimensionless spin parameter $a_{\star}\equiv{\rm
J}/(GM)^{2}\in[0,1)$, where ${\rm J}$ represents the BH angular momentum. To
analyse the angular distribution of particles, we adopt the Boyer-Lindquist
coordinate system $(t,r,\theta,\phi)$. In these coordinates, the line element
for the Kerr spacetime is given by [87]

$\displaystyle\differential s^{2}=\frac{\Delta}{\Sigma}(\differential
t-a\sin^{2}\theta\differential\phi)^{2}-\frac{\sin^{2}\theta}{\Sigma}(-a\differential
t+(r^{2}+a^{2})\differential\phi)^{2}-\frac{\Sigma}{\Delta}\differential
r^{2}-\Sigma\differential\theta^{2},$ (1)

where $\Delta\equiv r^{2}-2GMr+a^{2}$, $\Sigma=r^{2}+a^{2}\cos\theta$, with
$a=a_{\star}GM$. It has been demonstrated that the equations of motion for
bosons and fermions are separable, leading to the well-known Teukolsky master
equations [88, 89, 90]. Such a separability will allow us to analyse the
angular distribution of particles, and thus the possible determination of the
EPBH initial spin at the onset of a burst. Next, we will consider the primary
neutrino emission and its dependence on the polar angle $\theta$, as well as
the photon emission. Subsequently, we will examine the secondary spectra
resulting from the decay of unstable particles produced by the BH evaporation.

###  II.1 Primary Neutrino Emission and Asymmetry

Early studies focused on neutrinos as a prototype for understanding the
properties of massless fermions in the Kerr spacetime [80]. This choice was
influenced by the prevailing belief at that time that neutrinos were massless.
However, we now know that neutrinos are indeed massive particles that exhibit
mixing. Taking this into account, we have recently investigated the
implications of neutrino mass on the process of Hawking evaporation for
Schwarzschild BHs [91]. Here, we extend the discussion of massive neutrinos to
Kerr BHs, and examine their asymmetric emission as a function of the polar
angle.

The Dirac equation and its separability for massive fermions in the Kerr
background has been extensively studied in Refs. [92, 93, 94, 95]. For sake of
completeness, we briefly describe the properties of this equation in App. A.
After proposing an ansatz, we obtain equations for radial $R_{1,2}(r)$ and
angular functions $S_{1,2}(\theta)$ appearing in the spinor components,
depending on the energy $\omega$, and $l,m$ the total and axial angular
momentum quantum numbers [95] (for details see App. A)

$\displaystyle\sqrt{\Delta}(\partial_{r}-iK/\Delta)R_{1}$
$\displaystyle=(\\!\,_{\frac{1}{2}}A_{lm}+i\mu r)R_{2},$ (2a)
$\displaystyle\sqrt{\Delta}(\partial_{r}+iK/\Delta)R_{2}$
$\displaystyle=(\\!\,_{\frac{1}{2}}A_{lm}-i\mu r)R_{1},$ (2b)

and

$\displaystyle\left(\partial_{\theta}-\frac{1}{2}\cot\theta-m\csc\theta+a\omega\sin\theta\right)S_{1}(\theta)$
$\displaystyle=(+\\!\,_{\frac{1}{2}}A_{lm}+a\mu\cos\theta)S_{2}(\theta),$ (3a)
$\displaystyle\left(\partial_{\theta}-\frac{1}{2}\cot\theta+m\csc\theta-a\omega\sin\theta\right)S_{2}(\theta)$
$\displaystyle=(-\\!\,_{\frac{1}{2}}A_{lm}+a\mu\cos\theta)S_{1}(\theta),$ (3b)

where $K\equiv(r^{2}+a^{2})\omega-am$ and $\mu$ is the fermion mass. The
quantities $\\!\,{}_{\frac{1}{2}}A_{lm}$ represent the angular separation
constants, which are simply $\\!\,{}_{\frac{1}{2}}A_{lm}={\cal P}(l+1/2)$,
${\cal P}=\pm 1$, for a Schwarzschild BH.

The solutions to the radial equations, Eqs. (2), allow us to determine the
absorption probabilities $\\!\,{}_{\frac{1}{2}}\Gamma_{lm}$, which are crucial
for calculating the Hawking emission rate as they describe the effects of
centrifugal and gravitational potentials on particle production [7, 8, 70, 71,
96]. Numerical methods are required to find these solutions. We employ the
general procedure for Kerr-Newman BHs established in Ref. [96], where a
transformation of variables is applied to simplify the numerical treatment.
Once the solutions have been obtained, we calculate the absorption
probabilities $\\!\,{}_{\frac{1}{2}}\Gamma_{lm}$ using the expressions
provided in the appendix of Ref. [96].

To investigate the dependence of neutrino production on the polar angle
$\theta$ from an EPBH, it is also necessary to solve the angular equations,
Eqs. (3). Several numerical methods have been proposed in the literature to
solve these equations, see e.g. [97, 98, 94, 99]. In this work, we closely
follow the approach described in Ref. [94]. The method involves employing a
series expansion for the eigenfunctions $\\!\,{}_{\pm\frac{1}{2}}S_{lm}\equiv
S_{1,2}$ using associated Legendre polynomials, leading to a continued
fraction equation. By numerically solving this equation, we obtain the angular
eigenvalues $\\!\,{}_{\frac{1}{2}}A_{lm}$. Subsequently, the coefficients of
the series expansion for the eigenfunctions are computed using these angular
eigenvalues. Through this numerical procedure, we find the angular
eigenfunctions up to a normalization factor. To fix such a factor, we impose
the normalization condition,

$\displaystyle\int|\\!\,_{\pm\frac{1}{2}}S_{lm}(\theta)|^{2}\,\differential\Omega=1.$
(4)

The rates of neutrino emission as a function of time, energy, and solid angle
$\Omega$, are computed after considering the quantization of fermions in a
Kerr background. The quantum field theory for massless fermions has been
extensively studied in the literature, and interested readers can find
detailed discussions in Refs. [82, 80, 100, 101, 102]. However, to the best of
our knowledge, currently there is no available treatment for the quantization
of massive fermions in the Kerr geometry; for potential issues related to this
particular case, see Ref. [101]. Nevertheless, it is worth noting that the
neutrino mass $m_{\nu}\lesssim 1~{}\text{eV}$ is negligible compared to the
energies observed in neutrino telescopes, $E_{\nu}\gtrsim\text{100 GeV}$.
Consequently, we can reasonably _assume_ that the quantization procedure used
for massless fermions is applicable with a high level of accuracy to neutrinos
originating from an EPBH, a supposition that we will adopt from now on.

The main physical quantity of interest is the expectation value of the current
operator111Note that the commutator in the current definition should be taken
only with respect to the quantum operators and not the spinor structure.

$\displaystyle J^{\mu}=\frac{1}{2}[\overline{\nu},\gamma^{\mu}P_{L}\nu]$ (5)

where $\nu$ represents the neutrino field. Here,
$P_{L}=\frac{1}{2}(1-\gamma^{5})$ is the left chiral projector defined in
terms of the chirality matrix $\gamma^{5}$, see App. A. As discussed in Ref.
[91], neutrinos and antineutrinos are expected to be emitted as mass
eigenstates, $\nu_{1,2,3}$. The neutrino field $\nu$ will represent generally
any of the three fields associated to the mass eigenstates henceforth. The
radial component of the current $J^{r}$ describes the flow of neutrinos minus
antineutrinos [82, 101]. The net neutrino minus antineutrino flux, ${\cal
A}\equiv N_{\nu}-N_{\overline{\nu}}$, per unit time and solid angle in the
radial direction _away_ from the BH corresponds to the expectation value [82]

$\displaystyle\frac{\differential^{2}{\cal
A}}{r^{2}\differential\Omega\differential t}=\langle
U^{-}|J^{r}|U^{-}\rangle,$ (6)

where $|U^{-}\rangle$ represents the past Unruh vacuum [103]. This vacuum
state, defined as the state that does not contain any particles incoming from
the null past infinity ${\cal J}^{-}$, is the relevant for our purposes since
we are considering the evaporation of a BH originated from a gravitational
collapse event.

In terms of the angular functions and absorption probabilities defined
previously, the net neutrino emission flux can be expressed as follows [80,
81, 82, 83, 101]

$\displaystyle\frac{{\differential^{2}{\cal
A}}}{\differential\Omega\differential
t}=\frac{1}{4\pi}\sum_{l=1/2}\sum_{m=-l}^{l}\int_{0}^{\infty}\differential\omega\frac{\\!\,{}_{\frac{1}{2}}\Gamma_{lm}}{\exp(\varpi/T)+1}\\{|\\!\,_{-\frac{1}{2}}S_{lm}(\theta)|^{2}-|\\!\,_{+\frac{1}{2}}S_{lm}(\theta)|^{2}\\},$
(7)

where $\varpi=\omega-m\vartheta$ with
$\vartheta=a_{\star}/(2GM(1+\sqrt{1-a_{\star}^{2}}))$ representing the
horizon’s angular velocity. From the form of the emission asymmetry, we can
interpret the term proportional to $\\!\,{}_{-\frac{1}{2}}S_{lm}(\theta)$ as
the contribution coming from neutrinos, while
$\\!\,{}_{-\frac{1}{2}}S_{lm}(\theta)$ corresponds to the antineutrino term,
with the caveat that such interpretation is valid only far from the EPBH
[100]. In Eq. (7), it is necessary to write explicitly the sum over the
angular momentum quantum numbers as the PBH spin breaks the spherical
symmetry, given the explicit dependence of the Hawking rate on $m$. The
Hawking temperature $T$ is given by

$\displaystyle T=\frac{1}{4\pi
GM}\,\frac{\sqrt{1-a_{\star}^{2}}}{1+\sqrt{1-a_{\star}^{2}}}.$ (8)

Thus, we can define the net particle emission rate within time $\differential
t$, solid angle $\differential\Omega$ and energy
$[\omega,\omega+\differential\omega]$ as

$\displaystyle\frac{{\differential^{3}{\cal
A}}}{\differential\omega\differential
t\differential\Omega}=\frac{1}{4\pi}\sum_{l=1/2}\sum_{m=-l}^{l}\frac{\\!\,{}_{\frac{1}{2}}\Gamma_{lm}}{\exp(\varpi/T)+1}\\{|_{-\frac{1}{2}}S_{lm}(\theta)|^{2}-|_{+\frac{1}{2}}S_{lm}(\theta)|^{2}\\}.$
(9)

Notice that the integration of the rate over the solid angle produces a
vanishing net particle emission, as expected. Analogously, let us define the
total emission rate, i.e., neutrino plus antineutrino rate as

$\displaystyle\frac{\differential^{3}N_{\nu+\overline{\nu}}}{\differential\omega\differential
t\differential\Omega}=\frac{1}{4\pi}\sum_{l=1/2}\sum_{m=-l}^{l}\frac{\\!\,{}_{\frac{1}{2}}\Gamma_{lm}}{\exp(\varpi/T)+1}\\{|_{-\frac{1}{2}}S_{lm}(\theta)|^{2}+|_{+\frac{1}{2}}S_{lm}(\theta)|^{2}\\}.$
(10)

This total flux of neutrinos comes directly from the evaporation of the BH,
thus constituting what is referred to as the primary spectrum in the
literature.

Fig. 2:  (Top left) Instantaneous neutrino minus antineutrino net particle
flux, ${\cal A}=N_{\nu}-N_{\overline{\nu}}$, as function of energy for a BH
having a mass of $M=10^{9}~{}{\rm g}$ and spin parameter $a_{\star}=0.999$ for
different values of the polar angle $\theta\in[0,\pi]$. (Bottom left)
Instantaneous asymmetry ratio ${\cal R}$ between the net neutrino minus
antineutrino flux to the total emission rate
$\differential^{3}N_{\nu+\overline{\nu}}/\differential\omega\differential
t\differential\Omega$, cf. Eq. (10). (Right) Energy-integrated instantaneous
neutrino emission asymmetry as function of the polar angle $\theta$ for
different values of the spin parameter, $a_{\star}=0.999$ (purple full),
$a_{\star}=0.8$ (dashed blue), $a_{\star}=0.5$ (dotted dark yellow), and
$a_{\star}=0.1$ (dot-dashed red).

On the top left of Fig. 2, we present the instantaneous neutrino minus
antineutrino net flux $\differential^{3}{\cal
A}/\differential\omega\differential\Omega\differential t$ as function of the
neutrino energy for a BH having a mass of $M=10^{9}~{}{\rm g}$ and spin
parameter $a_{\star}=0.999$. The polar angle at which the observer is placed
varies from $0^{\circ}$ to $180^{\circ}$. When $\theta=0$, we observe that the
net rate is negative and presents a single peak. This distinctive feature
arises because only the angular eigenfunctions with $l=m=\frac{1}{2}$
contribute at the poles, while the remaining modes vanish there. Moreover, the
negative sign indicates that more antineutrinos than neutrinos are emitted. As
we increase the polar angle, the contribution from higher $l$ modes becomes
significant, resulting in additional peaks in the net flux. The maximum
emission asymmetry occurs at approximately $\theta\approx 60^{\circ}$, which
aligns with previous findings [82]. For polar angles within the range
$60^{\circ}\lesssim\theta<180^{\circ}$, the net flux is reduced, and at
$\theta=\pi$, the net flux completely vanishes due to a cancellation between
the angular eigenfunctions $\\!\,{}_{\pm\frac{1}{2}}S_{lm}(\theta)$ at the
equator. In the southern hemisphere, $90^{\circ}\leq\theta\leq 180^{\circ}$,
the net flux becomes positive, indicating a higher number of emitted neutrinos
compared to antineutrinos. Despite this sign change, the behavior in the
southern hemisphere is analogous to that in the northern hemisphere,
$0^{\circ}\leq\theta\leq 90^{\circ}$,. The net flux reaches its maximum value
around $\theta\approx 120^{\circ}$, while at $\theta=180^{\circ}$, only the
$l=m=\frac{1}{2}$ mode contributes.

To quantify the EPBH emission asymmetry with respect to the total
instantaneous neutrino flux, we introduce the asymmetry ratio ${\cal R}$ which
compares the net neutrino minus antineutrino flux to the total emission rate,

$\displaystyle{\cal R}\equiv\frac{\frac{{\differential^{3}{\cal
A}}}{\differential\omega\differential
t\differential\Omega}}{\frac{\differential^{3}N_{\nu+\overline{\nu}}}{\differential\omega\differential
t\differential\Omega}}.$ (11)

We present on the bottom left of Fig. 2 the asymmetry ratio ${\cal R}$ for the
same values for the polar angles $\theta$ as in the top pannel. The results
indicate that at the poles ($\theta=0^{\circ},180^{\circ}$), the asymmetry
ratio approaches -1 and 1, respectively. This implies that at these angles,
the flux is predominantly composed of antineutrinos and neutrinos. For angles
different from the poles, the absolute value of the ratio $|{\cal R}|$ remains
below 1, and it is constant for neutrino energies $E_{\nu}\lesssim 5\times
10^{4}$ GeV, where the dominating mode is $l=m=\frac{1}{2}$. However, at
higher energies, the contribution of higher-$l$ modes becomes evident,
resulting in ripples in the ratio. For the specific case of
$\theta=60^{\circ}$, the asymmetry ratio ${\cal R}\approx-0.5$ for energies
$E_{\nu}\lesssim 5\times 10^{4}$ GeV. This confirms that at this angle, there
are approximately $50\%$ more antineutrinos than neutrinos being emitted.
Conversely, at $\theta=120^{\circ}$, we observe the opposite behavior, as
expected from the previous discussion. Finally, at the equator,
$\theta=90^{\circ}$, the ratio ${\cal R}$ is precisely 0, indicating that the
number of emitted neutrinos and antineutrinos is equal at this angle.

To comprehend the dependence of the neutrino-antineutrino net flux on the BH
spin parameter, we present in the right panel of Fig. 2 the energy-integrated
emission asymmetry, $\differential^{2}{\cal
A}/\differential\Omega\differential t$, cf. Eq.(7), for varying spin parameter
values, $a_{\star}=0.999$ (purple full), $a_{\star}=0.8$ (dashed blue),
$a_{\star}=0.5$ (dotted dark yellow), and $a_{\star}=0.1$ (dot-dashed red).
For a close-to-maximally rotating BH with $a_{\star}=0.999$, we observe that
the net emission is maximal at a value $|\cos\theta|\approx 0.5$, whereas for
other spin parameter values, the emission asymmetry peaks at the poles. This
behavior is a consequence of higher-$l$ modes increasingly dominating for
$a_{\star}\gtrsim 0.8$, which in turn contribute away from the poles. On the
other hand, for $a_{\star}\lesssim 0.8$, the $l$ mode that predominantly
dominates is the $l=m=\frac{1}{2}$, resulting in a maximum emission at the
poles. As a result, the net emission asymmetry peaks at the poles for these
values of the spin parameter. For all spin parameter values presented, we
observe that the net flux is negative for $\cos\theta>0$, indicating a
preference for antineutrino emission over neutrinos in the northern
hemisphere. Conversely, for negative values of $\cos\theta$, the behavior is
opposite, with a preference for neutrino emission. At the equator, the net
emission asymmetry vanishes, ensuring that the numbers of emitted neutrinos
and antineutrinos coincide. This overall pattern demonstrates that
antineutrinos are preferentially emitted in the northern hemisphere, whereas
neutrinos are predominantly produced in the southern hemisphere.

Neutrino physics holds an unresolved question regarding the fundamental
fermionic nature of these particles. Neutrinos can be either their own
antiparticles, known as Majorana fermions, or as distinct entities from their
antiparticles, referred to as Dirac neutrinos. In Ref. [91], we investigated
the impact of each scenario on Hawking radiation from Schwarzschild black
holes. Now, the question arises as to the main effect for Kerr black holes. In
the case of Dirac neutrinos, additional dofs come into play, right-handed
neutrinos and left-handed antineutrinos. Consequently, the black hole would
emit these states asymmetrically, with right-handed neutrinos predominantly
emitted in the northern hemisphere and left-handed antineutrinos in the
southern hemisphere. This raises the possibility that the previously mentioned
overall asymmetry might be removed. However, the detection of these additional
states via weak interactions is hindered by a helicity factor of
$m_{\nu}/E_{\nu}\sim 10^{-11}$, for the energy range of interest. As a result,
the presence of the additional states becomes negligible, and the asymmetry
remains intact.

For Majorana neutrinos, only two dofs exist per active mass eigenstate,
corresponding to positive and negative helicities. While there is technically
no distinct state known as an antineutrino, weak interactions differentiate
between the helicity states due to parity violation. Consequently, a positive-
helicity state exhibits a different interaction compared to a negative-
helicity one in the ultra-relativistic limit, which is of interest to us here.
Thus, in such a limit, it is customary to denote the negative-helicity states
as “neutrinos” and the positive-helicity ones as “antineutrinos” [104, 105].
As a result, the asymmetric emission is also present for Majorana neutrinos.
Consequently, detecting a discrepancy between the numbers of neutrinos and
antineutrinos in a future detection of an EPBH could indicate the presence of
a nonzero angular momentum.

###  II.2 Primary Photon Emission

The emission of particles with higher spin, such as photons, is enhanced for
rotating BHs [70]. Consequently, we anticipate that photon emission will be
strongly influenced by both the EPBH spin and the polar angle $\theta$. To
investigate this dependence, we consider the equations for the radial
$\psi_{s}$ and the angular functions $S(\theta)$ for a massless field with
spin $s$,

$\displaystyle\frac{\differential^{2}\psi_{s}}{\differential
r_{\star}^{2}}+(\omega^{2}-V_{s}(r)\psi_{s}$ $\displaystyle=0,$ (12a)
$\displaystyle\frac{1}{\sin\theta}\frac{\differential}{\differential\theta}\left(\sin\theta\frac{\differential
S}{\differential\theta}\right)+[(c\cos\theta-s)^{2}-(m\csc\theta+s\cot\theta)^{2}-s(s-1)+\\!\,_{s}A_{lm}]\,S$
$\displaystyle=0,$ (12b)

where $\lambda_{s}\equiv\\!\,_{s}A_{lm}+c^{2}-2mc$, being $c=a\omega$, and
$\\!\,{}_{s}A_{lm}$ are the angular eigenvalues, and $r_{\star}$ is the
Eddington–Finkelstein radial coordinate

$\displaystyle\frac{\differential r_{\star}}{\differential
r}=\frac{r^{2}+a^{2}+am/\omega}{\Delta}.$ (13)

The Schrödinger-like equation Eq. (12a) is obtained through the Chandrashekar-
Detweiler method as presented in Refs. [106, 107, 108]. The potentials $V_{s}$
depend on the spin of the particle, which, for completeness, are provided in
App. B. To compute the absorption probabilities $\\!\,{}_{s}\Gamma_{lm}$ for
massless bosons, we solve the wave equation, Eq. (12a), and compute the
transmission coefficient for a purely ingoing wave222It is worth noting that
we have verified the consistency of our results for the absorption
probabilities $\\!\,{}_{s}\Gamma_{lm}$ with those obtained from the BlackHawk
code [78, 79].. On the other hand, for the angular eigenvalues and
eigenfunctions, we adopt the approach outlined in Refs. [109, 110], where the
angular separation constants are obtained by numerically solving a continued
fraction equation, similar to the case of massive fermions.

Similarly to the neutrino case, we define the total photon emission rate as
function of energy, time and solid angle as [111, 112]

$\displaystyle\frac{\differential^{3}N_{\gamma}}{\differential\omega\differential
t\differential\Omega}=\frac{1}{4\pi}\sum_{l=s}\sum_{m=-l}^{l}$
$\displaystyle\frac{\\!\,{}_{1}\Gamma_{lm}}{\exp(\varpi/T)-1}\times$
$\displaystyle\\{|\\!\,_{-1}S_{lm}(\theta)|^{2}+|\\!\,_{+1}S_{lm}(\theta)|^{2}\\}\,.$
(14)

Fig. 3 illustrates the instantaneous photon emission rate for various polar
angles of an EPBH with a mass of $M=10^{9}$ g and a spin of $a_{\star}=0.999$.
The overall emission rate behavior is similar to that of neutrino-antineutrino
emission asymmetry. At a polar angle of $\theta=0$ (black curve), only the
$l=1$ mode contributes to the emission. As the polar angle increases, the
contribution from other angular modes becomes noticeable, while the
contribution from the $l=1$ mode diminishes. An important observation is that
the emission rate exhibits a symmetry under the transformation
$\theta\to\pi-\theta$, as the angular modes satisfy [110, 112],

${}_{-s}S_{lm}(\theta)=(-1)^{-l-m}\\!\,_{s}S_{lm}(\pi-\theta).$

Consequently, this symmetry implies that the measurement of photons emitted by
an EPBH cannot determine, in principle, the hemisphere pointing towards Earth.
Upon careful examination of the total emission rate, Eq. (II.2), we observe
that the two different angular eigenfunctions possess distinct polarisations.
Therefore, in principle, measuring the polarisation of the gamma rays could
determine the value of $\theta$. Polarisation is commonly measured through the
Compton scattering angle of photons for energies below electron-positron pair
production, $E_{\gamma}\lesssim 1$ MeV [113]. For energies up to
$E_{\gamma}\lesssim 10$ MeV, the event distribution of electron-positron pairs
can be analysed [114, 115, 116]. However, at higher energies, these techniques
face limitations due to factors like multiple Coulomb scatterings. Some works
have proposed ideas for measuring polarisation up to energies of
$E_{\gamma}\sim$ 30 GeV [114, 115, 116, 113, 117]. Nonetheless, it remains
uncertain if these techniques are applicable to energies of interest here,
$E_{\gamma}\gtrsim 1$ TeV. Furthermore, Earth’s magnetic field may influence
the polarisation of incoming photons from the EPBH. Thus, in this study, we
adopt a conservative approach and assume that future experiments will not
measure gamma-ray polarisation. To resolve the ambiguity in measuring
$\theta$, we propose a multimessenger approach that incorporates the detection
of neutrinos.

Fig. 3:  Primary instantaneous photon emission spectrum from a BH having a
mass $M=10^{9}$ g and spin parameter $a_{\star}=0.999$ for different values of
the polar angle, $\theta=0.^{\circ}$ (black), $15^{\circ}$ (beige),
$30^{\circ}$ (light orange), $60^{\circ}$ (fucsia), $75^{\circ}$ (purple), and
$90^{\circ}$ (blue).

###  II.3 Secondary spectra

Since the EPBH also emits other SM dofs, most of which are unstable, there is
an additional production of neutrinos/photons arising from their subsequent
decay. This additional contribution is known as the neutrino/photon secondary
spectrum. To obtain this secondary spectrum, we convolve the primary spectrum
of each particle species, denoted by $i$, with the number of daughter
neutrinos/photons resulting from their decay as a function of energy and solid
angle. The expression for the secondary spectrum is given by

$\displaystyle\frac{\differential^{3}N_{\nu(\gamma)}^{\rm
sec}}{\differential\omega\differential
t\differential\Omega}=\int_{0}^{\infty}\differential\omega^{\prime}\int\differential\Omega^{\prime}\sum_{i}\frac{\differential^{3}N_{i}}{\differential\omega^{\prime}\differential
t\differential\Omega^{\prime}}\frac{\differential^{2}n_{i\to\nu(\gamma)}}{\differential\omega\differential\Omega}(\omega,\omega^{\prime}).$
(15)

where
$\differential^{2}n_{i\to\nu(\gamma)}/\differential\omega\differential\Omega$
denotes the energy and angular distribution of neutrinos (photons) resulting
from the decay of the $i$-th particle, and
$d^{3}N_{i}/d\omega^{\prime}dtd\Omega^{\prime}$ represents the primary
emission rate of the $i$-th particle species defined as

$\displaystyle\frac{\differential^{3}N_{i}}{\differential\omega\differential
t\differential\Omega}=\frac{g_{i}}{4\pi}\sum_{l=s_{i}}\sum_{m=-l}^{l}$
$\displaystyle\frac{\\!\,{}_{s_{i}}\Gamma_{lm}}{\exp(\varpi/T)-(-1)^{2s_{i}}}\times$
$\displaystyle\\{|\\!\,_{s_{i}}S_{lm}(\theta)|^{2}+|\\!\,_{-s_{i}}S_{lm}(\theta)|^{2}\\},$
(16)

being $s_{i}$ the spin and $g_{i}$ the internal dofs of the particle species
$i$. To compute the secondary spectrum, we thus need to determine the angular
and energy distributions of neutrinos resulting from the decay of each SM
particle. In our analysis, we have employed a similar approach to that used by
BlackHawk [78, 35] which uses the PYTHIA event generator [118] to calculate
these distributions. We use PYTHIA since it has the full information of the
four momenta of the produced particles, which is crucial for our purposes.

The numerical strategy employed to determine
$\differential^{2}n_{i\to\nu(\gamma)}/\differential\omega\differential\Omega$
can be summarized as follows. First, we set the center-of-mass energy of the
collision in PYTHIA to be twice the energy of the primary particle
$\omega^{\prime}$. We generate a large number of events, typically around
$10^{5}$, for a specific channel, $e^{+}+e^{-}\rightarrow
i+\overline{i}\rightarrow\cdots$, where $i$ represents the primary particle of
interest. For each event, we record the energy of the primary particle and
calculate the relative angle between each daughter particle and the primary
particle. To facilitate the analysis, we perform a three-dimensional rotation
that aligns the z-axis with the direction of the primary particle. This
rotation is then applied to all the final state particles, ensuring consistent
orientation. Next, we construct a two-dimensional histogram to capture the
distribution of energy and relative angle for each value of the center-of-mass
energy. Each entry in the histogram represents the frequency of occurrence for
a particular combination of energy and angle, normalised by the total number
of events. We repeat this process by systematically varying the center-of-mass
energy in the range of interest, corresponding to the range of $1~{}{\rm
GeV}-10^{6}~{}{\rm GeV}$. Finally, to convert the dependence on relative angle
to a dependence on the solid angle of the particle distributions, we perform a
rebinning of the histograms. This rebinning takes into account the different
possible orientations of the polar and azimuthal angles of the decay products
for a given angle $\theta$ of the primary particle.

Fig. 4:  Total neutrino (full) and antineutrino (dot-dashed) instantaneous
emission rates for three different values of the polar angle,
$\theta=15^{\circ}$ (purple), $\theta=90^{\circ}$ (red), and
$\theta=130^{\circ}$ (green) and for a BH having a mass of $M=10^{9}~{}{\rm
g}$ and spin parameter $a_{\star}=0.999$. Note that neutrino and antineutrino
lines for the secondary spectra and for the primary spectra for
$\theta=90^{\circ}$ lie on top of each other.

In Fig. 4, we present the total—summed over all neutrino states—primary and
secondary neutrino (full) and antineutrino (dot-dashed) emission rates
$d^{3}N_{\nu(\overline{\nu})}/d\omega d\Omega dt$, obtained by adding and
subtracting Eqs. (9) and (10), for $\theta=15^{\circ}$ (purple),
$\theta=90^{\circ}$ (red), and $\theta=120^{\circ}$ (green). The primary
spectra presents the peaks related to the $l$ modes as expected from the
neutrino-antineutrino emission rate in Fig. 2. We observe that for the polar
angle $\theta=15^{\circ}$, belonging to the upper hemisphere, the primary
emission rate of antineutrinos is $\sim 2$ orders of magnitude larger than the
one from neutrinos. Meanwhile, the rates at the equator, $\theta=90^{\circ}$
are equal. For $\theta=120^{\circ}$, the behaviour is the opposite, the
neutrino primary emission rate is larger, although only a factor of $\sim
3.5$. Regarding the secondary spectra, there are no significant differences
between the neutrino and antineutrino fluxes. This is expected since these
neutrinos originate from the decay of other particles that do not exhibit the
same emission asymmetry in their production that primary neutrinos present
[84]. The main distinction arises from the dependence on the polar angle.
Since the secondary spectrum arises from integrating the primary spectrum over
energy and angle, with weights determined by the decay distributions, higher
$l$ modes play a more significant role, especially at elevated energies.
Consequently, it becomes evident that the most substantial secondary spectrum
is produced from the equatorial region; meanwhile, for $\theta=15$, there is a
observable reduction in flux compared to other values.

##  III Time evolution

Once a PBH reaches its final stages of life, particle emission intensifies,
culminating in a potential observable burst for future observatories. These
observatories will track the temporal evolution of a PBH, with an initial mass
linked to the observed duration of the burst. Consequently, it becomes
essential to understand the underlying time evolution. To derive the evolution
equations for the PBH mass and spin, we multiply the integrated emission rate
in solid angle,

$\displaystyle\frac{\differential^{2}N_{i}}{\differential\omega\differential
t}$
$\displaystyle=\int\differential\Omega\,\frac{\differential^{3}N_{i}}{\differential\omega\differential
t\differential\Omega}$
$\displaystyle=\frac{g_{i}}{2\pi}\sum_{l=s_{i}}\sum_{m=-l}^{l}\frac{\\!\,{}_{s_{i}}\Gamma_{lm}}{\exp(\varpi/T)-(-1)^{2s_{i}}},$
(17)

by the total energy of a given particle $\omega$ or by the $m$ quantum number,
and then we integrate over the phase space . Defining the evaporation
functions for mass and angular momentum, $\epsilon_{i}(M,a_{\star})$ and
$\gamma_{i}(M,a_{\star})$ per particle $i$, respectively, as

$\displaystyle\varepsilon_{i}(M,a_{\star})$
$\displaystyle=\frac{g_{i}}{2\pi}\int_{0}^{\infty}\sum_{l=s_{i}}\sum_{m=-l}^{l}\frac{\omega\,_{s_{i}}\Gamma_{lm}}{\exp(\varpi/T)-(-1)^{2s_{i}}}\,d\omega\,,$
(18a) $\displaystyle\gamma_{i}(M,a_{\star})$
$\displaystyle=\frac{g_{i}}{2\pi}\int_{0}^{\infty}\sum_{l=s_{i}}\sum_{m=-l}^{l}\frac{m\,_{s_{i}}\Gamma_{lm}}{\exp(\varpi/T)-(-1)^{2s_{i}}}\,d\omega\,,$
(18b)

and summing over _all_ existing species, we obtain the following system of
coupled equations for the time evolution [70, 71, 14, 33]

$\displaystyle\frac{\differential M}{\differential t}$
$\displaystyle=-\epsilon(M,a_{\star})\frac{M_{p}^{4}}{M^{2}}\,,$ (19a)
$\displaystyle\frac{\differential a_{\star}}{\differential t}$
$\displaystyle=-a_{\star}[\gamma(M,a_{\star})-2\epsilon(M,a_{\star})]\frac{M_{p}^{4}}{M^{3}}\,.$
(19b)

Assuming only the presence of SM dofs, the numerical solution of these
equations reveals that a nearly maximal spinning BH loses its angular momentum
at a rate significantly faster than its mass. However, if there exists a large
sector of scalar particles, such as those in the string axiverse [72], the BH
spin does not completely evaporate but tends towards an asymptotic value [73,
74, 75]. This emphasizes the importance of determining the angular momentum of
the BH prior to its evaporation as a means to constrain these models. In what
follows, and in order to be model-independent, we will assume that the EPBH
follows the time evolution given the SM dofs at the onset of the final burst.
The dependence with BSM scenarios will be considered elsewhere.

Hence, it becomes crucial to investigate whether the neutrino-antineutrino
emission asymmetry could be observed in the final PBH burst after taking into
account the BH time evolution. We integrate the neutrino (antineutrino)
particle flux over time taking into account the time evolution of both mass
and angular momentum of the BH,

$\displaystyle\frac{\differential^{2}N_{\nu(\overline{\nu})}}{\differential\omega\differential\Omega}=\int_{0}^{\tau}\,\differential
t\,\frac{\differential^{3}N_{\nu(\overline{\nu})}}{\differential\omega\differential
t\differential\Omega}(M(t),a_{\star}(t)),$ (20)

where $\tau$ represents the remaining lifetime of the evaporating black hole,
assumed to be equal to the observed burst duration. The parameter $\tau$ also
determines the PBH mass at the onset of the burst, given a specific initial
$a_{\star}^{\rm in}$. For instance, if the burst duration is $\tau=100$ s and
we consider a nearly maximal spin case, $a_{\star}^{\rm in}=0.999$, the
initial mass is $M_{\rm in}\sim 8.3\times 10^{9}$ g. Conversely, for a non-
rotating Schwarzschild BH, $a_{\star}^{\rm in}=0$, the initial mass would be
smaller, approximately $M_{\rm in}\sim 6.3\times 10^{9}$ g, due to the
increased lifetime of a non-rotating BH.

Fig. 5:  Total time integrated neutrino (full) and antineutrino (dot-dashed)
emission rates for two different values of the polar angle,
$\theta=15^{\circ}$ (purple), $\theta=90^{\circ}$ (red), and
$\theta=130^{\circ}$ (green), and for a burst duration of $\tau=100$ s and
initial spin of $a_{\star}^{\rm in}=0.999$.

The total time integrated neutrino and antineutrino fluxes for a burst
duration of $\tau=100$ s and assuming an initial $a_{\star}^{\rm in}=0.999$
are presented in Fig. 5. Consistent with previous figures, we present the
fluxes for polar angles $\theta=15^{\circ}$ (purple), $\theta=90^{\circ}$
(red), and $\theta=120^{\circ}$ (green), distinguishing between the
contributions of the primary and secondary spectra. Despite the spin is
depleted faster than the mass, we observe an asymmetry in the emission for
neutrino and antineutrino fluences depending on the polar angle. For instance,
at $\theta=15^{\circ}$, the antineutrino fluence is greater than the neutrino
fluence by a factor of $\sim 5$ at an energy of $E_{\nu}\sim 8$ TeV. However,
we observe that the peak structure tends to vanish for this value of the polar
angle at higher energies, which in turn will modify the number of events in
neutrino telescopes. For other values of the polar angle the behaviour is
similar to the one encountered for the instantaneous spectra in Fig. 4;
however, the peak structure is less pronounced. Moreover, the secondary
spectra does not present any asymmetry between neutrinos and antineutrinos.
These fluence profiles suggest potential variations in the number of neutrino
vs antineutrino events in neutrino telescopes, which could enable the
determination of the EPBH’s spin during its final stages. This topic is
further explored in the subsequent section.

##  IV Detection Prospects

PBHs could have been formed in the Early Universe due to different mechanisms,
with initial masses $M_{\rm in}$ spanning from $\sim{\cal O}(1)$ g to values
several orders of magnitude larger than the Solar mass [10, 119]. The present-
day distribution of these PBHs depends not only on their formation mechanisms
but also on the subsequent evolution of the Universe, as PBHs may cluster
around galaxies, among other possibilities. Moreover, the PBH number density
is subject to stringent constraints contingent upon their initial masses,
given that they could be undergoing evaporation, producing observable particle
fluxes, such as $\gamma$-rays, today [10, 119]. Let us estimate the number
density of PBHs that could exist in the Solar neighborhood, neglecting the
effect of clustering, i.e., assuming a uniform PBH distribution. Denoting
$\beta^{\prime}$ as the ratio of PBH energy density at formation $\rho_{\rm
PBH}$ to the total energy density $\rho_{\rm tot}$,

$\displaystyle\beta^{\prime}=\alpha^{\frac{1}{2}}\left(\frac{g_{*,in}}{106.75}\right)^{-\frac{1}{4}}\left(\frac{\rho_{\rm
PBH}}{\rho_{\rm tot}}\right),$

with $\alpha$ the gravitational collapse factor and $g_{*,in}$ the
relativistic dofs, we obtain the PBH number density [10]

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

under the assumption of standard cosmological evolution. Consequently, on
average, $\sim 1.5$ PBHs with initial mass of $10^{15}~{}{\rm g}$ are
anticipated to exist within a spherical volume with radius of 1 parsec, for
$\beta^{\prime}=10^{-29}$, consistent with current constraints. Notably, this
number could be higher if clustering indeed influenced the distribution of
PBHs. This implies that the average distance between PBHs is

$\displaystyle d_{\rm sep}\sim
0.88\left(\frac{10^{-29}}{\beta^{\prime}}\right)^{\frac{1}{3}}\left(\frac{M_{\rm
in}}{10^{15}~{}{\rm g}}\right)^{\frac{1}{2}}~{}{\rm pc},$ (22)

If a PBH were to exist in the Solar neighborhood at a distance of
approximately $\sim 10^{-3}$ pc from Earth, the nearest neighboring PBH would
be situated at a distance of roughly $\sim 1$ pc. Consequently, the particle
flux originating from this additional PBH would be reduced by a factor of
$\sim 10^{-6}$, rendering its contribution to the observed flux effectively
negligible. Therefore, let us assume that an EPBH resides in close proximity
to the Solar System, and focus on the possibility of determining the angular
momentum and polar angle orientation of such an EPBH in both neutrino and
photon detectors.

###  IV.1 Neutrino Telescopes

Current and future neutrino telescopes hold the potential to observe the final
stages of PBH evaporation by detecting the flux of high-energy neutrinos
emitted during this process. Let us consider the observation of $\mu$-tracks
in the IceCube observatory resulting from an EPBH burst occurring in close
proximity to Earth. The exceptional angular resolution of IceCube for high-
energy tracks, achieving an impressive precision of $\lesssim 1^{\circ}$ for
energies in the TeV range [85, 86], will play a crucial role in determining
the origin of these events and distinguishing them from background signals.
Furthermore, IceCube has the capability to differentiate between neutrinos and
antineutrinos based on their characteristic inelasticity distributions,
particularly for neutrino energies below $\lesssim 10$ TeV [85]. This energy
range is of particular interest to us.

To accurately estimate the number of $\mu$-track events detected on Earth, we
must account for neutrino flavor oscillations that occur during their journey
from the source to the detector. As mentioned before, the primary neutrinos
and antineutrinos are assumed to be produced as mass eigenstates.
Consequently, the fluence of muon neutrinos (antineutrinos) originating from
the primary spectra at Earth can be expressed as

$\displaystyle F_{\nu_{\mu}(\overline{\nu}_{\mu})}^{\rm
primary}(a_{\star}^{\rm
in},\theta,d_{L})=\frac{1}{d_{L}^{2}}\sum_{i=1,2,3}|U_{\mu
i}|^{2}\frac{d^{2}N_{\nu_{i}(\overline{\nu}_{i})}}{d\omega d\Omega},$ (23)

where $U$ denotes the Pontecorvo-Maki-Nakagawa-Sakata mixing matrix,
$d^{2}N_{\nu_{i}(\overline{\nu}_{i})}/d\omega d\Omega$ represents the time-
integrated fluence of primary neutrinos (antineutrinos) for mass eigenstate
$i$, and $d_{L}$ corresponds to the distance between the PBH and Earth.
Notably, the absence of the usual $4\pi$ factor in the denominator is due to
the lack of spherical symmetry of particle emission from a Kerr EPBH, as
emphasized previously. Since the energies at which these primary neutrinos are
emitted are significantly higher than the neutrino masses, the primary spectra
for all three mass eigenstates are identical.

In contrast, secondary neutrinos are generated through weak interactions,
leading them to be produced in flavor eigenstates. Neutrino oscillations over
the PBH-Earth distances present decoherence since these distances are expected
to greatly exceed the standard oscillation lengths. Consequently, the fluence
for the secondary component can be expressed as

$\displaystyle F_{\nu_{\mu}(\overline{\nu}_{\mu})}^{\rm
secondary}(a_{\star}^{\rm
in},\theta,d_{L})=\frac{1}{d_{L}^{2}}\sum_{i=1}^{3}\sum_{\alpha=e}^{\tau}|U_{\mu
i}|^{2}|U_{\alpha
i}|^{2}\frac{d^{2}N_{\nu_{\alpha}(\overline{\nu}_{\alpha})}^{\rm sec}}{d\omega
d\Omega}.$ (24)

The number of muon neutrino(antineutrino) events in IceCube for a given zenith
angle $\zeta$, corresponding to the location of the EPBH, involves the
following expression

$\displaystyle N_{\nu_{\mu}(\overline{\nu}_{\mu})}(a_{\star}^{\rm
in},\theta,d_{L})=\int_{\omega_{\rm min}}^{\omega_{\rm
max}}F_{\nu_{\mu}(\overline{\nu}_{\mu})}\mathscr{A}_{\rm
eff}(\omega,\zeta)d\omega,$ (25)

where
$F_{\nu_{\mu}(\overline{\nu}_{\mu})}=F_{\nu_{\mu}(\overline{\nu}_{\mu})}^{\rm
primary}+F_{\nu_{\mu}(\overline{\nu}_{\mu})}^{\rm secondary}$ represents the
total neutrino fluence at the detector, and $\mathscr{A}_{\rm
eff}(\omega,\zeta)$ denotes the effective area of IceCube [86]. In principle,
the effective areas should be different for neutrinos and antineutrinos.
However, we use the publicly available effective area [86], which corresponds
to the averaged area for neutrinos and antineutrinos. The energy integration
is performed over the range between IceCube’s threshold energy of $\omega_{\rm
min}=100~{}{\rm GeV}$ and the maximum energy of the neutrino fluence, which,
in turn, depends on the duration of the burst [65]. The main background in
IceCube that could affect the measurement of the neutrinos from an EPBH
corresponds to high-energy atmospheric neutrinos creating observable tracks.
Nevertheless, the events from these neutrinos is of order $10^{-4}$ for a time
interval of $100$ s [65], making such a background negligible.

Figure 6 illustrates the variation of neutrino and antineutrino events in
IceCube as a function of the zenith angle (top) and the ratio of neutrinos to
antineutrinos (bottom) for an EPBH having an initial $a_{\star}^{\rm
in}=0.999$, a burst duration of $\tau=100$ s at distance of $d_{L}=10^{-4}$
pc. Similar to previous figures, we consider the polar angles
$\theta=15^{\circ}$ (purple), $\theta=90^{\circ}$ (red), and
$\theta=120^{\circ}$ (green), representing the orientation of the Earth
relative to the axis of rotation of the PBH. For the case of
$\theta=15^{\circ}$, it is observed that muon antineutrino events dominate,
while neutrino events are only $\sim 80\%$ of those for antineutrinos across
various zenith angles. Conversely, at $\theta=120^{\circ}$, where the Earth is
positioned in the southern hemisphere of the EPBH, neutrino events exceed
antineutrino events by $\sim 12\%$. In line with expectations, when the Earth
aligns with the EPBH’s equatorial plane, neutrino and antineutrino events
coincide.

Fig. 6:  Muon neutrino (full) and antineutrino (dot-dashed) events in IceCube
(top), and the neutrino-to-antineutrino ratio (bottom) as function of the
zenith angle for an EPBH at a distance of $d_{L}=10^{-4}$ pc and a burst
duration of $\tau=100$ s for different values of the polar angle:
$\theta=15^{\circ}$ (purple), $\theta=90^{\circ}$ (red), and
$\theta=120^{\circ}$ (green), representing the orientation of the Earth
relative to the axis of rotation of the PBH.

###  IV.2 Gamma Ray Experiments

Gamma rays, emitted abundantly during the final stages of PBH lifetime, can be
effectively detected by experiments specifically designed for high-energy
photon searches, such as the High Altitude Water Cherenkov observatory (HAWC)
[60, 61]. HAWC is a very-high-energy (VHE) air shower array located on the
slopes of the Sierra Negra volcano at an altitude of 4100 m above sea level.
HAWC’s main array comprises 300 cylindrical water tanks, each equipped with
photomultiplier tubes that detect the Cherenkov light produced by secondary
particles resulting from the interaction of VHE gamma rays. With a wide field
of view of approximately 2 sr, HAWC is capable of detecting photons with
energies ranging from $10^{2}$ to $10^{5}$ GeV [61]. Additionally, HAWC boasts
excellent angular resolution, ranging from approximately $\sim 0.2^{\circ}$ to
$1^{\circ}$, rendering it well-suited for capturing transient events such as
EPBH bursts [61].

To estimate the number of gamma-ray events from the burst of a EPBH, we follow
a similar approach as with neutrinos,

$\displaystyle N_{\gamma}(a_{\star}^{\rm in},\theta,d_{L})=\int_{\omega_{\rm
min}}^{\omega_{\rm max}}F_{\gamma}(a_{\star}^{\rm
in},\theta,d_{L})\mathscr{A}_{\rm eff}^{\rm HAWC}(\omega,\zeta)\,d\omega,$
(26)

where $F_{\gamma}(\theta)$ encompasses both primary and secondary photon
emissions, and $\mathscr{A}_{\rm eff}^{\rm HAWC}(\omega,\zeta)$ represents the
effective area of HAWC [60]. Given the larger cross-section of interaction for
photons, we anticipate that the number of photon events will surpass that of
neutrinos. However, as previously mentioned, the symmetric emission of photons
from both PBH hemispheres introduces an ambiguity in determining the
orientation of the axis of rotation relative to Earth. Consequently, a
simultaneous measurement of both neutrinos and photons would provide valuable
insights, not only regarding the initial angular momentum of the EPBH at the
onset of the burst but also enabling determination of the axis of rotation’s
orientation relative to Earth.

###  IV.3 Sensitivity

To demonstrate the feasibility of a multimessenger approach, we conduct a
sensitivity analysis of combined neutrino/antineutrino and photon
measurements. This analysis aims to understand the potential of such an
approach in determining the initial parameters of the EPBH at the beginning of
the burst. We employ a rate analysis utilizing the following test statistics

$\displaystyle\chi^{2}={\rm
min}_{\alpha}\left\\{\frac{(N_{\nu_{\mu}}(a_{\star}^{\rm
in},\theta,d_{L},\alpha)-N_{\nu_{\mu}}^{\rm b})^{2}}{N_{\nu_{\mu}}^{\rm
b}}+\frac{(N_{\overline{\nu}_{\mu}}(a_{\star}^{\rm
in},\theta,d_{L},\alpha)-N_{\overline{\nu}_{\mu}}^{\rm
b})^{2}}{N_{\overline{\nu}_{\mu}}^{\rm b}}+\frac{(N_{\gamma}(a_{\star}^{\rm
in},\theta,d_{L},\alpha)-N_{\gamma}^{\rm b})^{2}}{N_{\gamma}^{\rm
b}}+\frac{\alpha^{2}}{\sigma_{\alpha}^{2}}\right\\},$ (27)

where $N_{\nu_{\mu}}(a_{\star}^{\rm
in},\theta,d_{L},\alpha),N_{\overline{\nu}_{\mu}}(a_{\star}^{\rm
in},\theta,d_{L},\alpha),N_{\gamma}(a_{\star}^{\rm in},\theta,d_{L},\alpha)$
represent the predicted events of muon neutrinos, muon antineutrinos, and
photons, respectively, for a given set of initial parameters $a_{\star}^{\rm
in},\theta$ for a EPBH at a distance $d_{L}$, measured via some independent
parallax technique. We include a pull parameter $\alpha$ related to the
measurement of the EPBH-Earth distance by performing the substitution

$\displaystyle d_{L}\to(1+\alpha)d_{L},$

in the determination of the number of events, Eqs. (25) and (26). We assume an
uncertainty of $\sigma_{\alpha}=10\%$ on the measurement of the distance
$d_{L}$. $N_{\nu_{\mu}}^{\rm b}$, $N_{\overline{\nu}\mu}^{\rm b}$, and
$N{\gamma}^{\rm b}$ are the corresponding benchmark event values. For our
analysis, we assume benchmark EPBH parameter values of $a_{\star}^{\rm
in}=0.5$, $\theta=45^{\circ}$, and a distance of $d_{L}=10^{-4}$ pc, yielding
values of $N_{\nu_{\mu}}^{\rm b}=800.65$, $N_{\overline{\nu}\mu}^{\rm
b}=873.09$, and $N_{\gamma}^{\rm b}=2.81\times 10^{6}$ for $\tau=100$ s.
Additionally, we assume an ideal scenario with no backgrounds and perfect
discrimination between neutrinos and antineutrinos. While these assumptions
are quite optimistic, they provide a best-case scenario for assessing the
feasibility of determining the properties of an EPBH using neutrinos.

Fig. 7:  Sensitivity to the EPBH’s initial parameters at the onset of a burst
having three durations of $\tau=1$ s (left), $\tau=10$ s (center), and
$\tau=100$ s (right). The allowed regions are at the 95$\%$ CL from photons
(yellow), neutrinos (green), and their combination (purple). We assumed a
distance of $d_{L}=10^{-4}$ pc with a $10\%$ uncertainty, an initial
$a_{\star}^{\rm in}=0.5$ and the Earth placed at a polar angle of
$\theta=45^{\circ}$ with respect to the EPBH axis of rotation.

We perform a two-dimensional $\chi^{2}$ analysis varying the initial
characteristics of the EPBH at the onset of the burst, specifically its
angular momentum and orientation with respect to Earth, $a_{\star}^{\rm in}$
and $\cos\theta$. Moreover, we marginalize over the additional pull parameter
$\alpha$ to obtain the allowed regions. This analysis is performed for three
distinct durations, $\tau=1$ s (left panel), $\tau=10$ s (center panel), and
$\tau=100$ s (right panel). The outcomes of this analysis are visually
represented in Fig. 7. The contributions to the $\Delta\chi^{2}$ from
individual measurements of photons (yellow) and neutrinos (green), as well as
their combination (purple), are presented, with all regions corresponding to a
$95\%$ CL, while the input benchmark value is represented by the cyan star.
Analyzing the photon measurement alone, we observe that the large number of
photon events facilitates a precise determination of the initial EPBH spin at
the onset of the observed burst. However, the allowed regions exhibit a
degeneracy between the hemispheres, resulting in a symmetric region under the
transformation $\cos\theta\to-\cos\theta$, and exhibits a convex-like
structure for the larger burst durations, $\tau=10~{}{\rm s},100~{}{\rm s}$.
This structure arises due to the enhanced emission of photons at the EPBH
poles, resulting in similar gamma-ray event rates for EPBHs with lower spin
parameters in the polar region as compared to higher-spin BHs at different
polar angles pointing towards Earth. For the smaller burst time of $\tau=1$ s,
such a structure is absent since the number of events in the poles differ from
those at the equator at less than $\sim 1\%$, due to the short duration of the
burst.

Given the precision of the photon-only measurement in determining
$a_{\star}^{\rm in}$, one might wonder about the impact of uncertainties in
other parameters, such as the EPBH-Earth distance, $d_{L}$. We tested a range
of distance measurement uncertainties, ranging from $\sigma_{\alpha}=10\%$ to
$100\%$, and observed that even at the upper limit, the allowed regions were
affected by only approximately $10\%$.

In contrast, the neutrino measurement yields a broader region for all burst
durations primarily due to the lower number of neutrino events, which is
approximately $10^{-3}$ smaller than that of photons. However, we observe that
for the shortest burst duration ($\tau=1$ s), neutrinos alone can exclude an
initially close-to-maximally rotating EPBH ($a_{\star}^{\rm in}=0.999$) for
values of $\cos\theta\lesssim 1$ at more than 95$\%$ CL. For the same burst
duration, neutrinos could exclude $a_{\star}^{\rm in}\gtrsim 0.75$ for an EPBH
with its northern pole pointing towards Earth. This is because, in the
considered benchmark scenario, there are more antineutrino events than
neutrino events. Conversely, for $\cos\theta=-1$, the behavior is the
opposite, with neutrino events being approximately 14$\%$ larger than those
for the benchmark, while antineutrinos are reduced by about 5$\%$.

Longer observation times significantly contribute to better determination of
the initial EPBH parameters, as both the number of observed neutrinos and
photons increase with the burst duration. For a burst duration of $\tau=10$ s,
it becomes apparent that the asymmetry in neutrino-antineutrino emission
begins to impact the combined sensitivity, leading to a reduction in the
degeneracy associated with determining the EPBH hemisphere that is oriented
towards Earth. Notably, the previously degenerate solution at
$\theta=135^{\circ}$ is now disfavored with a $\Delta\chi^{2}=5.78$. Moreover,
the neutrino measurement alone excludes values of $a_{\star}^{\rm in}\gtrsim
0.6$ for $\cos\theta\lesssim-0.5$ at a 95 $\%$ confidence level, while for
$\cos\theta\gtrsim 0.5$, the exclusion occurs exclusively for highly rotating
EPBHs, with $a_{\star}^{\rm in}\gtrsim 0.8$ at the same confidence level. Such
a sensitivity is attained through the enhanced emission of fermions for higher
rotating BHs, consequently leading to an increased net particle flux.

The exclusion of the solution in the second quadrant ($\theta=135^{\circ}$) at
level exceeding 95$\%$ is ultimately achieved for the longer burst duration of
$\tau=100$ s. Furthermore, the neutrino measurement independently excludes a
significant portion of the parameter space, particularly for values
$\cos\theta\lesssim-0.1$ and $a_{\star}^{\rm in}\gtrsim 0.1$. For
$a_{\star}^{\rm in}\lesssim 0.1$, we observe that the sensitivity becomes
nearly independent of $\cos\theta$. This is because these EPBHs start to
resemble a Schwarzschild BH, which is spherically symmetric and does not
exhibit neutrino emission asymmetry. In this case, the sensitivity primarily
arises from the difference in events due to the enhanced emission for Kerr
BHs. A similar behavior occurs for $a_{\star}^{\rm in}\gtrsim 0.8$, where the
increased emission of neutrinos and antineutrinos results in a significantly
larger number of events than expected for the chosen benchmark. For instance,
at $a_{\star}^{\rm in}=0.9$ and $\theta=45^{\circ}$, the
$N_{\nu_{\mu}}=840.97$, and $N_{\overline{\nu}\mu}=975.45$, making this value
disfavored by approximately $\Delta\chi^{2}\approx 4\sigma$. Thus, by
combining photon and neutrino measurements, the degeneracy present in the
photon measurement is effectively resolved, leading to a more precise
determination of the EPBH’s initial characteristics.

Fig. 8:  Chi-squared difference between the two degenerate angles $\theta$ and
$\pi-\theta$ solutions in the photon measurement, after combination with the
neutrino measurement, as function of the distance to the EPBH for various
initial spin parameter values, $a_{\star}^{\rm in}=0.999$ (blue),
$a_{\star}^{\rm in}=0.5$ (orange), and $a_{\star}^{\rm in}=0.1$ (green). The
dashed vertical lines indicate the distances of the outer planets to the Sun
and a light-day distance for reference.

It is then natural to question how close the EPBH must be for neutrinos to
effectively break the degeneracy in determining the angle $\theta$. To gain
insights into this matter, we present in Fig. 8 the profiled $\Delta\chi^{2}$
between the two degenerate solutions, $\theta$ and $\pi-\theta$, as a function
of the EPBH-Earth distance for fixed values of the initial spin parameter and
a burst duration of $\tau=100$ s, assuming $\theta=45^{\circ}$. We observe
that, for an initially close-to-maximally rotating black hole, the neutrino
emission asymmetry is capable of breaking the degeneracy at a $3\sigma$
confidence level when the EPBH-Earth distance is approximately $2\times
10^{-4}$ pc. As expected, for lower spinning EPBHs, the required distances are
smaller due to the decreased asymmetry in emission. For $a_{\star}^{\rm
in}=0.5$, the necessary distance is roughly $1.1\times 10^{-4}$ pc, while for
$a_{\star}^{\rm in}=0.1$, the distance is approximately $2.3\times 10^{-5}$
pc. To provide a sense of scale, we have indicated in Fig. 8 the distance to
the Sun of the outer planets as vertical dashed lines. From this, we observe
that to determine the orientation of the EPBH with respect to Earth, a close-
to-maximally rotating EPBH would need to be closer than Pluto’s aphelion. In
contrast, an EPBH with an initial $a_{\star}^{\rm in}=0.1$ would need to be
much closer, approximately as close as Jupiter is to the Sun, in order to
measure the neutrino-antineutrino emission asymmetry and enable the
determination of the EPBH hemisphere facing Earth.

_Future neutrino sensitivity. —_ Our previous estimates for determining the
EPBH’s initial characteristics using neutrino asymmetric emission were based
on IceCube’s current capabilities, particularly its effective area. However,
it is important to note that IceCube is set to undergo an upgrade, which is
anticipated to increase its effective area by approximately five times,
leading to a substantial enhancement in the detection of high-energy neutrinos
[120]. Furthermore, the development of new techniques, such as deep learning,
holds promise for refining the measurement of cascades. This advancement has
already demonstrated improved angular resolution for neutrino observations, as
evidenced by recent findings from the Galactic Center [121]. As a result, we
can reasonably expect that these innovations will significantly bolster the
statistics for nearby EPBH events, especially if they aid in measuring
additional neutrino flavors.

In addition to IceCube’s upgrade, a plethora of new experiments, including
KM3Net [122], P-ONE [123], Trident [124], and Baikal-GVD [125], are either
currently being built or in the planning stages. These experiments are also
expected to measure neutrinos in the TeV scale, and in the event of an EPBH
burst occurring near Earth, they will undoubtedly contribute to a substantial
increase in the data. Furthermore, the distribution of these new facilities
across different locations on Earth introduces the possibility of a
independent neutrino measurement of the Earth-EPBH distance if a nearby event
is observed. Although we do not quantitatively assess the extent to which all
these neutrino telescopes would improve the measurement of the initial EPBH
characteristics at the onset of its burst, we can anticipate a substantial
improvement in the measurement of the initial characteristics of an EPBH,
should one be observed in the vicinity of Earth.

##  V Conclusions

The detection of an evaporating black hole near Earth would represent an
extraordinary triumph for theoretical physics, validating our understanding of
quantum fields in curved spacetimes and providing insights into the existing
dofs in nature. Moreover, the detection of an evaporating black hole with an
initial non-zero spin would hint the existence of physics beyond the Standard
Model.

Neutrino emission from Kerr black holes exhibits unique behavior, primarily
due to the fact that only left-handed neutrinos and right-handed antineutrinos
interact weakly. Since particles with positive helicity are preferentially
emitted along the rotation axis of a Kerr black hole, while those with
negative helicity are predominantly emitted in the opposite direction,
neutrinos present an asymmetric emission: neutrinos are preferentially emitted
in the southern hemisphere, whereas antineutrinos are mostly emitted in the
northern hemisphere. In this work, we have proposed to exploit this neutrino
emission asymmetry as a powerful tool to determine the initial properties of
an EPBH. After considering the distribution of neutrinos and photons, we have
also derived, for the first time, the angular distribution of secondary
neutrinos and photons. We have analysed the time evolution of EPBHs and
calculated the full net neutrino flux integrated over the observed burst
duration. Our analysis assumed that the EPBH mass and spin follow the standard
time evolution, considering the dofs in the SM. An interesting aspect to
consider is that this time evolution may be influenced by additional BSM dofs,
particularly if they manifest as scalars. This could potentially result in a
distinct integrated neutrino-antineutrino emission asymmetry. Consequently, we
might expect to observe specific correlations between the detected neutrino
and antineutrino events, contingent on the specific BSM scenario at play. This
will be investigated in detail in future work.

We have computed the expected events in the IceCube observatory, finding that,
for an initially close-to-maximally rotating EPBH, the number of antineutrino
events would be larger than the neutrino events by a factor of $\sim 1.2$
($\sim 0.9$), if Earth is located at an angle of $\theta=15^{\circ}$
($120^{\circ}$) with respect to the rotation axis for a distance of $10^{-4}$
pc. Furthermore, through simultaneous measurements of photons and neutrinos
emitted from an EPBH, we have explored the possibility of not only determining
the initial black hole angular momentum but also the orientation of its axis
of rotation with respect to Earth. While gamma-ray events are expected to
significantly outnumber neutrino events, the symmetry in the angular
eigenfunctions of photon emission under the transformation
$\theta\to\pi-\theta$ makes it challenging to determine the black hole
hemisphere facing Earth. However, the net neutrino-antineutrino flux displays
a definite dependence on the polar angle, aiding in breaking the degeneracy
present in the photon measurement.

Taking optimistic assumptions regarding backgrounds and discrimination between
neutrinos and antineutrinos in IceCube, we found that, depending on the burst
duration, a neutrino measurement could lift the degeneracy in determining the
polar angle quadrant. To exclude one of the degenerate angle solutions in the
photon measurement at more than 95$\%$ CL, we deduced that the EPBH should be
within a distance of approximately $1.78\times 10^{-4}$ pc for initial spin
parameters of $a_{\star}^{\rm in}=0.5$ or $3.6\times 10^{-5}$ pc for
$a_{\star}^{\rm in}=0.1$. These distances are smaller than the Neptune-Sun or
Saturn-Sun distances, respectively. Consequently, perhaps an EPBH exists
within our solar system, approaching its final stages of life, presenting us
with a unique opportunity to directly observe Hawking radiation through a
multimessenger approach and investigate its properties before the final burst.

## Acknowledgments

The author is immensely thankful to Jessica Turner for her encouragement to
pursue this work and for engaging in several insightful discussions on various
topics related to this work. Moreover, the author would like to thank Jessica
Turner and Pedro Machado for their meticulous review of the manuscript and
their valuable comments, which significantly improved the quality of this
paper. In addition, special thanks are owed to Ivan Martínez-Soler for
generously providing assistance with queries related to IceCube, and to Sam
Dolan for valuable correspondence on the determination of the angular
eigenvalues for fermions. Lastly, the author is deeply grateful for the warm
hospitality extended by the Particle and Astroparticle Division of the Max-
Planck-Institute für Kernphysik, where a portion of this research was
finalized. This work 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.
This work has made use of the Hamilton HPC Service of Durham University.

##  Appendix A Fermions in a Kerr Spacetime

Let us briefly describe the separation of variables in the Dirac equation for
massive fermions on a Kerr background [92, 93, 94, 95]. The Dirac equation for
a fermion with mass $\mu$ is

$\displaystyle(i\gamma^{\alpha}\hat{D}_{\alpha}-\mu)\Psi=0,$ (28)

where $\hat{D}_{\alpha}=\partial_{\alpha}-\Gamma_{\alpha}$, $\Gamma_{\alpha}$
being the spin connection matrices and $\gamma^{\mu}$ are the Dirac matrices
in the curved spacetime satisfying the algebra

$\displaystyle\\{\gamma^{\mu},\gamma^{\nu}\\}=2g^{\mu\nu}I_{4\times 4},$ (29)

with $g^{\mu\nu}$ the inverse metric tensor associated to the Kerr metric. An
appropriate choice for these matrices for the Kerr spacetime in the Boyer-
Lindquist coordinates is

$\displaystyle\gamma^{t}$
$\displaystyle=\frac{r^{2}+a^{2}}{\rho\sqrt{\Delta}}\tilde{\gamma}^{0}+\frac{a\sin\theta}{\rho}\tilde{\gamma}^{2},$
$\displaystyle\gamma^{r}$
$\displaystyle=\frac{\sqrt{\Delta}}{\rho}\tilde{\gamma}^{3}$
$\displaystyle\gamma^{\phi}$
$\displaystyle=\frac{a}{\rho\sqrt{\Delta}}\tilde{\gamma}^{0}+\frac{1}{\rho\sin\theta}\tilde{\gamma}^{2},$
$\displaystyle\gamma^{\theta}$
$\displaystyle=\frac{1}{\sqrt{\rho}}\tilde{\gamma}^{1},$ (30)

where $\rho^{2}=r^{2}+a^{2}\cos\theta$, and $\tilde{\gamma}^{a}$ are the usual
Dirac matrices in a flat spacetime. We choose here the standard chiral
representation for $\tilde{\gamma}^{a}$,

$\displaystyle\tilde{\gamma}^{0}=\begin{pmatrix}0&I\\\
I&0\end{pmatrix},\quad\tilde{\gamma}^{i}=\begin{pmatrix}0&\sigma_{i}\\\
-\sigma_{i}&0\end{pmatrix}.$ (31)

Here $\sigma_{i}$ are the Pauli matrices, and $I$ is the $2\times 2$ identity
matrix. The chirality matrix $\gamma^{5}$ is defined in a standard manner,

$\displaystyle\gamma^{5}=i\tilde{\gamma}^{0}\tilde{\gamma}^{1}\tilde{\gamma}^{2}\tilde{\gamma}^{3}=\begin{pmatrix}-1&0\\\
0&1\end{pmatrix}.$ (32)

The spin-connection matrices are

$\displaystyle\Gamma_{\alpha}=\frac{1}{4}\omega_{\alpha
bc}\tilde{\gamma}^{b}\tilde{\gamma}^{c},$ (33)

being $\omega_{\alpha bc}$ the spin connection. Explicitly, the spin
connection matrices are

$\displaystyle\Gamma_{t}$
$\displaystyle=\frac{GM}{2}\begin{pmatrix}\varrho^{-2}\sigma_{3}&0\\\
0&-(\varrho^{*})^{-2}\sigma_{3}\end{pmatrix},$ (34a) $\displaystyle\Gamma_{r}$
$\displaystyle=-\frac{1}{2}\frac{a\sin\theta}{\sqrt{\Delta}}\begin{pmatrix}\varrho^{-1}\sigma_{2}&0\\\
0&-(\varrho^{*})^{-1}\sigma_{2}\end{pmatrix},$ (34b)
$\displaystyle\Gamma_{\theta}$
$\displaystyle=-\frac{1}{2}\sqrt{\Delta}\begin{pmatrix}(i\varrho)^{-1}\sigma_{2}&0\\\
0&-(i\varrho^{*})^{-1}\sigma_{2}\end{pmatrix},$ (34c)
$\displaystyle\Gamma_{\phi}$
$\displaystyle=\frac{1}{2}\sqrt{\Delta}\sin\theta\begin{pmatrix}(i\varrho)^{-1}\sigma_{1}&0\\\
0&-(i\varrho^{*})^{-1}\sigma_{1}\end{pmatrix}+\frac{1}{2}\vartheta\begin{pmatrix}\vartheta\sigma_{3}&0\\\
0&-\vartheta^{*}\sigma_{3}\end{pmatrix},$ (34d)

where $\varrho=r+ia\cos\theta$, and
$\vartheta=i\cos\theta-a\varrho^{-2}(\varrho+GM)\sin^{2}\theta$. The
contraction of the spin-connection with the Dirac matrices
$\gamma^{\alpha}\Gamma_{\alpha}$, corresponding to the term which appears in
the Dirac equation, is explicitly given by

$\displaystyle\gamma^{\alpha}\Gamma_{\alpha}=\frac{1}{2}\begin{pmatrix}0&s_{\theta}^{*}\sigma_{1}+s_{r}^{*}\sigma_{3}\\\
-s_{\theta}\sigma_{1}-s_{r}^{*}\sigma_{3}&0\end{pmatrix},$ (35)

where

$\displaystyle
s_{r}=\frac{\rho}{\sqrt{\Delta}}\frac{1}{\varrho\sqrt{\Delta}}\frac{\partial}{\partial
r}(\varrho\sqrt{\Delta}),\quad
s_{\theta}=\rho\frac{1}{\varrho\sin\theta}\frac{\partial}{\partial\theta}(\varrho\sin\theta).$

Following Ref. [95], we propose an ansatz to separate the Dirac equation,

$\displaystyle\Psi=\Delta^{-\frac{1}{4}}e^{-i\omega
t+im\phi}\begin{pmatrix}\varrho^{-1/2}\eta_{-}(r,\theta)\\\
(\varrho^{*})^{-1/2}\eta_{+}(r,\theta)\\\ \end{pmatrix},$ (36)

with the spinors $\eta_{\pm}$

$\displaystyle\eta_{-}(r,\theta)$
$\displaystyle=-\begin{pmatrix}R_{2}(r)\,\\!\,_{-\frac{1}{2}}S_{lm}(\theta)\\\
R_{1}(r)\,\\!\,_{+\frac{1}{2}}S_{lm}(\theta)\end{pmatrix},$ (37a)
$\displaystyle\eta_{+}(r,\theta)$
$\displaystyle=~{}~{}\begin{pmatrix}R_{1}(r)\,\\!\,_{-\frac{1}{2}}S_{lm}(\theta)\\\
R_{2}(r)\,\\!\,_{+\frac{1}{2}}S_{lm}(\theta)\end{pmatrix}.$ (37b)

Here, $R_{1,2}(r)$, $\\!\,{}_{\pm\frac{1}{2}}S_{lm}(\theta)$ are the radial
and angular functions associated to massive fermions. Substituting in the
Dirac equation, one finds the Eqs. (2) and (3).

##  Appendix B Chandrashekar – Detweiler Potentials

In a series of works, Chandrashekar and Detweiler proposed a method to
transform the general Teukolsky equations into Schrödinger-like equations
[106, 107, 108]. They found explicit forms for the potentials depending on the
spin of the particle. We have the potentials for scalars and vectors

$\displaystyle V_{0}(r)$
$\displaystyle=\frac{\Delta}{\rho^{4}}\left(\lambda_{0}+\frac{\Delta+2r(r-GM)}{\rho^{2}}-\frac{3r^{2}\Delta}{\rho^{4}}\right),$
(38a) $\displaystyle V_{1}(r)$
$\displaystyle=\frac{\Delta}{\rho^{4}}\left[\lambda_{1}+2-\alpha^{2}\frac{\Delta}{\rho^{4}}\pm
i\alpha\rho^{2}\frac{\differential}{\differential
r}\left(\frac{\Delta}{\rho^{4}}\right)\right],$ (38b)

where $\lambda_{s}$ are the angular eigenvalues, as defined in the text, and
$\alpha^{2}=a^{2}+am/\omega$. For vectors, there are two different values of
the potential, depending on the sign of the last term. We have chosen in our
numerical code the minus sign.

## References

  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ [2017] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Akiyama  _et al._ [2019] K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 875, L1 (2019), arXiv:1906.11238 [astro-ph.GA] .
  * Zel’dovich and Novikov [1967] Y. B. Zel’dovich and I. D. Novikov, Soviet Astron. AJ (Engl. Transl. ), 10, 602 (1967).
  * Hawking [1971] S. Hawking, Mon. Not. Roy. Astron. Soc. 152, 75 (1971).
  * Carr and Hawking [1974] B. J. Carr and S. W. Hawking, Mon. Not. Roy. Astron. Soc. 168, 399 (1974).
  * Hawking [1974] S. W. Hawking, Nature 248, 30 (1974).
  * Hawking [1975] S. W. Hawking,  _Euclidean quantum gravity_ , Commun. Math. Phys. 43, 199 (1975), [,167(1975)].
  * Carr  _et al._ [2010] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D 81, 104019 (2010), arXiv:0912.5297 [astro-ph.CO] .
  * Carr  _et al._ [2020] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, (2020), arXiv:2002.12778 [astro-ph.CO] .
  * Keith  _et al._ [2020] C. Keith, D. Hooper, N. Blinov, and S. D. McDermott, (2020), arXiv:2006.03608 [astro-ph.CO] .
  * He and Fang [2002] P. He and L.-Z. Fang, Astrophys. J. Lett. 568, L1 (2002), arXiv:astro-ph/0202218 .
  * Mack and Wesley [2008] K. J. Mack and D. H. Wesley, (2008), arXiv:0805.1531 [astro-ph] .
  * Cheek  _et al._ [2022a] A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 105, 015022 (2022a), arXiv:2107.00013 [hep-ph] .
  * Cheek  _et al._ [2022b] A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 105, 015023 (2022b), arXiv:2107.00016 [hep-ph] .
  * Hooper  _et al._ [2019] D. Hooper, G. Krnjaic, and S. D. McDermott, JHEP 08, 001 (2019), arXiv:1905.01301 [hep-ph] .
  * Masina [2020] I. Masina, Eur. Phys. J. Plus 135, 552 (2020), arXiv:2004.04740 [hep-ph] .
  * Morrison  _et al._ [2019] L. Morrison, S. Profumo, and Y. Yu, JCAP 1905, 005 (2019), arXiv:1812.10606 [astro-ph.CO] .
  * Auffinger  _et al._ [2021] J. Auffinger, I. Masina, and G. Orlando, Eur. Phys. J. Plus 136, 261 (2021), arXiv:2012.09867 [hep-ph] .
  * Khlopov  _et al._ [2006] M. Yu. Khlopov, A. Barrau, and J. Grain, Class. Quant. Grav. 23, 1875 (2006), arXiv:astro-ph/0406621 [astro-ph] .
  * Allahverdi  _et al._ [2018] R. Allahverdi, J. Dent, and J. Osinski, Phys. Rev. D97, 055013 (2018), arXiv:1711.10511 [astro-ph.CO] .
  * Lennon  _et al._ [2018] O. Lennon, J. March-Russell, R. Petrossian-Byrne, and H. Tillim, JCAP 1804, 009 (2018), arXiv:1712.07664 [hep-ph] .
  * Gondolo  _et al._ [2020] P. Gondolo, P. Sandick, and B. Shams Es Haghi, Phys. Rev. D 102, 095018 (2020), arXiv:2009.02424 [hep-ph] .
  * Baldes  _et al._ [2020] I. Baldes, Q. Decant, D. C. Hooper, and L. Lopez-Honorez, JCAP 08, 045 (2020), arXiv:2004.14773 [astro-ph.CO] .
  * Bernal and Zapata [2021a] N. Bernal and O. Zapata, JCAP 03, 015 (2021a), arXiv:2011.12306 [astro-ph.CO] .
  * Bernal and Zapata [2021b] N. Bernal and O. Zapata, Phys. Lett. B 815, 136129 (2021b), arXiv:2011.02510 [hep-ph] .
  * Masina [2021] I. Masina, Grav. Cosmol. 27, 315 (2021), arXiv:2103.13825 [gr-qc] .
  * Kitabayashi [2021] T. Kitabayashi, (2021), arXiv:2101.01921 [hep-ph] .
  * Bernal  _et al._ [2021] N. Bernal, Y. F. Perez-Gonzalez, Y. Xu, and O. Zapata, Phys. Rev. D 104, 123536 (2021), arXiv:2110.04312 [hep-ph] .
  * Bernal and Zapata [2021c] N. Bernal and O. Zapata, JCAP 03, 007 (2021c), arXiv:2010.09725 [hep-ph] .
  * Bernal  _et al._ [2022a] N. Bernal, Y. F. Perez-Gonzalez, and Y. Xu, (2022a), arXiv:2205.11522 [hep-ph] .
  * Cheek  _et al._ [2023] A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 108, 015005 (2023), arXiv:2212.03878 [hep-ph] .
  * Cheek  _et al._ [2022c] A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 106, 103012 (2022c), arXiv:2207.09462 [astro-ph.CO] .
  * Hooper  _et al._ [2020] D. Hooper, G. Krnjaic, J. March-Russell, S. D. McDermott, and R. Petrossian-Byrne, (2020), arXiv:2004.00618 [astro-ph.CO] .
  * Arbey  _et al._ [2021] A. Arbey, J. Auffinger, P. Sandick, B. Shams Es Haghi, and K. Sinha, (2021), arXiv:2104.04051 [astro-ph.CO] .
  * Perez-Gonzalez and Turner [2020] Y. F. Perez-Gonzalez and J. Turner, (2020), arXiv:2010.03565 [hep-ph] .
  * Jyoti Das  _et al._ [2021] S. Jyoti Das, D. Mahanta, and D. Borah, (2021), arXiv:2104.14496 [hep-ph] .
  * Datta  _et al._ [2020] S. Datta, A. Ghosal, and R. Samanta, (2020), arXiv:2012.14981 [hep-ph] .
  * Fujita  _et al._ [2014] T. Fujita, M. Kawasaki, K. Harigaya, and R. Matsuda, Phys. Rev. D89, 103501 (2014), arXiv:1401.1909 [astro-ph.CO] .
  * Granelli  _et al._ [2021] A. Granelli, K. Moffat, Y. F. Perez-Gonzalez, H. Schulz, and J. Turner, Comput. Phys. Commun. 262, 107813 (2021), arXiv:2007.09150 [hep-ph] .
  * Hook [2014] A. Hook, Phys. Rev. D 90, 083535 (2014), arXiv:1404.0113 [hep-ph] .
  * Hamada and Iso [2017] Y. Hamada and S. Iso, PTEP 2017, 033B02 (2017), arXiv:1610.02586 [hep-ph] .
  * Hooper and Krnjaic [2021] D. Hooper and G. Krnjaic, Phys. Rev. D 103, 043504 (2021), arXiv:2010.01134 [hep-ph] .
  * Chaudhuri and Dolgov [2020] A. Chaudhuri and A. Dolgov, (2020), arXiv:2001.11219 [astro-ph.CO] .
  * Bernal  _et al._ [2022b] N. Bernal, C. S. Fong, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 106, 035019 (2022b), arXiv:2203.08823 [hep-ph] .
  * Barman  _et al._ [2022] B. Barman, D. Borah, S. Jyoti Das, and R. Roshan, (2022), arXiv:2212.00052 [hep-ph] .
  * Gehrman  _et al._ [2022] T. C. Gehrman, B. Shams Es Haghi, K. Sinha, and T. Xu, (2022), arXiv:2211.08431 [hep-ph] .
  * Agashe  _et al._ [2023] K. Agashe, J. H. Chang, S. J. Clark, B. Dutta, Y. Tsai, and T. Xu, Phys. Rev. D 108, 023014 (2023), arXiv:2212.11980 [hep-ph] .
  * Gehrman  _et al._ [2023] T. C. Gehrman, B. Shams Es Haghi, K. Sinha, and T. Xu, (2023), arXiv:2304.09194 [hep-ph] .
  * Friedlander  _et al._ [2023] A. Friedlander, N. Song, and A. C. Vincent, (2023), arXiv:2306.01520 [hep-ph] .
  * Carr  _et al._ [2016] B. Carr, F. Kuhnel, and M. Sandstad, Phys. Rev. D 94, 083504 (2016), arXiv:1607.06077 [astro-ph.CO] .
  * Green and Kavanagh [2021] A. M. Green and B. J. Kavanagh, J. Phys. G 48, 4 (2021), arXiv:2007.10722 [astro-ph.CO] .
  * Villanueva-Domingo  _et al._ [2021] P. Villanueva-Domingo, O. Mena, and S. Palomares-Ruiz, Front. Astron. Space Sci. 8, 87 (2021), arXiv:2103.12087 [astro-ph.CO] .
  * Friedlander  _et al._ [2022] A. Friedlander, K. J. Mack, S. Schon, N. Song, and A. C. Vincent, Phys. Rev. D 105, 103508 (2022), arXiv:2201.11761 [hep-ph] .
  * Auffinger [2022] J. Auffinger, (2022), arXiv:2206.02672 [astro-ph.CO] .
  * Glicenstein  _et al._ [2013] J.-F. Glicenstein, A. Barnacka, M. Vivier, and T. Herr (H.E.S.S.), in  _33rd International Cosmic Ray Conference_ (2013) p. 0930, arXiv:1307.4898 [astro-ph.HE] .
  * Tavernier  _et al._ [2020] T. Tavernier, J.-F. Glicenstein, and F. Brun, PoS ICRC2019, 804 (2020), arXiv:1909.01620 [astro-ph.HE] .
  * Abdo  _et al._ [2015] A. A. Abdo _et al._ , Astropart. Phys. 64, 4 (2015), arXiv:1407.1686 [astro-ph.HE] .
  * Archambault [2018] S. Archambault (VERITAS), PoS ICRC2017, 691 (2018), arXiv:1709.00307 [astro-ph.HE] .
  * Abeysekara  _et al._ [2013] A. U. Abeysekara _et al._ (HAWC), (2013), arXiv:1310.0073 [astro-ph.HE] .
  * Albert  _et al._ [2020] A. Albert _et al._ (HAWC), JCAP 04, 026 (2020), arXiv:1911.04356 [astro-ph.HE] .
  * Ackermann  _et al._ [2018] M. Ackermann _et al._ (Fermi-LAT), Astrophys. J. 857, 49 (2018), arXiv:1802.00100 [astro-ph.HE] .
  * Halzen  _et al._ [1995] F. Halzen, B. Keszthelyi, and E. Zas, Phys. Rev. D 52, 3239 (1995), arXiv:hep-ph/9502268 .
  * Dave and Taboada [2021] P. Dave and I. Taboada (IceCube), PoS ICRC2019, 863 (2021), arXiv:1908.05403 [astro-ph.HE] .
  * Capanema  _et al._ [2021] A. Capanema, A. Esmaeili, and A. Esmaili, JCAP 12, 051 (2021), arXiv:2110.05637 [hep-ph] .
  * Ukwatta  _et al._ [2016] T. N. Ukwatta, D. R. Stump, J. T. Linnemann, J. H. MacGibbon, S. S. Marinelli, T. Yapici, and K. Tollefson, Astropart. Phys. 80, 90 (2016), arXiv:1510.04372 [astro-ph.HE] .
  * Baker and Thamm [2022] M. J. Baker and A. Thamm, SciPost Phys. 12, 150 (2022), arXiv:2105.10506 [hep-ph] .
  * Baker and Thamm [2023] M. J. Baker and A. Thamm, JHEP 01, 063 (2023), arXiv:2210.02805 [hep-ph] .
  * Boluna  _et al._ [2023] X. Boluna, S. Profumo, J. Blé, and D. Hennings, (2023), arXiv:2307.06467 [astro-ph.HE] .
  * Page [1976a] D. N. Page, Phys. Rev. D 13, 198 (1976a).
  * Page [1976b] D. N. Page, Phys. Rev. D 14, 3260 (1976b).
  * Arvanitaki  _et al._ [2010] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, Phys. Rev. D 81, 123530 (2010), arXiv:0905.4720 [hep-th] .
  * Chambers  _et al._ [1997] C. M. Chambers, W. A. Hiscock, and B. Taylor, Phys. Rev. Lett. 78, 3249 (1997), arXiv:gr-qc/9703018 .
  * Taylor  _et al._ [1998] B. E. Taylor, C. M. Chambers, and W. A. Hiscock, Phys. Rev. D 58, 044012 (1998), arXiv:gr-qc/9801044 .
  * Calzà  _et al._ [2021] M. Calzà, J. March-Russell, and J. a. G. Rosa, (2021), arXiv:2110.13602 [astro-ph.CO] .
  * Calzà and Rosa [2022] M. Calzà and J. a. G. Rosa, JHEP 12, 090 (2022), arXiv:2210.06500 [gr-qc] .
  * NASA [1968] NASA, “Earthrise,”  (1968).
  * Arbey and Auffinger [2019] A. Arbey and J. Auffinger, Eur. Phys. J. C 79, 693 (2019), arXiv:1905.04268 [gr-qc] .
  * Arbey and Auffinger [2021] A. Arbey and J. Auffinger, Eur. Phys. J. C 81, 910 (2021), arXiv:2108.02737 [gr-qc] .
  * Unruh [1973] W. Unruh, Phys. Rev. Lett. 31, 1265 (1973).
  * Vilenkin [1978] A. Vilenkin, Phys. Rev. Lett. 41, 1575 (1978).
  * Leahy and Unruh [1979] D. A. Leahy and W. G. Unruh, Phys. Rev. D 19, 3509 (1979).
  * Vilenkin [1979] A. Vilenkin, Phys. Rev. D 20, 1807 (1979).
  * Vilenkin [1980] A. Vilenkin, Phys. Rev. D 22, 3067 (1980).
  * Aartsen  _et al._ [2019] M. G. Aartsen _et al._ (IceCube), Phys. Rev. D 99, 032004 (2019), arXiv:1808.07629 [hep-ex] .
  * Aartsen  _et al._ [2020] M. G. Aartsen _et al._ (IceCube), Phys. Rev. Lett. 124, 051103 (2020), arXiv:1910.08488 [astro-ph.HE] .
  * Boyer and Lindquist [1967] R. H. Boyer and R. W. Lindquist, J. Math. Phys. 8, 265 (1967).
  * Teukolsky [1973] S. A. Teukolsky, Astrophys. J. 185, 635 (1973).
  * Press and Teukolsky [1973] W. H. Press and S. A. Teukolsky, Astrophys. J. 185, 649 (1973).
  * Teukolsky and Press [1974] S. A. Teukolsky and W. H. Press, Astrophys. J. 193, 443 (1974).
  * Lunardini and Perez-Gonzalez [2020] C. Lunardini and Y. F. Perez-Gonzalez, JCAP 08, 014 (2020), arXiv:1910.07864 [hep-ph] .
  * Chandrasekhar [1976] S. Chandrasekhar, Proc. Roy. Soc. Lond. A 349, 571 (1976).
  * Page [1976c] D. N. Page, Phys. Rev. D 14, 1509 (1976c).
  * Dolan and Gair [2009] S. Dolan and J. Gair, Class. Quant. Grav. 26, 175020 (2009), arXiv:0905.2974 [gr-qc] .
  * Dolan and Dempsey [2015] S. R. Dolan and D. Dempsey, Class. Quant. Grav. 32, 184001 (2015), arXiv:1504.03190 [gr-qc] .
  * Page [1977] D. N. Page, Phys. Rev. D 16, 2402 (1977).
  * Batic  _et al._ [2005] D. Batic, H. Schmid, and M. Winklmeier, J. Math. Phys. 46, 012504 (2005), arXiv:math-ph/0402047 .
  * Batic and Schmid [2008] D. Batic and H. Schmid, Rev. Col. Mat. 42, 183 (2008), arXiv:gr-qc/0512112 .
  * Neznamov and Safronov [2016] V. P. Neznamov and I. I. Safronov, Int. J. Mod. Phys. D 25, 1650091 (2016), arXiv:1605.07450 [gr-qc] .
  * Casals  _et al._ [2007] M. Casals, S. R. Dolan, P. Kanti, and E. Winstanley, JHEP 03, 019 (2007), arXiv:hep-th/0608193 .
  * Casals  _et al._ [2013] M. Casals, S. R. Dolan, B. C. Nolan, A. C. Ottewill, and E. Winstanley, Phys. Rev. D 87, 064027 (2013), arXiv:1207.7089 [gr-qc] .
  * Winstanley [2016] E. Winstanley, Springer Proc. Phys. 170, 291 (2016), arXiv:1310.5064 [gr-qc] .
  * Unruh [1976] W. G. Unruh, Phys. Rev. D 14, 870 (1976).
  * Kayser [1982] B. Kayser, Phys. Rev. D 26, 1662 (1982).
  * Kayser and Shrock [1982] B. Kayser and R. E. Shrock, Phys. Lett. B 112, 137 (1982).
  * Chandrasekhar and Detweiler [1975] S. Chandrasekhar and S. L. Detweiler, Proc. Roy. Soc. Lond. A 345, 145 (1975).
  * Chandrasekhar and Detweiler [1976] S. Chandrasekhar and S. L. Detweiler, Proc. Roy. Soc. Lond. A 350, 165 (1976).
  * Chandrasekhar and Detweiler [1977] S. Chandrasekhar and S. L. Detweiler, Proc. Roy. Soc. Lond. A 352, 325 (1977).
  * Leaver [1985] E. W. Leaver, Proc. Roy. Soc. Lond. A 402, 285 (1985).
  * Berti  _et al._ [2006] E. Berti, V. Cardoso, and M. Casals, Phys. Rev. D 73, 024013 (2006), [Erratum: Phys.Rev.D 73, 109902 (2006)], arXiv:gr-qc/0511111 .
  * Duffy  _et al._ [2005] G. Duffy, C. Harris, P. Kanti, and E. Winstanley, JHEP 09, 049 (2005), arXiv:hep-th/0507274 .
  * Casals  _et al._ [2006] M. Casals, P. Kanti, and E. Winstanley, JHEP 02, 051 (2006), arXiv:hep-th/0511163 .
  * Eingorn  _et al._ [2018] M. Eingorn, L. Fernando, B. Vlahovic, C. Ilie, B. Wojtsekhowski, G. M. Urciuoli, F. D. Persio, F. Meddi, and V. Nelyubin, J. Astron. Telesc. Instrum. Syst. 4, 011006 (2018), arXiv:1501.05592 [astro-ph.HE] .
  * Depaola  _et al._ [1999] G. O. Depaola, C. N. Kozameh, and M. H. Tiglio, Astroparticle Physics 10, 175 (1999).
  * Morselli  _et al._ [2013] A. Morselli _et al._ , Nucl. Phys. B Proc. Suppl. 239-240, 193 (2013), arXiv:1406.1071 [astro-ph.IM] .
  * Bernard [2013] D. Bernard (HARPO), Nucl. Instrum. Meth. A 729, 765 (2013), arXiv:1307.3892 [astro-ph.IM] .
  * Gros and Bernard [2017] P. Gros and D. Bernard, Astroparticle Physics 88, 30 (2017).
  * Bierlich  _et al._ [2022] C. Bierlich _et al._ ,  (2022), 10.21468/SciPostPhysCodeb.8, arXiv:2203.11601 [hep-ph] .
  * Carr and Kuhnel [2020] B. Carr and F. Kuhnel, Ann. Rev. Nucl. Part. Sci. 70, 355 (2020), arXiv:2006.02838 [astro-ph.CO] .
  * Clark [2021] B. Clark (IceCube-Gen2), JINST 16, C10007 (2021), arXiv:2108.05292 [astro-ph.HE] .
  * Abbasi  _et al._ [2023] R. Abbasi _et al._ (IceCube), Science 380, 6652 (2023), arXiv:2307.04427 [astro-ph.HE] .
  * Adrian-Martinez  _et al._ [2016] S. Adrian-Martinez _et al._ (KM3Net), J. Phys. G 43, 084001 (2016), arXiv:1601.07459 [astro-ph.IM] .
  * Agostini  _et al._ [2020] M. Agostini _et al._ (P-ONE), Nature Astron. 4, 913 (2020), arXiv:2005.09493 [astro-ph.HE] .
  * Ye  _et al._ [2022] Z. P. Ye _et al._ , (2022), arXiv:2207.04519 [astro-ph.HE] .
  * Allakhverdyan  _et al._ [2023] V. A. Allakhverdyan _et al._ (Baikal-GVD), Phys. Rev. D 107, 042005 (2023), arXiv:2211.09447 [astro-ph.HE] .