# The Amaterasu particle: constraining the superheavy dark matter origin of
UHECRs

Prantik Sarmah<EMAIL_ADDRESS>Indian Institute of Technology Guwahati,
Guwahati, Assam-781039, India  Nayan Das<EMAIL_ADDRESS>Indian
Institute of Technology Guwahati, Guwahati, Assam-781039, India  Debasish
Borah<EMAIL_ADDRESS>Indian Institute of Technology Guwahati, Guwahati,
Assam-781039, India  Sovan Chakraborty<EMAIL_ADDRESS>Indian Institute of
Technology Guwahati, Guwahati, Assam-781039, India  Poonam Mehta<EMAIL_ADDRESS>School of Physical Sciences, Jawaharlal Nehru University, New Delhi-110067,
India

###### Abstract

Amaterasu, the second most energetic ($244$ EeV) cosmic ray particle has been
recently detected by the Telescope Array (TA) surface detector. The origin of
the TA Amaterasu event is puzzling, as its arrival direction points back to a
void in the local Universe, lacking conventional astrophysical ultra-high-
energy (UHE) cosmic ray (CR) sources. Hence, we explore the possibility if
this TA Amaterasu event could have originated from the decay of superheavy
dark matter (SHDM) in the Milky Way. Such an origin also opens up multi-
messenger detection channels in both UHE gamma-rays and UHE neutrinos. In this
present work, using the TA Amaterasu event and the multi-messenger
limits/sensitivities from various UHE telescopes, we place stringent
constraints on the lifetime and mass of the SHDM. We find that the non-
detection of the corresponding gamma-rays at the Pierre Auger Observatory
(PAO) and the TA is in severe tension with the SHDM parameter space required
to explain the TA Amaterasu event. Additionally, we extend the multi-messenger
analysis to the future UHE gamma-ray and UHE neutrino telescopes such as PAO
upgrade, GRAND 200k and IceCube-Gen2. We find that the bounds from the future
neutrino telescopes will be able to compete with the present UHECR bounds.
However, compared to the existing UHE gamma-ray bounds, the future PAO upgrade
and the GRAND 200k gamma-ray detectors will improve the bounds on SHDM
lifetime by at least one order of magnitude.

††preprint: APS/123-QED

##  I Introduction

The Telescope Array (TA) experiment has recently detected an ultra-high-energy
cosmic ray (UHECR) event (named Amaterasu) of energy $2.44\times 10^{11}$ GeV
which is the second highest energy UHECR ever detected [1]. The arrival
direction (R.A.: $255.9\pm 0.6^{\circ}$, Dec.: $16.1\pm 0.5^{\circ}$) of this
UHECR points back to a void in the local Universe rendering more questions
than answers about its origin. As UHECRs above energy $5\times 10^{10}$ GeV
propagating more than $\sim 10$ Mpc (mean free path) undergo energy losses due
to the Greisen–Zatsepin–Kuzmin (GZK) effect [2, 3, 4, 5, 6, 7], the likely
source of this event is expected to be located within a distance of $\sim 30$
Mpc [1]. The TA collaboration [1] analysed a catalogue of gamma-ray sources to
identify the probable source and found that the known sources were [8, 9]
located at distances beyond $30$ Mpc. One possibility could be that the
primary UHECR is an iron nucleus at source, which might have undergone
significant deflection in the extragalactic magnetic field (EGMF) during
propagation to Earth [10, 11, 12]. Due to this, the actual source direction
might be different than the arrival direction. Magnetic monopole has also been
explored in the literature as a possible origin of the TA Amaterasu event [13,
14, 15, 16, 17, 18]. Another possible explanation of the origin could be
attributed to decaying superheavy dark matter (SHDM) in the Milky Way or any
local galaxy along the arrival direction. In this work, we examine if such a
possibility of decaying SHDM in the Milky Way could give rise to the TA
Amaterasu event.

Contrary to the conventional production of ultra-high-energy cosmic rays
through the diffusive shock acceleration (bottom-up), dark matter (DM) decay
into standard model (SM) particles has also been hypothesised as a top-down
UHECR production mechanism (See review[7]). In order to produce the observed
energy of the TA Amaterasu event via DM decay, energy-momentum conservation
states that the DM particles should have mass larger than $\sim 5\times
10^{11}$ GeV and can extend up to the energy scale ($10^{15}-10^{16}$ GeV),
close to the typical grand unified theory (GUT) scale [19, 20, 21, 22, 23,
24]. Production of such SHDM may be possible via mechanisms such as vacuum or
gravitational fluctuations in the early Universe during or after inflation
[25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 22, 36, 37, 38, 39]. Note that
UHECRs can also be produced in principle from self-annihilation of SHDM
particles. However, the resultant flux is density squared suppressed [40]. Our
analysis is based on the Galactic SHDM decay possibility, the extra-galactic
component of the flux is found to be negligible in comparison to the the
dominant contribution from Galactic DM halo [41].

SHDM particle ($\chi$) in the Milky Way halo can produce UHECR protons ($p$),
ultra-high-energy (UHE) gamma-rays ($\gamma$) and UHE neutrinos ($\nu$) via
the decay mode $\chi\to X\bar{X}$. Here, $X$ and $\bar{X}$ are the SM
particle-antiparticle pairs [42]. For example, SHDM decaying via $b$ quark
channel can produce these UHE secondaries: $\chi\to b\bar{b}\to
p\bar{p}~{}\text{or}~{}\gamma\gamma~{}\text{or}~{}\nu\bar{\nu}$. The resulting
abundances of gamma-rays and neutrinos generally dominate over that of protons
for specific choice of the decay channel. This is in disagreement with the
UHECR composition measurements [43, 44, 45], that hint towards dominance of
nucleon over gamma-rays and neutrinos. Note that abundances of different
nuclei from protons to iron in the UHECR spectrum is still an unresolved issue
[43, 44, 45]. Because of the disagreement with observations, SHDM scenarios
for UHECRs are tightly constrained up to the maximum observed CR energy, $\sim
10^{11}$ GeV. The conventional UHECR spectrum is heavily suppressed around
$10^{11}$ GeV due to the GZK effect and falls with energy as $\sim E^{-5}$
[46]. Thus, the region above this energy becomes invaluable for probing
scenarios of SHDM.

The UHE gamma-ray and UHE neutrino fluxes associated with the conventional
astrophysical UHECR flux are generally much smaller than that of the UHECR
flux [4]. Whereas in SHDM scenarios, we expect them to be larger as mentioned
above. This distinct multi-messenger characteristic can be a useful probe of
UHECRs from SHDM and therefore, detection of UHE gamma-rays and UHE neutrinos
becomes crucial. There are mainly two different telescopes, Pierre Auger
Observatory (PAO) [47, 48, 49, 43, 46, 50, 51] and TA [52, 44, 53, 45, 54, 1]
currently operating and capable of detecting all the three components, UHECRs,
UHE gamma-rays, and UHE neutrinos. The Yakutsk Extensive Air Shower Array can
detect UHECRs and gamma-rays [55, 56], but its sensitivity is lower compared
to TA and PAO. On the other hand, the IceCube experiment is also capable of
detecting UHE neutrinos and holds sensitivity similar to that of PAO. The
future telescope proposals such as Giant Radio Array for Neutrino Detection
200k (GRAND 200k) [57] and IceCube-Gen2 [58] and PAO upgrade [48] having much
larger sensitivities compared to the current ones will provide us with
unprecedented opportunities to conduct multi-messenger analysis at these
extreme energies.

In this work, we for the first time consider the origin of the recent TA
Amaterasu event in light of SHDM decay. We put multi-messenger constraints on
the lifetime and mass of the SHDM using UHECR, gamma-ray and neutrino
telescopes. We also account for the low energy UHECR spectrum detected by the
TA. We consider different decay modes for SHDM as well as a wide range
($10^{11}-10^{15}$ GeV) of SHDM mass.  The constraints on the lifetime from
gamma-ray non-detection limit of TA and PAO are found to be stronger than that
of the UHECRs, discrediting the possibility of the TA Amaterasu event
originating from SHDM decay. Furthermore, we also analyse the potential of the
future telescopes mentioned above for probing SHDM lifetime through UHE gamma-
ray and UHE neutrino detection. While future UHE neutrino detection are found
to yield constraints comparable to that from UHECRs, UHE gamma-ray detection
will lead to a few order of magnitudes improved constraints than that from
UHECRs. Our constraints based on the TA Amaterasu event are consistent with
those existing in the literature [41, 40, 59, 60, 61, 62, 63, 64, 65, 66, 50,
67, 68, 69, 70, 71, 72, 73, 74, 31, 73, 75, 76]. These works are based on
various observations of UHECRs, UHE gamma-rays, and UHE neutrinos, but not the
recent TA Amaterasu event.

The paper is organised as follows. In Sec. II, we describe the model adopted
for flux calculation . In Sec. III, we present the results of our multi-
messenger analysis for both present and future telescopes. Finally we
summarize and discuss our main results in Sec. IV. The dependence of the
results on various DM density profiles is discussed in Appendix A.

##  II UHE flux from SHDM decay

In this section, we discuss the production of UHECRs, UHE gamma-rays and UHE
neutrinos from the decay of SHDM and estimate the flux of these daughter
particles. The SHDM particle, $\chi$ is considered to decay via the
$1\rightarrow 2$ process ($\chi\rightarrow X\bar{X}$) into various SM
particles-antiparticle pairs. These SM particles and antiparticles can further
produce other lighter particle-antiparticle pairs. In particular, we focus on
the secondaries convenient for detection and further phenomenology, i.e., the
UHECR protons (p), the UHE gamma-rays ($\gamma$) and the UHE neutrinos
($\nu$).

The differential energy spectra, $dN_{i}/dE$ (where, $i\equiv\text{p,
$\gamma$, and $\nu$}$) of these UHE final state particles depend on the SHDM
mass ($m_{\chi}$), lifetime ($\tau_{\chi}$) and the inclusive decay rate
($\Gamma$) of $\chi\to X\bar{X}\to i+...$ [42]. The flux of an observable
particle $i$ produced from DM decay depends on the differential energy spectra
given by,

$\displaystyle\dfrac{dN_{i}}{dx}=\dfrac{1}{\Gamma_{0}}\dfrac{d\Gamma}{dx}(\chi\to
i+\ldots)\ .$ (1)

Here, $x=2E/m_{\chi}$ and $\Gamma_{0}=1/\tau_{\chi}$ is the inverse lifetime
in the rest frame of $\chi$ [42]. The lifetime ($\tau_{\chi}$) and the mass
($m_{\chi}$) of SHDM are considered to be the free parameters of the model and
are constrained from the various observations of the UHE secondaries (p,
$\gamma$, and $\nu$).

In particular, the differential energy spectra ($dN_{i}/dE$) for different
decay modes ($\chi\rightarrow X\bar{X}\rightarrow i+...$) are computed with
the publicly available numerical code HDMSpectra [42]. HDMSpectra allows to
compute energy spectra for a wide $m_{\chi}$ range, from the electroweak scale
up to the GUT scale. The energy spectra are estimated by the fragmentation
function, $D_{a}^{b}(x;\mu_{Q},\mu_{0})$ that gives the probability of an
initial state $a$ (e.g., $X\bar{X}$) at energy scale $\mu_{Q}$ evolving to a
final state $b$ (e.g., $\rm p+\bar{p}$) at energy scale $\mu_{0}$ with
momentum fraction $x=2E/m_{\chi}$ (see [42] for details).

Figure 1: UHE flux fitted to different observations of TA for
$m_{\chi}=10^{12}$ GeV (solid) and $m_{\chi}=10^{14}$ GeV (dashed). The decay
channel is taken to be $\chi\to b\bar{b}$. Left plot: UHECR proton flux (blue
curves) from SHDM decay fitted to the TA Amaterasu event (red data point) and
the UHECR spectrum of TA (black data point). The UHECR flux above the energy
of the TA Amaterasu event is taken to be constant (purple dotted). The flux
fitting resulted in $\tau_{\chi}=5.5\times 10^{28}$ s and $7.5\times 10^{29}$
s for $m_{\chi}=10^{12}$ GeV and $m_{\chi}=10^{14}$ GeV, respectively. Middle
plot: Integral UHE gamma-ray flux (magenta curves) from SHDM decay for the two
$m_{\chi}$ values fitted to the upper limit on UHE gamma-ray flux from TA
(black downward arrow). The resultant values of $\tau_{\chi}$ are $3\times
10^{30}$ s and $8\times 10^{29}$ s, respectively. Right plot: UHE neutrino
flux (green curves) from SHDM decay for the two $m_{\chi}$ values fitted to
the upper limit on UHE neutrino flux from TA (black line). The corresponding
values of $\tau_{\chi}$ are $7\times 10^{25}$ s and $4.8\times 10^{25}$ s,
respectively. The UHE neutrino flux constraint infers significantly small
$\tau_{\chi}$ values compared to that of UHE protons and UHE gamma-rays.

Now, to estimate the total diffuse secondary particle flux from these SHDM
decays one needs to include the contribution from all the SHDM distributed all
over the Universe. In principle, the SHDM abundance over the Universe could be
crucial. However, it turns out that the extra-galactic contribution is
negligible in comparison to the secondary flux from SHDM decay in the Milky
Way [41]. Thus, we focus on the diffuse flux of UHE protons, UHE gamma-rays
and UHE neutrinos produced from the Galactic abundance of SHDM. The DM density
in the Milky Way is expected to be the largest at the Galactic center and
decreases as a function of the Galactocentric distance ($R_{\rm GC}$). There
are various models and simulations for the density profile of DM in the Milky
Way [77, 78, 79, 80, 81]. In this work, we consider the Navarro-Frenk-White
(NFW) profile [77] given by,

$\rho_{\rm NFW}(R_{\rm GC})=\frac{\rho_{\rm c}}{(R_{\rm GC}/R_{\rm
c})(1+R_{\rm GC}/R_{\rm c})^{2}}\ ,$ (2)

where, $\rho_{\rm c}$ is the characteristic density at the characteristic
scale $R_{\rm c}=11$ kpc and its value is determined by normalizing $\rho_{\rm
NFW}$ with the DM density in the solar neighbourhood, $\rho_{\odot}=0.43~{}\rm
GeV~{}cm^{-3}$.

Thus, given mass ($m_{\chi}$) and lifetime ($\tau_{\chi}$) of the SHDM
particle, the total flux, $\Phi_{i}(E,\theta)$ of the $i^{th}$ secondary UHE
particle ($i=\text{p, $\gamma$, and $\nu$}$) along the line of sight can be
expressed as,

$\begin{split}\Phi_{i}(E,\theta)&=\frac{1}{4\pi
m_{\chi}\tau_{\chi}}\frac{dN_{i}(E)}{dE}\int_{0}^{l_{\rm
max}(\theta)}\rho_{\chi}(l)dl\\\ &=\frac{\rho_{\odot}R_{\rm SGC}}{4\pi
m_{\chi}\tau_{\chi}}\frac{dN_{i}(E)}{dE}\mathcal{J}(\theta)\ ,\end{split}$ (3)

where,

$\mathcal{J}(\theta)=\frac{1}{\rho_{\odot}R_{\rm SGC}}\int_{0}^{l_{\rm
max}(\theta)}\rho_{\chi}(l)dl\ .$ (4)

Here, $R_{\rm SGC}$ is the distance between the Sun and the Galactic center
and is taken to be $8.34$ kpc. This allows changing the line of sight
coordinate ($l$) to the convenient Galactocentric coordinate ($R_{\rm GC}$).
The coordinate transformation from $l$ to $R_{\rm GC}$ is given by $R_{\rm
GC}=\sqrt{R_{\rm SGC}^{2}-2lR_{\rm SGC}\cos{\theta}+l^{2}}$.

The flux over a solid angle, $\Omega=2\pi(1-\cos{\theta})$ can be obtained by
integrating out Eq. (4) from $\theta=0$ to $\theta=\theta_{\rm max}$. Thus the
angle integrated total flux of the $i^{th}$ secondary is,

$\Phi_{i}(E)=\frac{\rho_{\odot}R_{\rm SGC}}{4\pi
m_{\chi}\tau_{\chi}}\frac{dN_{i}(E)}{dE}\mathcal{J}^{\theta_{\rm max}}\ ,$ (5)

where,

$\mathcal{J}^{\theta_{\rm max}}=\frac{2\pi}{\Omega}\int_{0}^{\theta_{\rm
max}}\sin{\theta}d\theta\mathcal{J}(\theta)\ .$ (6)

As we are interested in the diffuse flux of UHECRs, UHE neutrinos and UHE
gamma-rays, we integrate Eq. (6) over the entire solid angle, i.e.,
$\theta_{\rm max}=\pi$ and obtain $\mathcal{J}^{\theta_{\rm max}}\approx 2$.
The dependence of the fluxes of UHE secondaries on a specific DM density
profile can be estimated by calculating the $\mathcal{J}^{\theta_{\rm max}}$
factor. Note that, a different choice of DM density profile other than the NFW
profile leads to a negligible change in the UHE secondary fluxes, see Appendix
A.

Now, this diffuse flux of the secondary particles can be obtained for
different decay modes of $\chi$, i.e., for different quarks, leptons and
bosons for the $X$ and $\bar{X}$ combination (eg., $b\bar{b}$, $W^{+}W^{-}$,
$\mu^{+}\mu^{-}$, $\nu\bar{\nu}$) in the decay process $\chi\to X\bar{X}$. In
particular, the differential energy spectra ${dN_{i}(E)}/{dE}$ is estimated
for these different decay modes.

To demonstrate the flux estimation for the UHECRs, UHE gamma-rays and UHE
neutrinos, we consider the example of $\chi\to b\bar{b}$ mode using Eq. (5)
with $100\%$ branching fraction. The results are depicted in Fig. 1 together
with the observational data and constraints of these UHE fluxes from the TA
experiment. We consider two specific choices of $m_{\chi}$: $10^{12}$ GeV and
$10^{14}$ GeV. The fluxes of all these UHE particles are inversely
proportional (see Eq. (5) ) to $\tau_{\chi}$. Rather than fixing the parameter
$\tau_{\chi}$ arbitrarily, we constrain the value of $\tau_{\chi}$ in a
phenomenological manner. $\tau_{\chi}$ is estimated by normalising the SHDM
derived flux of UHECRs, UHE gamma-rays and UHE neutrinos with the observed
data or experimental upper limits on these fluxes. Since these fluxes as well
as the observed flux/limits are all different from each other, they give rise
to different values of $\tau_{\chi}$.

The left plot of Fig. 1 shows the observed UHECR spectrum data by TA in black
dots. The red data point represents the flux of the recently observed TA
Amaterasu event of energy $2.44\times 10^{11}$ GeV [1]. The flux for this TA
Amaterasu event is obtained by dividing the total number of events by the
total exposure of the telescope [1]. The uncertainties in the flux correspond
to the propagated uncertainties. Importantly, for the energies higher than
$244$ EeV, we consider a conservative estimate that the maximum possible CR
flux is not larger than the observed TA flux at the $244$ EeV event, shown by
the purple dashed line. This is contrary to the harder CR spectra suggestions
at the highest energies [43, 82]. For these maximum limits, we show the upper
limit of $p+\bar{p}$ flux for two different DM masses ($m_{\chi}=10^{12}$ GeV
and $m_{\chi}=10^{14}$ GeV) considering the decay of DM into $b\bar{b}$
channel. This establishes the lower bound on DM lifetime, as indicated by Eq.
(5). For the DM masses of $10^{12}$ GeV and $10^{14}$ GeV, the observational
constraints fix the $\tau_{\chi}$ at $5.5\times 10^{28}$ s and $7.5\times
10^{29}$ s, respectively. Note that the required energy for the TA Amaterasu
event can not be created from SHDM with $m_{\chi}\lesssim 10^{12}$ GeV. This
can be understood from the flux corresponding to $m_{\chi}=10^{12}$ GeV in the
left plot of Fig. 1 as the TA Amaterasu event lies on the high energy tail of
this flux. For $m_{\chi}\lesssim 10^{12}$ GeV, the flux will fall off at
energies smaller than the energy of the TA Amaterasu event.  Similarly, there
exists an upper limit for the SHDM mass. For example, the flux corresponding
to $m_{\chi}=10^{14}$ GeV is falling below the TA flux for our assumption of
constant UHECR limit beyond the energy of the TA Amaterasu event. In
particular, we find an upper limit on $m_{\chi}$ for this specific channel to
be $m_{\chi}\sim 4\times 10^{13}$ GeV.

In the middle plot of Fig. 1, we show the upper limit on the integral UHE
gamma-ray flux from the TA experiment [53, 54]. This upper limit comes from
the non-detection of UHE gamma-rays. The solid and dashed magenta lines
represent the gamma-ray flux for DM masses of $10^{12}$ GeV and $10^{14}$ GeV,
the resultant $\tau_{\chi}$ to maintain the TA gamma-ray upper limits are
$3\times 10^{30}$ s and $8\times 10^{29}$ s, respectively. Again, the decay of
DM is considered to be into $b\bar{b}$ final states. These values of
$\tau_{\chi}$ are significantly larger than those obtained from UHECR
constraints. This is due to the fact that production of gamma-rays from
$b\bar{b}$ annihilation is significantly enhanced compared to that of protons
[42]. Similarly, in the right panel plot of Fig. 1, we show the TA’s upper
limit (black continuous line) on the UHE neutrino flux from the non-detection
of UHE neutrinos in TA [83]. The solid and dashed green lines correspond to
the flux of UHE neutrinos for the DM masses of $10^{12}$ GeV and $10^{14}$
GeV, respectively. The lifetimes for these two masses respecting the TA
neutrino limits turn out to be $7\times 10^{25}$ s and $4.8\times 10^{25}$ s,
respectively. Due to the poorer sensitivity of the neutrino detectors, the
$\tau_{\chi}$ limiting values from UHE neutrinos are at least $3-4$ orders of
magnitude smaller than the gamma-ray and UHECR constraints. Thus, the three
UHE messengers give rise to different constraints on the SHDM lifetime and
mass.

In all the three plots of Fig. 1, we do not take into account the contribution
from the conventional astrophysical as well as extragalactic SHDM sources. For
UHECRs, the contribution of the extragalactic SHDM component is found to be
negligible [40, 41]. Whereas, the contribution of the conventional
astrophysical component might be significant in the UHECR spectral energy
range below the energy of the TA Amaterasu event, i.e., $\lesssim 10^{11}$ GeV
[41]. However, considering a dominant fraction of UHECRs from the conventional
astrophysical origin allows only a little fraction of UHECRs from SHDM. This
results in stringent constraints on the production of UHECRs via SHDM decay in
the DM mass range $\lesssim 10^{12}$ GeV and do not yield a conservative
description of the constraints on SHDM lifetime. In order to avoid this, we
consider that the UHECRs in this energy range are solely due to SHDM decay.
This assumption gives rise to a lower bound on the SHDM lifetime. For UHE
gamma-rays, the contribution of the extragalactic component is negligible in
the energy range of our interest in comparison to the Galactic SHDM component.
This is because UHE gamma-rays from extragalactic sources are attenuated
severely due to pair production losses in the cosmic microwave background
(CMB) and the extragalactic radio background during propagation to Earth [4,
41]. For UHE neutrinos, the contribution from extragalactic sources can be
important as unlike gamma-rays, neutrinos do not suffer any propagation
losses. Inclusion of extragalactic sources will result in lowering of the UHE
neutrino flux from SHDM decay in the Milky Way which, in turn, will give rise
to larger SHDM lifetime. In order to be consistent with the assumption for
UHECRs above, we do not include the extragalactic neutrino flux in our
analysis. Due to these assumptions, the constraints obtained on the SHDM
lifetime from different messengers can be comprehended as the most
conservative constraints. In the following, we discuss these constraints for
present and future telescopes in detail and expand the description to
different decay channels.

##  III Multi-messenger constraints on SHDM lifetime

Figure 2: Constraints on SHDM decay rate for all possible two-body decay
channels into SM particles derived from present UHECR data and UHE gamma-ray
non-observation limits. Constraints from UHE neutrino non-observation limits
are weak, hence not shown. The shaded regions above the curves in all the
plots represent the excluded part of the parameter space. The blue curves show
the constraints from the TA Amaterasu event together with the low energy UHECR
data. The magenta (TA) and red (PAO) curves depict the constraints from non-
observation of UHE gamma-rays. Constraints from different decay modes are
distinct. The gamma-ray constraints are in general stronger than that of
UHECRs. The vertically shaded grey regions depict the excluded range of
$m_{\chi}$ for the TA Amaterasu event.

In this section, we analyse the possibility that the TA Amaterasu event was
produced via the decay of SHDM in the Milky Way. Following this, we place
constraints on the lifetime/decay rate of SHDM. In order to derive the
constraints in the energy band lower to the detected energy of the TA
Amaterasu event, we consider the UHECR data detected by the TA. We also carry
out a multi-messenger analysis on the SHDM decay rate considering the non-
detection of UHE gamma-rays and neutrinos at various telescopes such as TA and
PAO [84, 49]. Furthermore, we look into the SHDM lifetime probing capabilities
of the various proposed future UHE telescopes. Below we present a brief
discussion on the observational capabilities of different past and present
telescopes as well as the capabilities of various future telescopes for
detection of UHECRs, UHE gamma-rays and UHE neutrinos.

UHECR telescopes and measurements: UHECRs were discovered by the Akeno Giant
Air Shower Array (AGASA) [85] and the High Resolution Fly’s Eye (HiRes) [86,
87] telescopes about a couple of decades ago. In fact, HiRes had detected the
most energetic ($320$ EeV) UHECRs till date [88]. Recent measurements of the
UHECR spectrum above $0.1$ EeV by the PAO [47] and TA [89] experiments agree
well with these old measurements except for some disagreements regarding mass
composition and spectral index at extreme energies. In this work, we focus on
the UHECR spectrum measured by TA and also on the recent TA Amaterasu event.
In future, the Giant Radio Array for Neutrino Detection (GRAND) [57] will
complement these existing experiments and may provide us with useful insights
into UHECR source properties.

Figure 3: Projected constraints on SHDM decay rate from future observation of
UHE gamma-rays and UHE neutrinos for all possible two-body decay channels into
SM particles. The shaded regions above the curves depict the parameter space
which can be probed by the future UHE gamma-ray (PAO upgrade: red and GRAND
200k: purple) and UHE neutrino (IceCube-Gen2: cyan and GRAND 200k: green)
telescopes. UHE gamma-ray observations will put stringent constraints on the
SHDM lifetime. The vertically shaded grey regions depict the excluded range of
$m_{\chi}$ for the TA Amaterasu event.

UHE gamma-ray telescopes: There has been no detection of UHE gamma-rays in the
PeV-EeV energy band till date, except for some Galactic sources slightly above
PeV energies [90]. Various experiments such as Haverah Park [91], AGASA [92],
Yakutsk [93, 94], AGASA+Yakutsk [95, 96], TA [53, 54], PAO [84, 49] have
placed different upper limits on the UHE gamma-ray flux. Here, we focus on the
recent upper limits from TA and PAO. In future, the GRAND 200k telescopes and
the proposed upgrade of PAO surface detector (SD) [43] will either detect or
put more stringent constraints on the UHE gamma-ray flux.
UHE neutrino telescopes: Similar to UHE gamma-rays, no UHE neutrinos have yet
been detected in the past and presently running detectors such as Antarctic
Impulsive Transient Antenna (ANITA) [97], IceCube [58], TA [83] and PAO [98].
Most of the theoretical [4, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108,
109] and experimental [110, 98, 111, 112] upper limits of the UHE neutrino
(astrophysical origin) flux are even below the sensitivities of some of the
future telescopes such as Askaryan Radio Array (ARA) [113], Radio Ice
Cherenkov Experiment Observatory (RICE) [114], and Square Kilometer Array
(SKA) [115, 116]. However, the upcoming IceCube-Gen2 [58] and GRAND 200k [117,
57] experiments have significantly larger sensitivities and hold excellent
potential of UHE neutrino detection. For the analysis regarding UHE neutrinos,
we mainly focus on these two future experiments.

###  III.1 Present UHE telescopes

In this sub-section, we present the constraints on SHDM lifetime from
different observations of UHECRs, UHE gamma-rays and UHE neutrinos and analyse
the origin of the TA Amaterasu event. The constraints are obtained by
normalising the flux of these particles with the experimental data or limits
from non-observation. For UHECRs, we consider the spectrum measured by TA
together with the recent TA Amaterasu event as depicted in the left plot of
Fig. 1. Including the UHECR data from other experiments such as PAO, AGASA and
HiRes will yield similar values of DM lifetime as the CR spectral shapes from
these experiments are also very similar to that of TA. For the UHE gamma-rays,
we take the upper limits on the flux from TA [53, 54] and PAO [84, 49].
Whereas for the UHE neutrinos, we consider the upper limits on the flux from
PAO [98] and IceCube [58]. We consider a wide range of
$m_{\chi}=(10^{11}-10^{15})$ GeV and different two-body decay modes. For the
decay modes, we take all the quarks (u, d, s, b, t, c), the gauge bosons (W,
Z, $\gamma$, gluon), and the Higgs boson (h). The lower limit of $m_{\chi}$ is
chosen to be $10^{11}$ GeV as our primary goal is to analyse the origin of the
TA Amaterasu event at $2.44\times 10^{11}$ GeV. On the other hand, the upper
limit of $m_{\chi}=10^{15}$ GeV is taken to be close to the GUT scale. For the
given $m_{\chi}$, we constrain the lifetime such that the UHE secondary
particle observations, limits and sensitivities are satisfied. In particular,
we follow convention of exclusion plots and constrain the inverse of the
lifetime, i.e., $\tau_{\chi}^{-1}$ rather than the lifetime ($\tau_{\chi}$).
The constraints on $\tau_{\chi}^{-1}$ as a function of $m_{\chi}$ for
different decay modes are depicted in Fig. 2. The shaded regions above the
curves depict the excluded part of the parameter space. The vertically shaded
grey regions depict the excluded range of $m_{\chi}$ for the TA Amaterasu
event as discussed in Sec. II.

In Fig. 2, the constraints on $\tau_{\chi}^{-1}$ from TA’s UHECR data, TA’s
gamma-ray flux upper limit and PAO’s gamma-ray flux upper limit are shown by
the blue, magenta and red curves, respectively. As the bounds from non-
observation of UHE neutrinos are orders of magnitude smaller than the bounds
from UHECR and gamma-rays, the neutrino bounds are not shown in Fig. 2. In the
top panel, the left plot shows the bounds on SHDM lifetime obtained from
$\chi\to hh,~{}W^{+}W^{-},~{}ZZ$, the middle plot shows the bounds from
$\chi\to b\bar{b},~{}t\bar{t},~{}c\bar{c}$, and the right plot shows the
bounds from $\chi\to u\bar{u},~{}d\bar{d},~{}s\bar{s}$. Similarly, in the
bottom panel, the left plot displays the bounds from the decay modes $\chi\to
e^{+}e^{-},~{}\mu^{+}\mu^{-},~{}\tau^{+}\tau^{-}$, the middle plot from
$\chi\to\nu_{e}\bar{\nu}_{e},~{}\nu_{\mu}\bar{\nu}_{\mu},~{}\nu_{\tau}\bar{\nu}_{\tau}$,
and the right plot from $\chi\to\gamma\gamma,~{}gg$. The excluded regions of
the $\tau_{\chi}^{-1}$ vs $m_{\chi}$ parameter in the plots correspond to flux
larger than the experimental data or limits and hence disallowed.

There are some general features of all these exclusion plots coming for
different decay channels. The bounds on $\tau_{\chi}^{-1}$ from UHECR data for
$m_{\chi}\lesssim 10^{13}$ GeV mimic the spectral shape of the UHECR spectrum
for all decay modes, whereas above this energy they remain almost constant due
to the constant flux assumption above energy of the TA Amaterasu event. If we
consider the TA Amaterasu event only, the resulting constraints remain almost
constant for all values of $m_{\chi}$ and are same as the constant values of
$\tau_{\chi}^{-1}$ for $m_{\chi}\gtrsim 10^{13}$ GeV.

The different channels also show some distinctions. The constraints from
leptons are in general weaker than that of quarks or the bosons by about an
order of magnitude. Therefore, the bounds on SHDM lifetime from the TA
Amaterasu event can vary depending on the decay mode. However, within a
specific family (quarks, leptons, bosons), the bounds do not differ
significantly.

Regarding the constraints from UHE neutrinos, the present UHE neutrino
telescopes such as IceCube, TA and PAO are not sensitive to the flux
corresponding to $\tau_{\chi}$ inferred by the TA Amaterasu event and
therefore remains consistent with non-detection of UHE neutrino counterpart to
the TA Amaterasu event. Interestingly, the gamma-ray sensitivities of the TA
and PAO are close to the TA Amaterasu constrained SHDM gamma-ray flux.
However, no such correlated gamma-ray flux has been detected by the TA and PAO
collaboration. This results in much stronger constraints from the gamma-ray
upper limits of these telescopes than that of the UHECR data. Nevertheless, it
is noteworthy that at the extreme values of $m_{\chi}$, close to $10^{15}$ GeV
for the quark channels the UHECR (Fig. 2, top middle and right plot)
constraint is the strongest.

In all decay modes, for values of $m_{\chi}$ above $\sim 10^{13}$ GeV, the TA
UHECR constraint curves cross the UHE gamma-ray constraint ones, opening a
tiny parameter space supporting the SHDM origin of UHECRs. However, this
inference is in conflict with the gamma-ray constraints from PAO except for
the top left plot (bosonic decay modes) in Fig. 2. In this case, a small
region above $m_{\chi}>10^{14}$ GeV is still allowed. However, the allowed
range of $m_{\chi}$ for the TA Amaterasu event is found to be between
$\sim(10^{11.8}-10^{13.6})$ GeV. Note the slight differences in the allowed
mass ranges for different decay modes in the vertically shaded grey bands.
Outside this range the SHDM decay can not explain both the TA Amaterasu event
and the UHECR limit. Because, in this $m_{\chi}$ range, the constraints from
the gamma-ray limits of TA and PAO disfavour the SHDM origin of TA Amaterasu
event.

While considering the PAO limits one needs to be careful with the field of
view (FOV) of the PAO (RA: $0~{}\text{to}~{}2\pi$, Dec:
$-\pi/2~{}\text{to}~{}+\pi/4$). Indeed, the possible coordinate of the recent
TA Amaterasu event (RA: $256^{\circ}$, Dec: $16^{\circ}$) is outside the FOV
of PAO and can explain the non-detection of any correlated gamma-ray in PAO.
In this scenario, the gamma-ray non-detection constraints from the TA turn out
to be more crucial than the PAO constraint. Nevertheless, the values of
$\tau_{\chi}$ inferred by the TA Amaterasu event produce UHE gamma-ray flux
larger than the present upper limit from the PAO. Hence, taking into account
these multi-messenger constraints, we argue that the SHDM decay origin of the
TA Amaterasu event is in severe tension.

The constraints extrapolated to higher $m_{\chi}$ values from the TA Amaterasu
event allow a parameter space of
$10^{-30}s^{-1}\lesssim\tau_{\chi}^{-1}\lesssim 10^{-29}~{}\rm s^{-1}$ for
$m_{\chi}\gtrsim 10^{14}$ GeV, for SHDM originated UHECRs at the extreme
energies of the spectrum. However, this UHE secondary flux being a diffuse
flux, the PAO should also observe gamma-rays at similar energies. For
instance, if we consider decay via the $b\bar{b}$ channel, the PAO can be
expected to detect a UHE gamma-ray flux of about $(10^{-3}-10^{-4})$ $\rm
km^{-2}~{}yr^{-1}~{}sr^{-1}$ at $\sim 10^{20}$ GeV. Hence, future detection of
UHECRs at or above this energy will be crucial for constraining UHECRs from
SHDM decay. Finally, we also study the dependence of these constraints on the
choice of Galactic DM density profile and find that the results are nearly
independent of density profile choice (see Appendix A).

###  III.2 Future UHE telescopes

The above analysis shows that the observed UHECR data together with the
experimental upper limits on the UHE gamma-ray and UHE neutrino flux put a
strong constraint on the SHDM origin for the recent TA Amaterasu event.
However, these stringent constraints do not completely rule out the
possibility of a SHDM component in the UHECR spectrum. It might be possible
that the SHDM component of the UHECR flux is smaller than that of the observed
spectrum, eradicating detection in current telescopes and, therefore motivates
us to look at the possibilities of the future telescopes. Detection of such a
small flux requires more sensitive telescopes. The upcoming telescopes such as
GRAND 200k, IceCube-Gen2 and PAO upgrade having larger sensitivities compared
to the existing ones provide us with better opportunities to probe the SHDM
through detection of UHE gamma-rays and UHE neutrinos. In the following, we
discuss the potential of these future telescopes for constraining
$\tau_{\chi}$ and $m_{\chi}$ for the various two-body decay modes as
considered above.

Figure 4: Constraints on DM lifetime/decay rate from both current and future
experiments considering the DM decay via $b\bar{b}$ channel. The blue curve
depicts the constraint from TA’s UHECR data. The cyan and green dashed curves
represent the constraints from UHE neutrino sensitivities of IceCube-Gen2 and
GRAND 200k, respectively. The magenta and red solid curves show the present
constraints from UHE gamma-ray non-observation in TA and PAO, respectively,
whereas, the red and purple dashed curves depict the projected reach of the
upcoming telescopes PAO upgrade and GRAND 200k, respectively. While the
constraints from current observations (blue, magenta, red solid curves)
suggest $\tau_{\chi}\gtrsim 10^{30}$ s, future gamma-ray observations (red and
purple dashed curves) will be able to probe $\tau_{\chi}$ up to $10^{32}$ s.
The vertically shaded grey regions depict the excluded range of $m_{\chi}$ for
the TA Amaterasu event.

The constraints on the $\tau_{\chi}$ for various decay modes are obtained by
normalising the computed fluxes from SHDM at the level of sensitivities
(gamma-ray and neutrino) of the said telescopes. Fig. 3 shows the gamma-ray
and neutrino bounds on the DM lifetime and mass for the different decay
channels. For the gamma-ray constraints, we consider the sensitivities of the
proposed PAO upgrade [43] and the GRAND 200k detector. The neutrino bounds are
based on the UHE neutrino flux sensitivities of the IceCube-Gen 2 and the
GRAND 200k detectors. The red and purple curves in Fig. 3 depict the
constraints from gamma-rays in PAO upgrade and GRAND 200k, respectively. The
constraints from neutrinos in IceCube-Gen2 and GRAND 200k are shown by the
cyan and green curves, respectively. Overall, the constraints from the gamma-
ray sensitivities are about two to three order of magnitudes stringent than
that of the ones from neutrinos. Similar to bounds from the present
experiments, the more stringent limits come from the quarks and gauge bosons
decay channels as depicted in Fig. 3 for future telescopes. For gamma-rays,
the constraints from GRAND 200k are found to be stronger than that of PAO for
all decay modes. For neutrinos, the constraints from quarks and bosons
channels are found to be similar for both the telescopes. However, for the
leptonic channels the constraints from IceCube-Gen2 are stronger than that of
GRAND 200k. This is because the energy spectra of the leptonic channels give
rise to a peak at the energy close to the value of $m_{\chi}$ [42]. The
constraints anticipated from UHE neutrino detection in both the future
telescopes will be comparable to that from current UHECR data for
$m_{\chi}\gtrsim 2\times 10^{12}$ GeV. Therefore, the results from the future
experiments will be crucial to put further and stronger multi-messenger
constraints on the parameter space of SHDM.

##  IV Discussion and Conclusion

In this work, we have analysed the possibility of the recently detected TA
Amaterasu event, at $244$ EeV to have originated from the decay of SHDM in the
Milky Way. The arrival direction of this UHECR being a void in the local
Universe makes it hard to explain its origin from the conventional
astrophysical shock acceleration and opens up possibility of a SHDM origin.
The analysis is carried out with a multi-messenger approach considering UHECR
data from TA and upper limits on UHE gamma-ray and UHE neutrino flux from
different experiments such as TA and PAO. For this, we have estimated the
UHECR, UHE gamma-ray and UHE neutrino flux produced by SHDM decay with the NFW
DM density profile for various SHDM decay channels and a wide range of SHDM
mass. The constraints on these fluxes from the experimental data/limits have
resulted in constraints on the SHDM lifetime or decay rates. While we have
found an allowed SHDM mass range of $\sim(10^{11.8}-10^{13.6})$ GeV for the TA
Amaterasu event, the multi-messenger gamma-ray constraints from TA and PAO are
in severe tension with the SHDM origin for the said TA Amaterasu event.
However, these constraints indicate the possibility of UHECR flux at higher
energies originating from heavier SHDM particles, i.e., $m_{\chi}\gtrsim
10^{14}$ GeV. In order to test the viability of such a scenario one would need
more sensitive future telescopes such as PAO upgrade, IceCube-Gen2 and GRAND
200k.

Therefore, we have also looked into the potential of these future telescopes
to constrain the lifetime and mass of SHDM. All the results of our work are
summarised in Fig. 4 for the decay mode $\chi\to b\bar{b}$. We show the bounds
on $\tau_{\chi}^{-1}$ from UHECRs (blue curve), UHE gamma-rays from TA
(magenta curve), PAO (solid red curve), PAO upgrade (red dashed curve), GRAND
200k (purple dashed curve) and UHE neutrinos from IceCube-Gen2 (cyan dashed
curve), GRAND 200k (green dashed curve). The constraints on lifetime from the
TA UHECR data are found to be of the order of $10^{30}$ s for SHDM mass,
$m_{\chi}\gtrsim 10^{13}$ GeV, whereas at lower masses the lifetime can be
smaller. The constraints from gamma-ray upper limit of the TA data vary within
an order of magnitude and found to be maximum of $\sim 3\times 10^{30}$ s at
around $m_{\chi}\sim 10^{12}$ GeV. The TA constraints (UHECR and UHE gamma-
ray) alone allow for a small parameter space above $10^{14}$ GeV for the SHDM.
However, the gamma-ray upper limit of PAO is about a factor of $3$-$4$ larger
and creates tension for the SHDM origin of the TA Amaterasu event. The future
telescopes will be able to put much stronger constraints ($\gtrsim 10^{32}$ s)
on the SHDM lifetime from gamma-ray observations as depicted by the purple
dashed (GRAND 200k) and red dashed (PAO SD upgrade) curves. However, future
UHE neutrino detectors such as IceCube-Gen2 (cyan dashed) and GRAND 200k
(green dashed) will place constraints similar to the present UHECR derived
limits (blue solid). These results vary substantially between different decay
modes, we have done a detailed analysis for the different decay channels and
found that the lepton channels to have the weakest constraints (see, Fig. 2
and Fig. 3).

To conclude, we have provided an updated estimate of multi-messenger
constraints on the lifetime of SHDM by considering the recent TA Amaterasu
event in conjunction with the TA UHECR data and non-observation of UHE gamma-
rays and UHE neutrinos. While doing so, we have ensured consistency with the
existing constraints in the literature [41, 62, 61, 59, 73, 40, 118]. Note
that the constraints for heavier DM mass, i.e., $m_{\chi}>10^{12}$ GeV are not
very robust due to the lack of data (UHECRs, UHE gamma-rays and UHE neutrinos)
at energies above $\sim 10^{11}$ GeV. For harder UHECR flux anticipated from
SHDM [41] or conventional diffusive acceleration [43, 82], the resulting
constraints on the lifetime for $m_{\chi}\gtrsim 10^{12}$ GeV will be smaller
than the ones obtained in this work. Nevertheless, such shorter lifetimes for
the SHDM are not supported by UHE gamma-ray observations and do not impact the
conclusion of this work.

###### Acknowledgements.

This research (SC and PM) was supported in part by the International Centre
for Theoretical Sciences (ICTS) for participating in the program -
Understanding the Universe Through Neutrinos (code: ICTS/Neus2024/04). The
work of D.B. is supported by the Science and Engineering Research Board
(SERB), Government of India grants MTR/2022/000575 and CRG/2022/000603. D.B.
also acknowledges the support from the Simons Foundation (Award
Number:1023171-RC) to visit the International Institute of Physics, Natal,
Brazil in May 2024 when part of this work was completed. SC and PS have
received funding from SERB projects CRG/2021/002961 and MTR/2021/000540. PS
thanks Toshihiro Fujii for providing crucial information on the TA observation
and Pritam Das for useful discussion. The work of ND is supported by the
Ministry of Education, Government of India via the Prime Minister’s Research
Fellowship (PMRF) December 2021 scheme.

##  Appendix A Dependence on DM density profile

Figure A.1: Top plot: DM energy density vs Galactocentic distance for three
different DM profiles. Bottom plot: Constraints on DM lifetime from UHECR vs
DM mass for different DM profiles considering $b\bar{b}$ channel. The
constraints are nearly independent of the DM profiles.

We have performed the analysis on SHDM decay considering the NFW DM density
profile in the Milky Way throughout the paper. For completeness of the
analysis, it is crucial to understand the variation of the results for
different DM profiles. Therefore, in this appendix, we present an analysis on
the dependence of our results on different DM density profiles. Other than the
NFW profile, here we take two other DM profiles - Einesto [79] and Burkert
[81]. The equations for DM energy density for these two profiles are

$\rho_{\rm Ein}(R_{\rm GC})=\rho_{c}\exp\left[-12.5\left(\left(\frac{R_{\rm
GC}}{R_{C}}\right)^{0.16}-1\right)\right],$ (7)

$\rho_{\rm Bur}(R_{\rm GC})=\frac{\rho_{c}R^{3}_{c}}{(R_{\rm
GC}+R_{c})(R^{2}_{\rm GC}+R^{2}_{c})},$ (8)

respectively.

To illustrate the effect of different DM profiles, we compute the constraints
on SHDM lifetime from TA’s UHECR observation for the DM profiles discussed
above. The results are shown in Fig. A.1. In the top panel of A.1, we show the
DM energy density in the Milky Way as a function of Galactocentric distance
for the different DM profiles. These DM profiles give different
$\mathcal{J}^{\theta_{\rm max}}$ values in equation (5). While for NFW
profile, $\mathcal{J}^{\theta_{\rm max}}$ takes value of $1.99$, for Einesto
and Burkert profiles, its value is $2.08$ and $1.95$ respectively. This
variation of $\mathcal{J}^{\theta_{\rm max}}$ gives a maximum of around $7\%$
changes in SHDM lifetime depending on the DM profile as shown in the bottom
plot of figure A.1. Hence, the dependence of the constraints on the DM profile
is negligible.

## References

  * Abbasi  _et al._ [2023a] R. U. Abbasi _et al._ (Telescope Array), Science 382, abo5095 (2023a), arXiv:2311.14231 [astro-ph.HE] .
  * Greisen [1966] K. Greisen, Phys. Rev. Lett. 16, 748 (1966).
  * Zatsepin and Kuz’min [1966] G. T. Zatsepin and V. A. Kuz’min, Soviet Journal of Experimental and Theoretical Physics Letters 4, 78 (1966).
  * Chakraborty  _et al._ [2024] S. Chakraborty, P. Mehta, and P. Sarmah, JCAP 01, 058 (2024), arXiv:2307.15667 [astro-ph.HE] .
  * Gelmini  _et al._ [2008] G. Gelmini, O. E. Kalashev, and D. V. Semikoz, J. Exp. Theor. Phys. 106, 1061 (2008), arXiv:astro-ph/0506128 .
  * Gelmini  _et al._ [2022] G. B. Gelmini, O. Kalashev, and D. Semikoz, Universe 8, 402 (2022), arXiv:2206.00408 [astro-ph.HE] .
  * Bhattacharjee and Sigl [2000] P. Bhattacharjee and G. Sigl, Phys. Rept. 327, 109 (2000), arXiv:astro-ph/9811011 .
  * Abdollahi  _et al._ [2020] S. Abdollahi _et al._ (Fermi-LAT), Astrophys. J. Suppl. 247, 33 (2020), arXiv:1902.10045 [astro-ph.HE] .
  * Farrar and Gruzinov [2009] G. R. Farrar and A. Gruzinov, Astrophys. J. 693, 329 (2009), arXiv:0802.1074 [astro-ph] .
  * Unger and Farrar [2024] M. Unger and G. R. Farrar, Astrophys. J. Lett. 962, L5 (2024), arXiv:2312.13273 [astro-ph.HE] .
  * Kuznetsov [2024] M. Y. Kuznetsov, JCAP 04, 042 (2024), arXiv:2311.14628 [astro-ph.HE] .
  * Zhang  _et al._ [2024a] B. T. Zhang, K. Murase, N. Ekanger, M. Bhattacharya, and S. Horiuchi, (2024a), arXiv:2405.17409 [astro-ph.HE] .
  * Zou  _et al._ [2023] L. Zou, P. Zhang, and Y. Cho, (2023), arXiv:2311.17412 [hep-ph] .
  * Cho and Cho [2024] Y. M. Cho and F. H. Cho, Phys. Lett. B 851, 138598 (2024), arXiv:2312.08115 [hep-ph] .
  * Perri  _et al._ [2023] D. Perri, K. Bondarenko, M. Doro, and T. Kobayashi, (2023), arXiv:2401.00560 [hep-ph] .
  * Frampton and Kephart [2024] P. H. Frampton and T. W. Kephart, (2024), arXiv:2403.12322 [hep-ph] .
  * Zhang  _et al._ [2024b] C. Zhang, S.-H. Zhang, B. Fu, J.-F. Zhang, and X. Zhang, (2024b), arXiv:2404.04926 [hep-ph] .
  * Lang [2024] R. G. Lang, (2024), arXiv:2405.03528 [astro-ph.HE] .
  * Medina-Tanco and Watson [1999] G. A. Medina-Tanco and A. A. Watson, Astropart. Phys. 12, 25 (1999), arXiv:astro-ph/9903182 .
  * Dubrovich  _et al._ [2004] V. K. Dubrovich, D. Fargion, and M. Y. Khlopov, Astropart. Phys. 22, 183 (2004), arXiv:hep-ph/0312105 .
  * Kalashev  _et al._ [2008] O. E. Kalashev, B. A. Khrenov, P. Klimov, S. Sharakin, and S. V. Troitsky, JCAP 03, 003 (2008), arXiv:0710.1382 [astro-ph] .
  * Kalashev and Kuznetsov [2017] O. E. Kalashev and M. Y. Kuznetsov, JETP Lett. 106, 73 (2017), arXiv:1704.05300 [astro-ph.HE] .
  * Aloisio and Tortorici [2008] R. Aloisio and F. Tortorici, Astropart. Phys. 29, 307 (2008), arXiv:0706.3196 [astro-ph] .
  * Aloisio  _et al._ [2015] R. Aloisio, S. Matarrese, and A. V. Olinto, JCAP 08, 024 (2015), arXiv:1504.01319 [astro-ph.HE] .
  * Kuzmin and Tkachev [1999] V. Kuzmin and I. Tkachev, Phys. Rev. D 59, 123006 (1999), arXiv:hep-ph/9809547 .
  * Ling and Long [2021] S. Ling and A. J. Long, Phys. Rev. D 103, 103532 (2021), arXiv:2101.11621 [astro-ph.CO] .
  * Kramer  _et al._ [2021] E. D. Kramer, E. Kuflik, N. Levi, N. J. Outmezguine, and J. T. Ruderman, Phys. Rev. Lett. 126, 081802 (2021), arXiv:2003.04900 [hep-ph] .
  * Kim and Kuflik [2019] H. Kim and E. Kuflik, Phys. Rev. Lett. 123, 191801 (2019), arXiv:1906.00981 [hep-ph] .
  * Baker  _et al._ [2020] M. J. Baker, J. Kopp, and A. J. Long, Phys. Rev. Lett. 125, 151102 (2020), arXiv:1912.02830 [hep-ph] .
  * Babichev  _et al._ [2020] E. Babichev, D. Gorbunov, S. Ramazanov, and L. Reverberi, JCAP 09, 059 (2020), arXiv:2006.02225 [hep-ph] .
  * Chianese  _et al._ [2020] M. Chianese, B. Fu, and S. F. King, JCAP 06, 019 (2020), arXiv:2003.07366 [hep-ph] .
  * Supanitsky and Medina-Tanco [2019] A. D. Supanitsky and G. Medina-Tanco, JCAP 11, 036 (2019), arXiv:1909.09191 [astro-ph.HE] .
  * Alcantara  _et al._ [2019] E. Alcantara, L. A. Anchordoqui, and J. F. Soriano, Phys. Rev. D 99, 103016 (2019), arXiv:1903.05429 [hep-ph] .
  * Babichev  _et al._ [2019] E. Babichev, D. Gorbunov, and S. Ramazanov, Phys. Lett. B 794, 69 (2019), arXiv:1812.03516 [hep-ph] .
  * Kolb and Long [2017] E. W. Kolb and A. J. Long, Phys. Rev. D 96, 103540 (2017), arXiv:1708.04293 [astro-ph.CO] .
  * Marzola and Urban [2017] L. Marzola and F. R. Urban, Astropart. Phys. 93, 56 (2017), arXiv:1611.07180 [astro-ph.HE] .
  * Kannike  _et al._ [2017] K. Kannike, A. Racioppi, and M. Raidal, Nucl. Phys. B 918, 162 (2017), arXiv:1605.09378 [hep-ph] .
  * Santillan [2012] O. P. Santillan, (2012), arXiv:1202.4182 [hep-th] .
  * Arbuzova  _et al._ [2024] E. V. Arbuzova, A. D. Dolgov, and A. A. Nikitenko, (2024), arXiv:2405.12560 [hep-ph] .
  * Murase and Beacom [2012] K. Murase and J. F. Beacom, JCAP 10, 043 (2012), arXiv:1206.2595 [hep-ph] .
  * Das  _et al._ [2023] S. Das, K. Murase, and T. Fujii, Phys. Rev. D 107, 103013 (2023), arXiv:2302.02993 [astro-ph.HE] .
  * Bauer  _et al._ [2021] C. W. Bauer, N. L. Rodd, and B. R. Webber, JHEP 06, 121 (2021), arXiv:2007.15001 [hep-ph] .
  * Aab  _et al._ [2017] A. Aab _et al._ (Pierre Auger), JCAP 04, 038 (2017), [Erratum: JCAP 03, E02 (2018)], arXiv:1612.07155 [astro-ph.HE] .
  * Abbasi  _et al._ [2019a] R. U. Abbasi _et al._ (Telescope Array), Phys. Rev. D 99, 022002 (2019a), arXiv:1808.03680 [astro-ph.HE] .
  * Abbasi  _et al._ [2021a] R. U. Abbasi _et al._ (Telescope Array), Astrophys. J. 909, 178 (2021a), arXiv:2012.10372 [astro-ph.HE] .
  * Aab  _et al._ [2020] A. Aab _et al._ (Pierre Auger), Phys. Rev. D 102, 062005 (2020), arXiv:2008.06486 [astro-ph.HE] .
  * Abreu  _et al._ [2021] P. Abreu _et al._ (Pierre Auger), Eur. Phys. J. C 81, 966 (2021), arXiv:2109.13400 [astro-ph.HE] .
  * Aab  _et al._ [2016] A. Aab _et al._ (Pierre Auger), (2016), arXiv:1604.03637 [astro-ph.IM] .
  * Abraham  _et al._ [2008] J. Abraham _et al._ (Pierre Auger), Astropart. Phys. 29, 243 (2008), arXiv:0712.1147 [astro-ph] .
  * Abreu  _et al._ [2023] P. Abreu _et al._ (Pierre Auger), Phys. Rev. Lett. 130, 061001 (2023), arXiv:2203.08854 [astro-ph.HE] .
  * Pierre Auger Collaboration [2015] Pierre Auger Collaboration, Nuclear Instruments and Methods in Physics Research A 798, 172 (2015).
  * Teshima [1997] M. Teshima (Telescope Array), in  _32nd Rencontres de Moriond: High-Energy Phenomena in Astrophysics_ (1997) pp. 217–222.
  * Abbasi  _et al._ [2019b] R. U. Abbasi _et al._ (Telescope Array), Astropart. Phys. 110, 8 (2019b), arXiv:1811.03920 [astro-ph.HE] .
  * Abbasi  _et al._ [2021b] R. Abbasi _et al._ (Telescope Array), PoS ICRC2021, 864 (2021b).
  * Dedenko  _et al._ [2013] L. Dedenko, S. Knurenko, A. Makarov, I. Makarov, M. Pravdin, I. Sleptsov, A. Glushkov, G. Fedorova, T. Roganova, and A. Sabourov, in  _33rd International Cosmic Ray Conference_ (2013) p. 0045.
  * Knurenko and Petrov [2018] S. P. Knurenko and I. S. Petrov, JETP Lett. 107, 676 (2018).
  * Álvarez-Muñiz  _et al._ [2020] J. Álvarez-Muñiz _et al._ (GRAND), Sci. China Phys. Mech. Astron. 63, 219501 (2020), arXiv:1810.09994 [astro-ph.HE] .
  * Aartsen  _et al._ [2021] M. G. Aartsen _et al._ (IceCube-Gen2), J. Phys. G 48, 060501 (2021), arXiv:2008.04323 [astro-ph.HE] .
  * Chianese  _et al._ [2021a] M. Chianese, D. F. G. Fiorillo, R. Hajjar, G. Miele, S. Morisi, and N. Saviano, JCAP 05, 074 (2021a), arXiv:2103.03254 [hep-ph] .
  * Bérat  _et al._ [2022] C. Bérat, C. Bleve, O. Deligny, F. Montanet, P. Savina, and Z. Torrès, Astrophys. J. 929, 55 (2022), arXiv:2203.08751 [astro-ph.HE] .
  * Chianese  _et al._ [2021b] M. Chianese, D. F. G. Fiorillo, R. Hajjar, G. Miele, and N. Saviano, JCAP 11, 035 (2021b), arXiv:2108.01678 [hep-ph] .
  * Ishiwata  _et al._ [2020] K. Ishiwata, O. Macias, S. Ando, and M. Arimoto, JCAP 01, 003 (2020), arXiv:1907.11671 [astro-ph.HE] .
  * Argüelles  _et al._ [2023] C. A. Argüelles, D. Delgado, A. Friedlander, A. Kheirandish, I. Safa, A. C. Vincent, and H. White, Phys. Rev. D 108, 123021 (2023), arXiv:2210.01303 [hep-ph] .
  * Esmaili  _et al._ [2012] A. Esmaili, A. Ibarra, and O. L. G. Peres, JCAP 11, 034 (2012), arXiv:1205.5281 [hep-ph] .
  * Munbodh and Profumo [2024] P. Munbodh and S. Profumo, (2024), arXiv:2405.00798 [hep-ph] .
  * Muzio  _et al._ [2019] M. S. Muzio, M. Unger, and G. R. Farrar, Phys. Rev. D 100, 103008 (2019), arXiv:1906.06233 [astro-ph.HE] .
  * Ando and Ishiwata [2015] S. Ando and K. Ishiwata, JCAP 05, 024 (2015), arXiv:1502.02007 [astro-ph.CO] .
  * Cohen  _et al._ [2017] T. Cohen, K. Murase, N. L. Rodd, B. R. Safdi, and Y. Soreq, Phys. Rev. Lett. 119, 021102 (2017), arXiv:1612.05638 [hep-ph] .
  * Blanco and Hooper [2019] C. Blanco and D. Hooper, JCAP 03, 019 (2019), arXiv:1811.05988 [astro-ph.HE] .
  * Cao  _et al._ [2022] Z. Cao _et al._ (LHAASO), Phys. Rev. Lett. 129, 261103 (2022), arXiv:2210.15989 [astro-ph.HE] .
  * Kachelriess  _et al._ [2018] M. Kachelriess, O. E. Kalashev, and M. Y. Kuznetsov, Phys. Rev. D 98, 083016 (2018), arXiv:1805.04500 [astro-ph.HE] .
  * Abbasi  _et al._ [2023b] R. Abbasi _et al._ (IceCube), JCAP 10, 003 (2023b), arXiv:2205.12950 [hep-ex] .
  * Song  _et al._ [2024a] D. Song, K. Murase, and A. Kheirandish, JCAP 03, 024 (2024a), arXiv:2308.00589 [astro-ph.HE] .
  * Song  _et al._ [2024b] D. Song, N. Hiroshima, and K. Murase, (2024b), arXiv:2401.15606 [astro-ph.HE] .
  * Murase and Beacom [2013] K. Murase and J. F. Beacom, JCAP 02, 028 (2013), arXiv:1209.0225 [astro-ph.HE] .
  * Zhang  _et al._ [2010] L. Zhang, C. Weniger, L. Maccione, J. Redondo, and G. Sigl, JCAP 06, 027 (2010), arXiv:0912.4504 [astro-ph.HE] .
  * Navarro  _et al._ [1997] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490, 493 (1997), arXiv:astro-ph/9611107 .
  * Navarro  _et al._ [2004] J. F. Navarro, E. Hayashi, C. Power, A. Jenkins, C. S. Frenk, S. D. M. White, V. Springel, J. Stadel, and T. R. Quinn, Mon. Not. Roy. Astron. Soc. 349, 1039 (2004), arXiv:astro-ph/0311231 .
  * Einasto and Haud [1989] J. Einasto and U. Haud, aap 223, 89 (1989).
  * Graham  _et al._ [2006] A. W. Graham, D. Merritt, B. Moore, J. Diemand, and B. Terzic, Astron. J. 132, 2685 (2006), arXiv:astro-ph/0509417 .
  * Burkert [1995] A. Burkert, apjl 447, L25 (1995), arXiv:astro-ph/9504041 [astro-ph] .
  * Halim  _et al._ [2023] A. A. Halim _et al._ (Pierre Auger), JCAP 05, 024 (2023), arXiv:2211.02857 [astro-ph.HE] .
  * Abbasi  _et al._ [2020] R. U. Abbasi _et al._ (Telescope Array), J. Exp. Theor. Phys. 131, 255 (2020), arXiv:1905.03738 [astro-ph.HE] .
  * Rautenberg [2021] J. Rautenberg (Pierre Auger), PoS ICRC2019, 398 (2021).
  * Hayashida  _et al._ [2000] N. Hayashida _et al._ , arXiv e-prints , astro-ph/0008102 (2000), arXiv:astro-ph/0008102 [astro-ph] .
  * Abbasi  _et al._ [2004] R. U. Abbasi _et al._ (HiRes), Phys. Rev. Lett. 92, 151101 (2004), arXiv:astro-ph/0208243 .
  * Abbasi  _et al._ [2008] R. U. Abbasi _et al._ (HiRes), Phys. Rev. Lett. 100, 101101 (2008), arXiv:astro-ph/0703099 .
  * Bird  _et al._ [1995] D. J. Bird _et al._ (HIRES), Astrophys. J. 441, 144 (1995), arXiv:astro-ph/9410067 .
  * Abbasi  _et al._ [2023c] R. U. Abbasi _et al._ , Astropart. Phys. 151, 102864 (2023c).
  * Cao  _et al._ [2021] Z. Cao _et al._ , Nature (London) 594, 33 (2021).
  * Ave  _et al._ [2000] M. Ave, J. A. Hinton, R. A. Vazquez, A. A. Watson, and E. Zas, Phys. Rev. Lett. 85, 2244 (2000), arXiv:astro-ph/0007386 .
  * Shinozaki  _et al._ [2002] K. Shinozaki _et al._ , The Astrophysical Journal 571, L117 (2002).
  * Glushkov  _et al._ [2007] A. V. Glushkov, D. S. Gorbunov, I. T. Makarov, M. I. Pravdin, G. I. Rubtsov, I. E. Sleptsov, and S. V. Troitsky, Soviet Journal of Experimental and Theoretical Physics Letters 85, 131 (2007), arXiv:astro-ph/0701245 [astro-ph] .
  * Glushkov  _et al._ [2010] A. V. Glushkov, I. T. Makarov, M. I. Pravdin, I. E. Sleptsov, D. S. Gorbunov, G. I. Rubtsov, and S. V. Troitsky (Yakutsk EAS Array), Phys. Rev. D 82, 041101 (2010).
  * Rubtsov  _et al._ [2006] G. I. Rubtsov, L. G. Dedenko, G. F. Fedorova, E. Y. Fedunin, A. V. Glushkov, D. S. Gorbunov, I. T. Makarov, M. I. Pravdin, T. M. Roganova, I. E. Sleptsov, and S. V. Troitsky, Phys. Rev. D 73, 063009 (2006).
  * Rubtsov  _et al._ [2013] G. I. Rubtsov, M. Fukushima, D. Ivanov, B. Stokes, G. Thomson, and S. V. Troitsky, in  _International Cosmic Ray Conference_ , International Cosmic Ray Conference, Vol. 33 (2013) p. 526.
  * Gorham  _et al._ [2009] P. W. Gorham _et al._ (ANITA), Astropart. Phys. 32, 10 (2009), arXiv:0812.1920 [astro-ph] .
  * Aab  _et al._ [2019] A. Aab _et al._ (Pierre Auger), JCAP 10, 022 (2019), arXiv:1906.07422 [astro-ph.HE] .
  * Kalashev  _et al._ [2002] O. E. Kalashev, V. A. Kuzmin, D. V. Semikoz, and G. Sigl, Phys. Rev. D 66, 063004 (2002), arXiv:hep-ph/0205050 .
  * Kusenko and Postma [2001] A. Kusenko and M. Postma, Phys. Rev. Lett. 86, 1430 (2001), arXiv:hep-ph/0007246 .
  * Mannheim  _et al._ [2001] K. Mannheim, R. J. Protheroe, and J. P. Rachen, Phys. Rev. D 63, 023003 (2001), arXiv:astro-ph/9812398 .
  * Madsen [2019] J. Madsen (IceCube), in  _38th International Symposium on Physics in Collision_ (2019) arXiv:1901.02528 [astro-ph.HE] .
  * Esmaeili  _et al._ [2022] A. Esmaeili, A. Capanema, A. Esmaili, and P. D. Serpico, Phys. Rev. D 106, 123016 (2022), arXiv:2208.06440 [hep-ph] .
  * Li and Waxman [2007] Z. Li and E. Waxman, in  _32nd International Cosmic Ray Conference_, Vol. 8 (2007) p. 262, arXiv:0711.4969 [hep-ph] .
  * Lee [1998] S. Lee, Phys. Rev. D 58, 043004 (1998), arXiv:astro-ph/9604098 .
  * Murase  _et al._ [2014] K. Murase, Y. Inoue, and C. D. Dermer, Phys. Rev. D 90, 023007 (2014), arXiv:1403.4089 [astro-ph.HE] .
  * van Vliet  _et al._ [2019] A. van Vliet, R. Alves Batista, and J. R. Hörandel, Phys. Rev. D 100, 021302 (2019), arXiv:1901.01899 [astro-ph.HE] .
  * Heinze  _et al._ [2019] J. Heinze, A. Fedynitch, D. Boncioli, and W. Winter, Astrophys. J. 873, 88 (2019), arXiv:1901.03338 [astro-ph.HE] .
  * Waxman and Bahcall [1999] E. Waxman and J. N. Bahcall, Phys. Rev. D 59, 023002 (1999), arXiv:hep-ph/9807282 .
  * Aartsen  _et al._ [2018] M. G. Aartsen _et al._ (IceCube), Phys. Rev. D 98, 062003 (2018), arXiv:1807.01820 [astro-ph.HE] .
  * Gorham  _et al._ [2019] P. W. Gorham _et al._ (ANITA), Phys. Rev. D 99, 122001 (2019), arXiv:1902.04005 [astro-ph.HE] .
  * Anker  _et al._ [2020] A. Anker _et al._ , JCAP 03, 053 (2020), arXiv:1909.00840 [astro-ph.IM] .
  * Allison  _et al._ [2016] P. Allison _et al._ (ARA), Phys. Rev. D 93, 082003 (2016), arXiv:1507.08991 [astro-ph.HE] .
  * Kravchenko  _et al._ [2012] I. Kravchenko, S. Hussain, D. Seckel, D. Besson, E. Fensholt, J. Ralston, J. Taylor, K. Ratzlaff, and R. Young, prd 85, 062004 (2012), arXiv:1106.1164 [astro-ph.HE] .
  * Bray  _et al._ [2015] J. Bray, J. Alvarez-Muniz, S. Buitink, R. Dagkesamanskii, R. D. Ekers, H. D. E. Falcke, K. Gayley, T. Huege, C. W. James, M. Mevius, R. Mutel, R. J. Protheroe, O. Scholten, F. Schroeder, R. E. Spencer, and S. ter Veen, in  _Advancing Astrophysics with the Square Kilometre Array (AASKA14)_ (2015) p. 144, arXiv:1408.6069 [astro-ph.IM] .
  * Labate  _et al._ [2021] M. G. Labate, J. Wagg, B. Alachkar, S. Breen, G. Swart, and A. van Es, in  _2021 15th European Conference on Antennas and Propagation (EuCAP)_ (2021) pp. 1–5.
  * Kotera [2021] K. Kotera (GRAND), PoS ICRC2021, 1181 (2021), arXiv:2108.00032 [astro-ph.HE] .
  * Das  _et al._ [2024] S. Das, J. A. Carpio, and K. Murase, (2024), arXiv:2405.06382 [hep-ph] .
  * Berezinsky and Zatsepin [1969] V. S. Berezinsky and G. T. Zatsepin, Phys. Lett. B 28, 423 (1969).