Probing heavy dark matter decays with multi-messenger astrophysical data
Koji Ishiwata Oscar Macias Shin'ichiro Ando Makoto Arimoto
astro-ph.HE astro-ph.CO hep-ph




a]Koji Ishiwata,
b,c]Oscar Macias,
c,b]Shin'ichiro Ando
d]and Makoto Arimoto


[a]Institute for Theoretical Physics, Kanazawa
  University, Kanazawa 920-1192, Japan

[b]Kavli Institute for the Physics and Mathematics of the
  Universe (WPI), University of Tokyo, Kashiwa, Chiba 277-8583, Japan

[c]GRAPPA Institute, University of Amsterdam, 1098 XH
  Amsterdam, The Netherlands

[d]Faculty of Mathematics and Physics, Institute of
  Science and Engineering, Kanazawa University, Kanazawa 920-1192,
  Japan





s.ando@uva.nl
arimoto@se.kanazawa-u.ac.jp
ishiwata@hep.s.kanazawa-u.ac.jp
oscar.macias@ipmu.jp



We set conservative constraints on decaying dark matter
  particles with masses spanning a very wide range (10^4-10^16
  GeV). For this we use multimessenger observations of cosmic-ray (CR)
  protons/antiprotons, electrons/positrons, neutrinos/antineutrinos
  and gamma rays. Focusing on decays into the b̅b channel, we
  simulate the spectra of dark matter yields by using the
  Dokshitzer-Gribov-Lipatov-Altarelli-Parisi equations and the
   package. We then propagate the CRs of dark matter
  origin till Earth by using the state-of-the-art numerical frameworks
  ,  and  for the
  solution of the CR transport equation in the extragalactic, Galactic
  region and the heliosphere, respectively. Conservative limits are
  obtained by requiring that the predicted dark matter spectra at
  Earth be less than the observed CR spectra. Overall, we exclude dark
  matter lifetimes of 10^28 s or shorter for all the masses
  investigated in this work. The most stringent constraints reach
  10^30 s for very heavy dark matter particles with masses in the
  range 10^11–10^14 GeV. 


Probing heavy dark matter decays with multi-messenger astrophysical data
    [
    December 30, 2023
========================================================================



§ INTRODUCTION

    


The standard model of cosmology is very successful in
explaining the history and evolution of the Universe.
Precise observations of the cosmic microwave background (e.g.,
<cit.>), primordial abundances of heavy isotopes
created by Big Bang nucleosynthesis (e.g., <cit.>), and
other measurements at smaller scales indicate that the energy budget
of the Universe is dominated by dark matter (DM) and dark
energy. However, their fundamental nature has not been identified yet.
In the case of DM, many different attempts have been made to unravel
its particle nature. For example, direct detection experiments are
aiming to detect nuclear recoil events caused by a specific type of DM
candidate called weakly interacting massive
particles <cit.>. However, it is difficult to measure a
positive signal with this search method if DM is feebly interacting
with standard model particles or is extremely heavy. It is also
possible that DM has a lifetime that is longer than the age of the
Universe, in which case the DM lifetime cannot be measured at direct
detection experiments or collider searches. These difficulties could
potentially be overcome by detecting yields of DM annihilation or
decay with the use of cosmic ray (CR) detectors. Indirect detection
experiments could thus play a complementary role to other search
strategies in our quest for the discovery of DM particles.


In this paper we search for potential DM signatures in a variety of
archival CR data. We focus on heavy DM candidates whose mass ranges
between ∼ 10^4 and 10^16  GeV assuming a finite DM
lifetime.  Such (ultra)heavy DM was proposed in the literature <cit.>. An interesting candidate is decaying gravitino in supergravity model. The CRs from decaying dark matter with TeV scale mass have been studied (see e.g., Refs. <cit.> for earlier works), and recently Ref. <cit.> have extended the study for heavier gravitino whose mass is around EeV. When the DM mass is much larger than ∼ 1 TeV, various
particles are produced as the result of fragmentation processes,
including electroweak cascades. This leads to the production of stable
particles such as p, p̅, γ, e^±, ν and ν̅ that in turn diffuse out from
their sources to our detectors. While propagating, CRs undergo several
interactions in the Galactic and extragalactic regions. For example,
Galactic CRs interact with the interstellar gas, ambient photons and
magnetic fields in the interstellar medium.  In addition,
extragalactic CR protons and anti-protons (photons, electrons and
positrons) experience additional photo-hadronic processes
(electromagnetic cascades) by interacting with the background photon
fields, including the cosmic microwave background (CMB) and
extragalactic background light (EBL).  It will be shown that each CR
species from DM has a characteristic spectrum in the energy range of
10^-3 to 10^16  GeV that could in principle be detected
in archival CR data. There are some works that have a similar aim to our current study (see e.g.,
Refs. <cit.>). However, to the best of our
knowledge, self-consistent simulations of the propagation of all the
stable particles in the energy range of 10^-3 to 10^16 
  GeV in both the Galactic and extragalactic regions have not been
attempted yet.






Here we simulate the production and propagation of DM decay yields,
including p, p̅, γ, e^±,
ν and ν̅, in the Galactic and extragalactic
regions.  Various types of CRs have been observed in a wide energy
range; MeV–TeV γ, p̅ and e^+ with
Fermi-LAT <cit.> and
AMS-02 <cit.>, respectively;
in the PeV energy range, photons are observed or
constrained with, e.g., KASCADE <cit.>,
KASCADE-Grande <cit.>,
CASA-MIA <cit.>,
CASA-BLANCA <cit.>, and DICE <cit.>.
Furthermore, for energies in the EeV range, photon flux upper limits
have been obtained by (PAO) <cit.> and
Telescope Array
(TA) <cit.>. Astrophysical ν have been observed/constrained by
IceCube <cit.>, Pierre Auger
Observatory (PAO) <cit.>, and ANITA <cit.>.
The unprecedented high quality of the publicly available
multi-messenger data described above will allow us to impose robust
constraints on the DM lifetime in a very wide DM mass range. A list
with the CR particles assumed in our analysis is given in
Table <ref> along with the corresponding references that we
use to extract the data.

This paper is organized as follows. Section <ref> presents the computation of the DM decaying spectra for the different CR species and the model frameworks for the solution of the CR transport equation in the extragalactic, Galactic region and Heliosphere, respectively. In Sec. <ref> we show the predicted CR spectra after propagation and the resulting limits on the DM lifetime. Finally, we conclude in Sec. <ref>.     












§ COSMIC RAYS FROM HEAVY DARK MATTER

    


The predicted CR spectrum from DM decays (at source) is given by the
product of two factors: one that encapsulates the particle physics
properties of the DM candidate and another that gives account of the
abundance and distribution of the DM. This is written as

    dΦ_X(E_X,ψ)/dE_X =
      ( 1/4πτ_ dm m_ dmdN_X/dE_X)
      (1/ΔΩ∫_ΔΩ dΩ∫_ l.o.s
      dl ρ_ dm(r(l,ψ))),

where m_ dm is the DM mass, τ_ dm the DM lifetime, r(l,ψ) is the Galactocentric distance, l and ψ are the distance and direction measured along the line of sight, respectively. dN_X/dE_X is the CR spectrum of stable particle X at source, with X=p, p̅, e^±, γ, and ν,
ν̅. 

For local DM energy density ρ_ dm, we adopt the spherically
symmetric Navarro-Frenk-White (NFW) profile:

    ρ(r) = ρ_s/(r/r_s)(r/r_s+1)^2,

where we select r_s = 11 kpc, ρ_⊙ = 0.43 GeV/cm^3 and
R_⊙ = 8.34 kpc for the scale radius and local DM density. We
extract these parameters by inspection of Fig. 6 in
Ref. <cit.>.[We have checked
    that the gamma-ray intensity does not change significantly if
    other halo profile is adopted.  For example, we find a ∼ 10 %
    difference if a Burkert profile <cit.> is used. ]

For definiteness, here we consider a scenario in which DM
decays into b b̅ final states with a branching ratio of 100%.
Our simulations are performed in two steps: First, we compute the CR
spectra at source for p, p̅, e^±, γ and ν,
ν̅ from prompt DM decays. Second, we propagate these
particles in the Galactic and extragalactic medium to derive
observable spectra. A flowchart of our simulations is displayed in
Fig. <ref>.





 §.§ Computation of Cosmic Ray spectra at source





The energy spectra at source of the stable particles resulting from
decaying DM can be computed using the  <cit.>. This is the standard method followed
by most studies in the literature. However, this method can be highly
computationally expensive, specially when the DM mass is
larger than ∼ 10 PeV (in the case of b b̅ final state
particles). Due to this technical limitation, in this work we predict
the CR spectra at source using a hybrid approach.  Namely, we use the
Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations for the
quantum chromodynamics (QCD) calculations involving DM yields with
m_ dm≥ 100 PeV, otherwise we use the  package. The
procedure to solve the DGLAP equations consists of two parts;
calculation of fragmentation functions (FFs) of hadrons h=π^±,
π^0, K^±, K^0, K̅^0, n, n̅ and p,
p̅ (using DGLAP equations); and calculation of the energy
spectra of stable particles resulting from the decay of unstable
hadrons (using ). Similar attempts have been made in earlier
studies using  <cit.>, the
QCD event generator <cit.>,  and FFs in
(supersymmetric) QCD  <cit.>, Monte Carlo simulation and
FFs in (supersymmetric) QCD <cit.>
and FFs in the (supersymmetric) standard
model <cit.>.








The fragmentation functions D^h_i(z,Q^2) of the hadron h for a given 
parton i with energy fraction z are calculated by solving the DGLAP
equations. Currently the next-to-leading order (NLO) results in the
MS scheme are
available in <cit.>, and the
uncertainties of parton distribution functions are provided by <cit.>. In
our study we use the code made available in Refs. <cit.>.

The energy spectra f^I_h(x) of stable particles (I=p, e,
γ, ν) from unstable hadrons h with energy fraction x
are calculated with . Here f^I_h(x) is normalized to
single hadron decay, and both particles and anti-particles are
counted. We have checked that the results agree with the analytical
results given, e.g., in Ref. <cit.>, for pion
decay products.

Consequently, the energy spectra of stable particles from DM
decays (in the bb̅ channel) are given by

    dN_I/dE_I=
      2/m_ dmdN_I/dz ,

where z=2E_I/m_ dm and 

    dN_I/dz=
      2∑_h ∫^1_z dy/y D^h_b(y,m_ dm^2)f^I_h(z/y) .

The factor of 2 included in the right-hand side of
Eq. (<ref>) results from taking into account contributions
from both b and b̅ final states. Figure <ref> shows
the predicted spectra from DM decaying into the b
b̅. For clarity, all panels display the quantity
dN_I/dx=2dN_I/dz|_z=2x (where x=E_I/m_ dm). As it can be
seen, the spectra present asymptotic behavior as the m_ dm
increases. Also, the spectral shape and normalizations are different
for each species under consideration. We have checked that the results
obtained with our hybrid method for m_ dm≤ 10TeV agrees well
with those produced by the PPPC4 <cit.> package with or
without electroweak corrections.[It has been pointed out that
  electroweak corrections become important for DM masses ≳
  10 TeV <cit.>. This effect is essential in the
  simulation of stable particles. Specially for leptophilic decaying
  DM because p, p̅, e^± and γ
  are produced even when DM decays to, for example, neutrino
  pairs. However, for hadronic decays, the electroweak corrections
  have minor effects on the spectra.] For larger DM mass values, our
results for γ can be compared with, e.g.,
Refs. <cit.>. We have found that the
γ spectra shown in the published version of Fig. S12 in
Ref. <cit.> are quantitatively different from our
results.[In private communication with the authors of that
  paper the apparent disagreement has been resolved. They have updated
  their Fig. S12 on the arXiv with results that now match our own.
] On the other hand, the γ spectra shown in Fig. 1 of
Ref. <cit.> are almost consistent with ours. We
noticed that they are harder in the low x regime. To explore this
further, we have compared our results using -only versus those
using DGLAP+ for m_ dm<100 PeV and found that the results
obtained with -only gave softer spectra in the x≲
10^-4 region. Although it would be interesting to run a more
in-depth investigation of this discrepancy from a viewpoint of Monte
Carlo simulations versus DGLAP evolution, this is beyond the scope of
our current study. In addition, it is expected that the CR particles
at such a low x will have a minor effect on the observable CR fluxes
at Earth after propagation.












 §.§ Propagation of Cosmic Rays in the Galaxy



In this section we describe the methods used in this work to model the
propagation of CRs in the Galactic region. The propagation of CR
particles in the Galaxy can be studied with various numerical packages
like  <cit.> or  <cit.>.
However, in the very high energy regime, gamma rays of DM
origin can be attenuated via pair production. As pointed out in
Ref. <cit.>, gamma rays with energies in the range
0.1–100 PeV tend to be absorbed by the the interstellar radiation
field (ISRF) through interactions of the form γγ→
e^+e^-. Since one of the primary science goals of a propagation code
such as  has been the study of lower energy gamma ray
observations with Fermi-LAT, it does not contain specific routines
designed to output attenuated gamma ray maps. In the very high energy
range, we follow the prescription given in
Ref. <cit.> which we implement outside the
 framework, but using the ISRF data that comes with
that package. This is also outlined in the flowchart of
Fig. <ref>.


For CR particles of energies ≲ 10^8 GeV, our method consists
of using the propagation packages [For
  this part of the analysis we use a customized 
  version explained in Ref. <cit.>.] and  for the solution of the transport equation in the interstellar
medium and the heliosphere, respectively. At its core,
 consists of a suite of routines that solve the
particle transport equation via numerical methods. Given a certain CR
source distribution, injection spectrum, boundary conditions and
Galactic structure (e.g. interstellar gas, radiation and magnetic
fields),  makes detailed predictions of relevant
observables for all CR species. The processes accounted for by
 include pure diffusion, convection (Galactic winds),
diffusive re-acceleration (diffusion in energy space), energy losses
(ionization, Coulomb scattering bremmstrahlung, inverse Compton
scattering and synchrotron radiation), nuclear fragmentation, and
radioactive decay <cit.>. Measurements of CR
isotopes and spectra of primary and secondary CR species made by
Voyager 1, PAMELA, AMS-02, BESS and other balloon experiments allow
the estimate of some of the most important CR propagation
parameters. For example, the ratio of the CR halo size to the
diffusion coefficient can be obtained from measurements of stable
secondary particles such as Boron. The resulting degeneracy between
the CR halo size and the diffusion coefficient can be alleviated with
the observed abundances of radioactive isotopes such as ^10_4Be,
^26_13Al, ^36_17Cl and
^54_25Mn <cit.>.

Except for Voyager 1 that since 2012 is streaking through space
outside of the heliosphere, all other indirect or direct CR detectors
reside well within its boundaries. While the 
framework allows for detailed studies of CR propagation through the
Galaxy, it does not contain tools for the solution of the particle
transport equation in the heliosphere. The spectrum of charged CRs
measured at Earth vary with time according to the solar activity. In
particular, solar modulation effects are expected to be important
mainly for CRs of moderate energies (≲ 30–50
GeV).[We note that at this energy level the AMS-02 detector
  has made very precise CR observations which we put to use in the
  present work.]  The  package contains dedicated
routines to robustly model the solar modulation on the Galactic CR
spectra. As Galactic CRs enter the heliosphere their trajectories are
affected by solar wind outflows and corresponding magnetic-field
irregularities.  considers both a macroscopic and small
scale heliospheric magnetic field. The former is given by an
Archimedean spiral and the latter by the irregularities originated in
the solar wind. In particular,  uses Monte Carlo
methods to solve the two-dimensional Parker equation for CR transport
through the heliosphere <cit.>. For rigidities
greater than 1 GV, it assumes a parallel component to the magnetic
field of the diffusion tensor given by:

    K_||=β/3 K_0[ P/1GV
        +g_ low] (1+r/1 AU),

where β=v/c with v the particle velocity c the speed of
light, K_0 is the diffusion parameter, P=qc/|Z|e is the CR
particle rigidity, r is the heliocentric distance from the Sun and,
g_ low represents the level of solar activity. It also assumes
that the perpendicular diffusion coefficient is proportional to
K_||, with their ratio denoted by K_⊥,i/K_||=ρ_i and
i refers to Cartesian coordinate index. 

In order to propagate energetic Galactic CRs to the Earth, we first
use  to obtain the local interstellar spectra (LIS)
and its output is subsequently fed into  which allows
us to calculate modulated CR spectra for the particular time periods
in which the AMS-02 observations were carried out. In
Ref. <cit.> the two packages were combined to
self-consistently model the LIS for protons, helium and anti-protons
assuming different modulation levels and both polarities of the solar
magnetic field. In that work, a propagation parameter scan was carried
out by optimization of a likelihood function constructed using data
taken by AMS-02, BESS, and PAMELA as well as the predicted spectra of
corresponding CR species. Table <ref> displays the
best-fit main  propagation parameters obtained in that
reference which we adopt as our baseline propagation model
setup.[We use the default setting for the
    magnetic fields. It produces synchrotron emissions that can be
    another signal of DM, e.g., <cit.>. Ref. <cit.>
    shows the conservative (progressive) bounds for DM decaying to
    bb̅, τ_ dm≳ 10^24 (10^26)  s for
    m_ dm≥ 10  TeV. It will be seen that this
    constraint (even progressive one) is much weaker that the
    constraints obtained from the gamma-ray observations.  ] These
are the parameters that were found to produce the largest effect on
the propagated CR spectrum. Namely, the CR halo height z_h (in
Galactocentric coordinates), diffusion coefficient D_0 at reference
rigidity R_D=4.5 GV, diffusion slope δ, Alfvén velocity
V_ Alf, convection velocity V_ conv and convection
velocity gradient dV_ conv/dz. Other 
propagation parameters used in our analysis are as given in Tab. 2 and
3 of Ref. <cit.>. In turn, the 
propagation parameter setup assumed in our simulations is displayed in
Tab. <ref>.

Of particular relevance to this study is the production and
propagation of CR p̅ and e^±.  classifies
the p̅ produced by our Galaxy as “secondary” and
“tertiary” depending on their origin. Namely, “secondary p̅” are produced through the inelastic interactions given
by pp, pA, and AA (where A refers to the atomic number of
heavy nuclei) while “tertiary p̅” result from inelastic
scattering of p̅ at propagation. In the case of
e^±,  also considers primary e^-
which are accelerated in CR sources (e.g., supernova remnants) as well as
secondary e^± from the collisions of nuclei with the interstellar media. We use the same parameter
setup shown in Tab. <ref> for the computation of the
fore/background p̅ as well as those of DM origin. Also,
Ref. <cit.> did not include e^± data in their
MCMC scans. In this sense, we do not expect to obtain a suitable
astrophysical background model for e^± using
Tab. <ref>. In light of this we have opted for using
the same propagation setup for p̅ as for e^± when
propagating e^± of DM origin but only simulated an astrophysical
background model for p̅ particles. It will be detailed in a
later section that this is a conservative assumption as in the
e^± case we will compute DM constraints by imposing that our DM
predicted fluxes do not saturate e^± measurements.
 





We note that p and p̅ of energies ≳
10^8 GeV propagate just like neutral particles
and thus we could apply the same propagation methods as for photons
and neutrinos. In this case, we can safely neglect the diffusion
effects. As such, we compute the flux of these CRs at Earth by
computing a line-of-sight integral as is done in
Ref. <cit.>.







 §.§ Propagation of Cosmic Rays in the extragalactic region



Decay products from DM undergo cascading processes in the
extragalactic region during the propagation to Earth. We use
 <cit.> for the
simulation of such processes. Within the  framework,
 <cit.> and
 <cit.> are assumed for the computation of
photo-hadronic processes and electromagnetic cascades,
respectively. We have customized the original code to include CR
particles from decaying DM.  is specially
suitable to study the propagation of CR nuclei, photons and
electrons/positrons.

In the case of p and p̅, two photo-hadronic processes are
relevant (see also `/photo-hadronic' in
Fig. <ref>):


  * Photo-pion production: p+γ_ bg→ p+π,

  * Pair production (Bethe-Heitler): p+γ_ bg→ p+e^++e^-.

Here γ_ bg refers to the background photons present in the
extragalactic region. For this component we take into account the CMB
and EBL (using the default model in Kneiske
2004 <cit.>). Through the two processes mentioned
above, e^±, γ and ν, ν̅ are produced as
secondary CRs. The threshold energies for photo-pion production and
pair production are estimated as 6.8× 10^10 (
  meV/E_γ_ bg)  GeV and 4.8× 10^8 (
  meV/E_γ_ bg)  GeV,
respectively <cit.>, with E_ bg being the energy
of the background photons. For p with energies above ∼
10^11  GeV, photo-pion productions becomes the dominant
dissipation process with mean energy-loss length of
10 Mpc <cit.>. This is the main process of the
Greisen–Zatsepin–Kuzmin (GZK)
effect <cit.>.  Photodisintegration and
elastic scattering processes, on the other hand, are irrelevant for
p, and we have checked that nuclear decays produce negligible
effects on the p propagation.

As for e^± and γ case, four different electromagnetic
cascading effects need to be taken into account (see also
`/EM cascades' in Fig. <ref>),


  * Inverse Compton scattering (ICS): e^±+γ_ bg→
  e^±+γ_ bg,

  * Triplet pair production (TPP): e^±+γ_ bg→ e^± + e^+
  + e^-,

  * Pair production (PP): γ + γ_ bg→ e^+ + e^-,

  * Double pair production (DPP): γ + γ_ bg→ e^+ +
  e^- +e^+ + e^-.

For the photon background fields, we assume the default setting in
: CMB, EBL (Stecker 2006 model <cit.>),
and radio background (Protheroe 1996
model <cit.>).[Sometimes the EBL and radio
  background are called IRB and URB in ,
  respectively.]  The impact of each process can be seen in Fig. 5 of
Ref. <cit.>. Regarding e^±-γ_ bg
scattering, ICS (TPP) with the CMB is dominant for energy of e^±
smaller (larger) than 10^8  GeV. As for γ-γ_
  bg scattering, DPP is subdominant compared to PP.  The later is
most relevant in the energy range of 10^5  GeV≲ E
≲ 10^11  GeV where the main photon background is again
the CMB. It is clear that interactions with the CMB is the most
relevant process in a wide energy range. The
  inter-galactic magnetic fields, on the other hand, have large
  uncertainties. A lower bound is obtained, e.g.,
  Ref. <cit.>, that is round 10^-19  G. On
  the other hand, it is shown in Ref. <cit.> that the
  synchrotron process becomes subdominant when the magnetic fields are
  smaller than 0.1 nG. Therefore, we conservatively ignore the effects
  of the magnetic fields in our evaluation.

Finally, ν and ν̅ are produced via the photo-hadronic
interactions in addition to the prompt DM decay. Such high-energy
neutrinos may suffer from resonant absorption
processes <cit.>. However, we have found that this has
a negligible effect on the neutrino propagation. Therefore, neutrinos
produced via both the photo-hadronic interactions and the prompt dark
matter decay only get redshifted when they reach Earth.





§ RESULTS



















In this section we present our procedure to derive conservative
constraints on the DM lifetime. Except for a few cases detailed below,
we do not subtract background/foreground models and only require that
any putative DM signal does not overshoot the observed CR flux at any
given energy bin. In practice, this means that our lower limits on the
DM lifetime are calculated by varying the dark matter lifetime until
the observed CR flux is saturated. However, when using γ-rays
(Fermi-LAT) and p̅ (AMS-02) data, we will run our lower limits
pipeline by taking into account the respective background/foreground
models. This is because our understating of the astrophysical
background for these particular channels has increasingly improved
recently, and thus, we can be less conservative when using these data
sets.

For the two exceptions mentioned above we use the F-test to compute
the 95% CL lower limits. This is done by comparing the null model
(background-only hypothesis), with the alternate model (background
plus DM hypothesis), where the DM flux norm is fixed to a specific
value. This value is then changed until the difference between
χ^2_ Null and χ^2_ Alternate cannot be explained
by the loss of a degree of freedom within a 95% confidence. In
particular,

    F(l_ fixed,N-l_ Null)   =
      (χ^2_ Alternate/χ^2_ Null-1 )
      N-l_ Null/l_ fixed,

where l_ Null is the number of parameters of the null model,
N is the number of data points and l_ fixed is the difference
of number of parameters between the null and alternate hypotheses. We
solve Eq. (<ref>) with the non-linear least-squares
minimization package
.[<https://lmfit.github.io/lmfit-py/intro.html>]



 §.§ Cosmic ray fluxes




Using Eq. (<ref>) and the propagation methods
explained in the previous section we compute predicted CR fluxes at
Earth originating from DM masses larger than 10^3 GeV. Several
particular examples of our predictions along with their respective
data sets (that will be used to impose constraints) are shown
below. Specifically, here we show observable fluxes for CR p,
p̅, γ-rays, and ν. All the CR spectra shown in
this section assume a DM lifetime of 10^27 s.

Figure <ref> displays the p+p̅ fluxes for DM masses
of 10^10, 10^12, 10^14, and 10^16 GeV (from top to
bottom, left to right). As it can be seen, the Galactic components are
comparable to the extragalactic ones for m_ dm≲
10^11  GeV, however the later become dominant for larger DM
masses. We anticipate that more stringent bounds on DM lifetime will
be obtained by using the predicted Galactic CR spectra.  Furthermore,
while the extragalactic contributions are suppressed for m_
  dm≳ 10^11  GeV, its overall intensity remains
unchanged up to m_ dm∼ 10^11  GeV. This behavior is
a result of the GZK effect.  Namely, p (p̅) lose their
energies due to photo-pion production process which is relevant for
p energy over 10^11  GeV.  Then part of that lost
energy is converted into pions, whose decay products emit a given
amount of γ, e^± and ν, ν̅.  Although their
fluxes are suppressed for E≳ 10^11  GeV, these are
nonetheless comparable to the observed CR fluxes at Earth. Thus,
models of new physics predicting DM particles with τ_
  dm≲ 10^27  s and m_ dm≳ 10^10 
  GeV are expected to be constrained by observations.


We show the p̅ spectra for m_ dm=10^3 and 10^4 GeV
in Fig. <ref>. In this figure, the astrophysical
background is also shown. As explained in the previous section, the
astrophysical background used in this work reproduces the one explored
in Ref. <cit.>.[The antiproton background
  computed with the GALPROP-HelMod method in
  Ref. <cit.> slightly overpredicts the AMS-02
  measurements at  10 GeV. However, no such discrepancy is observed
  when the predictions are compared to PAMELA
  data <cit.>. It should be mentioned that the MCMC
  scan procedure performed in that study included antiproton data from
  PAMELA, AMS-02 and BESS-Polar II. So systematic uncertainties in
  that energy range explain any apparent discrepancy between
  background model predictions and observations. In our work we
  confirm such results and set out to impose constraints on decaying
  DM particles.] In this case we find that the extragalactic flux
spectra is negligibly small for this energy range. In addition, it can
be noticed that the p̅ flux gets suppressed as the DM mass
increases. It will be shown in the next section, that the resulting
constraints for this channel (using AMS-02 data) are stringent around
m_ dm∼ 1  TeV but become weaker for larger DM
masses. Using the same propagation parameter setup as for other CR
species, the e^+ spectra is computed. It turns out that the flux is
much smaller than the AMS-02 e^+ data for τ_
  dm=10^27  s and that it is suppressed when the DM mass
gets large. We have found the constraints from the AMS-02 e^+ data
is irrelevant.



Figure <ref> shows γ fluxes for the same mass values
assumed in Fig. <ref>. The spectral bump seen in the high
energy regime corresponds to the contribution from the Galactic
component. We find that the γ rays due to the ICS and
bremsstrahlung in the Galaxy are subdominant in the total flux. The
extragalactic component, on the other hand, exhibits two spectral
peaks; one at low energies and another one at high energies. The
former originates in the cascades from prompt DM decays, while the
later arises from electromagnetic cascades of γ and e^±
coming from photo-hadronic processes. In all the panels we observe an
energy range (10^5  GeV≲ E ≲ 10^10 
  GeV) where the emission of γ is suppressed. This is because
the PP process is so effective that photons with these energies lose
most of their energy producing lower energy γ and e^±. This
explains how even for very high DM masses a fair amount of photons
with energies of MeV to TeV exist.  We note that this fact makes it
possible to constrain decaying DM particles of very high masses using
Fermi-LAT observations.[We found that there is a few factor
  uncertainties in the γ-ray flux in the Fermi-LAT energy range
  using , which was also stated in
  Ref. <cit.>. As it will be shown later, however, these
  uncertainties are irrelevant for the constraints on the DM
  lifetime.] Furthermore, as can be seen specially in the bottom row
of Fig. <ref>, γ with energies larger than
10^11  GeV also survive. Consequently, the CR fluxes
observed by PAO and TA can be used to constrain such γ fluxes.

Figure <ref> shows the integrated gamma flux. In this
energy range, the flux is dominated by Galactic contributions.
It is seen that the lifetime of DM is expected to be
constrained by CASA-MIA, KASCADE, and KASCADE-Grande for m_ dm≳ 10^9  GeV and by TA and PAO for m_ dm≳
10^12  GeV.

Finally ν+ν̅ fluxes are displayed in
Fig. <ref>. Here the Galactic contributions are shown
separately. As can be seen, the Galactic component is subdominant
compared to the extragalactic one.  As what happened in the photon
channel, neutrino fluxes in the extragalactic region are composed of
two components; prompt neutrinos from DM and secondary ones resulting
from photo-hadronic processes. We find that the secondary neutrinos
contribute much less than the prompt component. We see that the prompt
component starts to surpass observed flux or the upper bounds for DM
masses of 10^6  GeV≲ m_ dm≲ 10^12 
  GeV. As such, this observations (upper limits) can be used to
constrain the DM lifetime in this mass range.




 §.§ Constraints on dark matter lifetime








Using the observational data and our flux predictions, we set
conservative and robust constraints on the DM lifetime as a function
of its mass.  Figure <ref> shows the main results of our
study. To demonstrate the impact of the Galactic and extragalactic CRs
from DM, we construct lower limits on the lifetime by using
both components separately. In that figure, we derive 95% CL limits
from Fermi-LAT and AMS-02 data while the limits from other
observations are given at the CL of each observation as shown in
Table <ref>. Fig. <ref> shows a
combination of extragalactic and Galactic limits together and include
a comparison with previous results in the
literature <cit.>.


The left panel of Fig. <ref> shows lower limits for the
lifetime obtained by using the Galactic fluxes, while right one is
given by the extragalactic fluxes. PAO and KASCADE-Grande give the
most stringent constraints on τ_ dm due to Galactic γ
rays. The PAO data bounds the DM lifetime τ_ dm≳
10^30s for 10^10  GeV≲ m_ dm≲
10^15  GeV. In the mass range 10^8  GeV≲
m_ dm≲ 10^10  GeV, KASCADE-Grande gives the most
stringent constraints, i.e., τ_ dm≳
10^29s. We note that these results are consistent with
Ref. <cit.>. Finally, Fermi-LAT constrains the DM
lifetime at roughly τ_ dm≳ 10^28s for
10^3  GeV≲ m_ dm≲ 10^6  GeV.
However, constrains obtained using p+p̅ and p̅ spectra
and PAO and AMS-02 data, respectively, are found to be weaker than
those obtained using gamma-ray observations. It is also found that the
constraints from e^+ flux data by AMS-02 is so weak that it is out
of the range of the plot.

The constraints obtained using extragalactic CRs are shown in the
right panel of Fig. <ref>.  It turns out that these are
weaker compared to those obtained with Galactic ones for most of the
DM mass range. The exception being the mass range of 10^6 
  GeV≲ m_ dm≲ 10^8  GeV where we find
that the constrains on the neutrino flux using IceCube observations
are the most stringent, i.e., τ_ dm≳
10^28s. This is consistent with limits reported in
Ref. <cit.>. It is worth noticing that Fermi-LAT gives
a constraint on the DM lifetime in the entire DM mass range.  This is
a consequence of cascading processes during the propagation CRs in the
extragalactic region. This is also in agreement with results shown in
Fig. 3 of Ref. <cit.> in m_ dm≤ 10 TeV
obtained through analytic modelling. This is an important consistency
check of our methods given that in this study we simulate CR particles
by using  instead of analytic methods described in that
reference. Reference <cit.> reports a qualitatively
similar result for 10^7  GeV≲ m_ dm≲
10^12  GeV, except that their bound is a factor of a few
weaker. In addition, we find that for 10^13  GeV≲
m_ dm≲ 10^16  GeV the PAO constraints are
comparable to those obtained with Fermi-LAT. Although the constraints
obtained with our extragalactic predictions are found to be weaker
than those using the Galactic component, our simulations could
potentially be used in future analyses of all sky gamma-ray analyses
of, for example, tomographic cross-correlation using the local galaxy
distributions <cit.>.





§ CONCLUSIONS


Using all the multi-messenger astrophysics probes — photons,
protons, anti-protons, and neutrinos, we set constraints on the
lifetime of heavy dark matter particles in the mass ranges between
10^4 and 10^16 GeV.  We computed the fluxes of all the
multi-messenger probes from dark matter decays in both the Galaxy and
extragalactic halos.

The lower limits on heavy dark matter particles that we found are
summarized in Figs. <ref>.  Dark
matter less massive than 10^8 GeV is most stringently constrained by
unresolved diffuse gamma-ray emission measured by Fermi-LAT.  For dark
matter with much heavier masses above ∼10^10 GeV, both gamma
rays and protons of ultrahigh energies measured with Pierre Auger
Observatory are best used to place very stringent lower limits on the
order of 10^30 s.  For masses between 10^8 and 10^10 GeV,
stringent constraints are set with KASCADE-Grande using the predicted
Galactic gamma-ray flux component.

We also found that dark matter decay yields originating in
extragalactic halos produce gamma-ray signals of GeV energies nearly
independent of dark matter mass, and hence, the Fermi-LAT diffuse
gamma-ray background are used to place constraints on the order of
10^28 s throughout the wide mass range between 10^4 and
10^16 GeV.  Yet, in general, the extragalactic constraints are
found to be weaker than those obtained with the Galactic component.
The only exception is the constraints obtained with the IceCube
neutrino data, which provide the best constraints on dark matter decay
in a narrow mass range around 10^7–10^8 GeV.

Overall, we exclude dark matter lifetime (into bb̅ final state)
of 10^28 s or shorter for all the masses investigated in this
work, while the most stringent constraints reach 10^30 s for very
heavy dark matter of 10^11–10^14 GeV.  On the other hand,
studies on decay modes into a final state that involves leptons are
slated for a future study given that the electroweak corrections have
to be carefully assessed, which is a nontrivial problem especially for
dark matter with very heavy masses.

Although the limits derived in this work are comparable with other
existing limits in the literature, our self-consistent simulations
including extragalactic and Galactic propagation effects and all CR
species serve as an important consistency check of previous studies
and at the same time clarifies which components or modelling
assumptions have the greatest impact on the final results.




We are grateful to Daisuke Yonetoku for fruitful discussions in the early
stage of this project, Rafael Alves Batista, Günter Sigl and
Tobias Winchen for useful discussions about the use of
, Shunzo Kumano and Masanori Hirai for providing us
the codes for solving the DGLAP equations and valuable discussions,
Timothy Cohen for private communication regarding
Fig. <ref>, and Kohta Murase. This work was supported by JSPS KAKENHI Grant
Numbers JP17H05402, JP17K14278, JP17H02875, JP18H05542 and Sakigake 2018 Project of Kanazawa University (KI). SA and OM were supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan and by JSPS KAKENHI Grant
Numbers JP17H04836, JP18H04340 and JP18H04578. MA was JSPS KAKENHI Grant Numbers JP17H06362 and 
the JSPS Leading Initiative for Excellent Young Researchers program.




99



Aghanim:2018eyx 
  N. Aghanim et al. [Planck Collaboration],
  
  arXiv:1807.06209 [astro-ph.CO].
  
  

Cyburt:2015mya 
  R. H. Cyburt, B. D. Fields, K. A. Olive and T. H. Yeh,
  
  Rev. Mod. Phys.  88, 015004 (2016)
  doi:10.1103/RevModPhys.88.015004
  [arXiv:1505.01076 [astro-ph.CO]].
  
  
  

Aprile:2018dbl 
  E. Aprile et al. [XENON Collaboration],
  
  Phys. Rev. Lett.  121, no. 11, 111302 (2018)
  doi:10.1103/PhysRevLett.121.111302
  [arXiv:1805.12562 [astro-ph.CO]].
  
  
  
Chung:1998zb 
  D. J. H. Chung, E. W. Kolb and A. Riotto,
  
  Phys. Rev. D 59, 023501 (1998)
  doi:10.1103/PhysRevD.59.023501
  [hep-ph/9802238].
  
  
  
Kuzmin:1998uv 
  V. Kuzmin and I. Tkachev,
  
  JETP Lett.  68, 271 (1998)
  [Pisma Zh. Eksp. Teor. Fiz.  68, 255 (1998)]
  doi:10.1134/1.567858
  [hep-ph/9802304].
  
  
  
Chung:1999ve 
  D. J. H. Chung, E. W. Kolb, A. Riotto and I. I. Tkachev,
  
  Phys. Rev. D 62, 043508 (2000)
  doi:10.1103/PhysRevD.62.043508
  [hep-ph/9910437].
  
  
  
Chung:2001cb 
  D. J. H. Chung, P. Crotty, E. W. Kolb and A. Riotto,
  
  Phys. Rev. D 64, 043503 (2001)
  doi:10.1103/PhysRevD.64.043503
  [hep-ph/0104100].
  
  
 
Kolb:2007vd 
  E. W. Kolb, A. A. Starobinsky and I. I. Tkachev,
  
  JCAP 0707, 005 (2007)
  doi:10.1088/1475-7516/2007/07/005
  [hep-th/0702143].
  
  
 
 Fedderke:2014ura 
  M. A. Fedderke, E. W. Kolb and M. Wyman,
  
  Phys. Rev. D 91, no. 6, 063505 (2015)
  doi:10.1103/PhysRevD.91.063505
  [arXiv:1409.1584 [astro-ph.CO]].
  
  
 
 
  








Takayama:2000uz 
  F. Takayama and M. Yamaguchi,
  
  Phys. Lett. B 485, 388 (2000)
  doi:10.1016/S0370-2693(00)00726-7
  [hep-ph/0005214].
  
  

Ibarra:2007wg 
  A. Ibarra and D. Tran,
  
  Phys. Rev. Lett.  100, 061301 (2008)
  doi:10.1103/PhysRevLett.100.061301
  [arXiv:0709.4593 [astro-ph]].
  
  

Ibarra:2008qg 
  A. Ibarra and D. Tran,
  
  JCAP 0807, 002 (2008)
  doi:10.1088/1475-7516/2008/07/002
  [arXiv:0804.4596 [astro-ph]].
  
  

Ishiwata:2008cu 
  K. Ishiwata, S. Matsumoto and T. Moroi,
  
  Phys. Rev. D 78, 063505 (2008)
  doi:10.1103/PhysRevD.78.063505
  [arXiv:0805.1133 [hep-ph]].
  
  
 
 Covi:2008jy 
  L. Covi, M. Grefe, A. Ibarra and D. Tran,
  
  JCAP 0901, 029 (2009)
  doi:10.1088/1475-7516/2009/01/029
  [arXiv:0809.5030 [hep-ph]].
  
  
 
 Ishiwata:2008cv 
  K. Ishiwata, S. Matsumoto and T. Moroi,
  
  Phys. Lett. B 675, 446 (2009)
  doi:10.1016/j.physletb.2009.04.049
  [arXiv:0811.0250 [hep-ph]].
  
  

Ishiwata:2009vx 
  K. Ishiwata, S. Matsumoto and T. Moroi,
  
  JHEP 0905, 110 (2009)
  doi:10.1088/1126-6708/2009/05/110
  [arXiv:0903.0242 [hep-ph]].
  
  
  
Ishiwata:2009dk 
  K. Ishiwata, S. Matsumoto and T. Moroi,
  
  Phys. Lett. B 679, 1 (2009)
  doi:10.1016/j.physletb.2009.07.004
  [arXiv:0905.4593 [astro-ph.CO]].
  
  
  
Buchmuller:2009xv 
  W. Buchmuller, A. Ibarra, T. Shindou, F. Takayama and D. Tran,
  
  JCAP 0909, 021 (2009)
  doi:10.1088/1475-7516/2009/09/021
  [arXiv:0906.1187 [hep-ph]].
  
  
 
Ishiwata:2010am 
  K. Ishiwata, S. Matsumoto and T. Moroi,
  
  JHEP 1012, 006 (2010)
  doi:10.1007/JHEP12(2010)006
  [arXiv:1008.3636 [hep-ph]].
  
  
  
Dudas:2018npp 
  E. Dudas, T. Gherghetta, K. Kaneta, Y. Mambrini and K. A. Olive,
  
  Phys. Rev. D 98, no. 1, 015030 (2018)
  doi:10.1103/PhysRevD.98.015030
  [arXiv:1805.07342 [hep-ph]].
  
  






































Esmaili:2012us 
  A. Esmaili, A. Ibarra and O. L. G. Peres,
  
  JCAP 1211, 034 (2012)
  doi:10.1088/1475-7516/2012/11/034
  [arXiv:1205.5281 [hep-ph]].
  
  

Murase:2012xs 
  K. Murase and J. F. Beacom,
  
  JCAP 1210, 043 (2012)
  doi:10.1088/1475-7516/2012/10/043
  [arXiv:1206.2595 [hep-ph]].
  
  
  
 Murase:2015gea 
  K. Murase, R. Laha, S. Ando and M. Ahlers,
  
  Phys. Rev. Lett.  115, no. 7, 071301 (2015)
  doi:10.1103/PhysRevLett.115.071301
  [arXiv:1503.04663 [hep-ph]].
  
  
  
Ahlers:2013xia 
  M. Ahlers and K. Murase,
  
  Phys. Rev. D 90, no. 2, 023010 (2014)
  doi:10.1103/PhysRevD.90.023010
  [arXiv:1309.4077 [astro-ph.HE]].
  
  
  
Cohen:2016uyg 
  T. Cohen, K. Murase, N. L. Rodd, B. R. Safdi and Y. Soreq,
  
  Phys. Rev. Lett.  119, no. 2, 021102 (2017)
  doi:10.1103/PhysRevLett.119.021102
  [arXiv:1612.05638 [hep-ph]].
  
  

Aloisio:2015lva 
  R. Aloisio, S. Matarrese and A. V. Olinto,
  
  JCAP 1508, no. 08, 024 (2015)
  doi:10.1088/1475-7516/2015/08/024
  [arXiv:1504.01319 [astro-ph.HE]].
  
  

Kalashev:2016cre 
  O. K. Kalashev and M. Y. Kuznetsov,
  
  Phys. Rev. D 94, no. 6, 063535 (2016)
  doi:10.1103/PhysRevD.94.063535
  [arXiv:1606.07354 [astro-ph.HE]].
  
  

Kachelriess:2018rty 
  M. Kachelriess, O. E. Kalashev and M. Y. Kuznetsov,
  
  Phys. Rev. D 98, no. 8, 083016 (2018)
  doi:10.1103/PhysRevD.98.083016
  [arXiv:1805.04500 [astro-ph.HE]].
  
  

Sui:2018bbh 
  Y. Sui and P. S. Bhupal Dev,
  
  JCAP 1807, no. 07, 020 (2018)
  doi:10.1088/1475-7516/2018/07/020
  [arXiv:1804.04919 [hep-ph]].
  
  

Ackermann:2014usa 
  M. Ackermann et al. [Fermi-LAT Collaboration],
  
  Astrophys. J.  799, 86 (2015)
  doi:10.1088/0004-637X/799/1/86
  [arXiv:1410.3696 [astro-ph.HE]].
  
  

Aguilar:2016kjl 
  M. Aguilar et al. [AMS Collaboration],
  
  Phys. Rev. Lett.  117, no. 9, 091103 (2016).
  doi:10.1103/PhysRevLett.117.091103
  
  

Aguilar:2019owu 
  M. Aguilar et al. [AMS Collaboration],
  
  Phys. Rev. Lett.  122, no. 4, 041102 (2019).
  doi:10.1103/PhysRevLett.122.041102
  
  

Antoni:2005wq 
  T. Antoni et al. [KASCADE Collaboration],
  
  Astropart. Phys.  24, 1 (2005)
  doi:10.1016/j.astropartphys.2005.04.001
  [astro-ph/0505413].
  
  

Apel:2013uni 
  W. D. Apel et al.,
  
  Astropart. Phys.  47, 54 (2013)
  doi:10.1016/j.astropartphys.2013.06.004
  [arXiv:1306.6283 [astro-ph.HE]].
  
  

Apel:2017ocm 
  W. D. Apel et al. [KASCADE Grande Collaboration],
  
  Astrophys. J.  848, no. 1, 1 (2017)
  doi:10.3847/1538-4357/aa8bb7
  [arXiv:1710.02889 [astro-ph.HE]].
  
  
  
Chantell:1997gs 
  M. C. Chantell et al. [CASA-MIA Collaboration],
  
  Phys. Rev. Lett.  79, 1805 (1997)
  doi:10.1103/PhysRevLett.79.1805
  [astro-ph/9705246].
  
  

Glasmacher:1999id 
  M. A. K. Glasmacher et al.,
  
  Astropart. Phys.  10, 291 (1999).
  doi:10.1016/S0927-6505(98)00070-X
  
  

Fowler:2000si 
  J. W. Fowler, L. F. Fortson, C. C. H. Jui, D. B. Kieda, R. A. Ong, C. L. Pryke and P. Sommers,
  
  Astropart. Phys.  15, 49 (2001)
  doi:10.1016/S0927-6505(00)00139-0
  [astro-ph/0003190].
  
  

Swordy:1999um 
  S. P. Swordy and D. B. Kieda,
  
  Astropart. Phys.  13, 137 (2000)
  doi:10.1016/S0927-6505(99)00117-6
  [astro-ph/9909381].
  
  

Aab:2015bza 
  A. Aab et al. [Pierre Auger Collaboration],
  
  arXiv:1509.03732 [astro-ph.HE].
  
  

Aab:2016agp
  A. Aab et al. [Pierre Auger Collaboration],
  
  JCAP 1704 (2017) no.04,  009
  doi:10.1088/1475-7516/2017/04/009
  [arXiv:1612.01517 [astro-ph.HE]].
  
  
  


Abu-Zayyad:2013qwa 
  T. Abu-Zayyad et al. [Telescope Array Collaboration],
  
  Astropart. Phys.  61, 93 (2015)
  doi:10.1016/j.astropartphys.2014.05.002
  [arXiv:1305.7273 [astro-ph.HE]].
  
  
  
Tsunesada:2017aaq 
  Y. Tsunesada, T. Abuzayyad, D. Ivanov, G. Thomson, T. Fujii and D. Ikeda,
  
  PoS ICRC 2017, 535 (2018).
  doi:10.22323/1.301.0535
  
  

Abbasi:2018ywn 
  R. U. Abbasi et al. [Telescope Array Collaboration],
  
  arXiv:1811.03920 [astro-ph.HE].
  









  
Kopper:2017zzm 
  C. Kopper [IceCube Collaboration],
  
  PoS ICRC 2017, 981 (2018).
  doi:10.22323/1.301.0981
  
  










  
Aartsen:2016ngq 
  M. G. Aartsen et al. [IceCube Collaboration],
  
  Phys. Rev. Lett.  117, no. 24, 241101 (2016)
  Erratum: [Phys. Rev. Lett.  119, no. 25, 259902 (2017)]
  doi:10.1103/PhysRevLett.117.241101, 10.1103/PhysRevLett.119.259902
  [arXiv:1607.05886 [astro-ph.HE]].
  
  

  
Aab:2017njo 
  D. Veberič [Pierre Auger Collaboration],
  
  arXiv:1708.06592 [astro-ph.HE].
  
  
  
Gorham:2019guw 
  P. W. Gorham et al. [ANITA Collaboration],
  
  arXiv:1902.04005 [astro-ph.HE].
  
  
  
Karukes:2019jxv 
  E. V. Karukes, M. Benito, F. Iocco, R. Trotta and A. Geringer-Sameth,
  
  arXiv:1901.02463 [astro-ph.GA].
  
  

 Burkert:1995yz 
  A. Burkert,
  
  IAU Symp.  171, 175 (1996)
  [Astrophys. J.  447, L25 (1995)]
  doi:10.1086/309560
  [astro-ph/9504041].
  
  

 Sjostrand:2014zea 
  T. Sjöstrand et al.,
  
  Comput. Phys. Commun.  191, 159 (2015)
  doi:10.1016/j.cpc.2015.01.024
  [arXiv:1410.3012 [hep-ph]].
  
  
  
Bahr:2008pv 
  M. Bahr et al.,
  
  Eur. Phys. J. C 58, 639 (2008)
  doi:10.1140/epjc/s10052-008-0798-9
  [arXiv:0803.0883 [hep-ph]].
  
  

Bellm:2015jjp 
  J. Bellm et al.,
  
  Eur. Phys. J. C 76, no. 4, 196 (2016)
  doi:10.1140/epjc/s10052-016-4018-8
  [arXiv:1512.01178 [hep-ph]].
  
  

Birkel:1998nx 
  M. Birkel and S. Sarkar,
  
  Astropart. Phys.  9, 297 (1998)
  doi:10.1016/S0927-6505(98)00028-0
  [hep-ph/9804285].
  
  

Sarkar:2001se 
  S. Sarkar and R. Toldra,
  
  Nucl. Phys. B 621, 495 (2002)
  doi:10.1016/S0550-3213(01)00565-X
  [hep-ph/0108098].
  
  



Aloisio:2003xj 
  R. Aloisio, V. Berezinsky and M. Kachelriess,
  
  Phys. Rev. D 69, 094023 (2004)
  doi:10.1103/PhysRevD.69.094023
  [hep-ph/0307279].
  
  



Berezinsky:2000up
  V. Berezinsky and M. Kachelriess,
  
  Phys. Rev. D 63 (2001) 034007
  doi:10.1103/PhysRevD.63.034007
  [hep-ph/0009053].
  
  



Barbot:2002ep 
  C. Barbot and M. Drees,
  
  Phys. Lett. B 533, 107 (2002)
  doi:10.1016/S0370-2693(02)01621-0
  [hep-ph/0202072].
  
  
  
Barbot:2002gt 
  C. Barbot and M. Drees,
  
  Astropart. Phys.  20, 5 (2003)
  doi:10.1016/S0927-6505(03)00134-8
  [ ].
  
  
  
  
  
  
Kniehl:2000fe 
  B. A. Kniehl, G. Kramer and B. Potter,
  
  Nucl. Phys. B 582, 514 (2000)
  doi:10.1016/S0550-3213(00)00303-5
  [hep-ph/0010289].
  
  

Kretzer:2000yf 
  S. Kretzer,
  
  Phys. Rev. D 62, 054001 (2000)
  doi:10.1103/PhysRevD.62.054001
  [hep-ph/0003177].
  
  
  
  
Albino:2005me 
  S. Albino, B. A. Kniehl and G. Kramer,
  
  Nucl. Phys. B 725, 181 (2005)
  doi:10.1016/j.nuclphysb.2005.07.010
  [hep-ph/0502188].
  
  

 Pumplin:2001ct 
  J. Pumplin, D. Stump, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai and W. K. Tung,
  
  Phys. Rev. D 65, 014013 (2001)
  doi:10.1103/PhysRevD.65.014013
  [hep-ph/0101032].
  
  

Pumplin:2002vw 
  J. Pumplin, D. R. Stump, J. Huston, H. L. Lai, P. M. Nadolsky and W. K. Tung,
  
  JHEP 0207, 012 (2002)
  doi:10.1088/1126-6708/2002/07/012
  [hep-ph/0201195].
  
  

Martin:2002aw 
  A. D. Martin, R. G. Roberts, W. J. Stirling and R. S. Thorne,
  
  Eur. Phys. J. C 28, 455 (2003)
  doi:10.1140/epjc/s2003-01196-2
  [hep-ph/0211080].
  
  
  
Martin:2003sk 
  A. D. Martin, R. G. Roberts, W. J. Stirling and R. S. Thorne,
  
  Eur. Phys. J. C 35, 325 (2004)
  doi:10.1140/epjc/s2004-01825-2
  [hep-ph/0308087].
  
  

Bluemlein:2002be 
  J. Blumlein and H. Bottcher,
  
  Nucl. Phys. B 636, 225 (2002)
  doi:10.1016/S0550-3213(02)00342-5
  [hep-ph/0203155].
  
  


Hirai:2003pm 
  M. Hirai et al. [Asymmetry Analysis Collaboration],
  
  Phys. Rev. D 69, 054021 (2004)
  doi:10.1103/PhysRevD.69.054021
  [hep-ph/0312112].
  
  

Leader:2005ci 
  E. Leader, A. V. Sidorov and D. B. Stamenov,
  
  Phys. Rev. D 73, 034023 (2006)
  doi:10.1103/PhysRevD.73.034023
  [hep-ph/0512114].
  
  

deFlorian:2005mw 
  D. de Florian, G. A. Navarro and R. Sassot,
  
  Phys. Rev. D 71, 094018 (2005)
  doi:10.1103/PhysRevD.71.094018
  [hep-ph/0504155].
  
  

Hirai:2004wq 
  M. Hirai, S. Kumano and T.-H. Nagai,
  
  Phys. Rev. C 70, 044905 (2004)
  doi:10.1103/PhysRevC.70.044905
  [hep-ph/0404093].
  
  

Hirai:2007cx 
  M. Hirai, S. Kumano, T.-H. Nagai and K. Sudoh,
  
  Phys. Rev. D 75, 094009 (2007)
  doi:10.1103/PhysRevD.75.094009
  [hep-ph/0702250].
  
  
  
Hirai:2011si 
  M. Hirai and S. Kumano,
  
  Comput. Phys. Commun.  183, 1002 (2012)
  doi:10.1016/j.cpc.2011.12.022
  [arXiv:1106.1553 [hep-ph]].
  
  

Kelner:2006tc 
  S. R. Kelner, F. A. Aharonian and V. V. Bugayov,
  
  Phys. Rev. D 74, 034018 (2006)
  Erratum: [Phys. Rev. D 79, 039901 (2009)]
  doi:10.1103/PhysRevD.74.034018, 10.1103/PhysRevD.79.039901
  [astro-ph/0606058].
  
  

Cirelli:2010xx 
  M. Cirelli et al.,
  
  JCAP 1103, 051 (2011)
  Erratum: [JCAP 1210, E01 (2012)]
  doi:10.1088/1475-7516/2012/10/E01, 10.1088/1475-7516/2011/03/051
  [arXiv:1012.4515 [hep-ph]].
  
  
  
Ciafaloni:2010ti 
  P. Ciafaloni, D. Comelli, A. Riotto, F. Sala, A. Strumia and A. Urbano,
  
  JCAP 1103, 019 (2011)
  doi:10.1088/1475-7516/2011/03/019
  [arXiv:1009.0224 [hep-ph]].
  
  

Strong:1998fr 
  A. W. Strong, I. V. Moskalenko and O. Reimer,
  
  Astrophys. J.  537, 763 (2000)
  Erratum: [Astrophys. J.  541, 1109 (2000)]
  doi:10.1086/309038
  [astro-ph/9811296].
  
  

Evoli:2008dv 
  C. Evoli, D. Gaggero, D. Grasso and L. Maccione,
  
  JCAP 0810, 018 (2008)
  Erratum: [JCAP 1604, no. 04, E01 (2016)]
  doi:10.1088/1475-7516/2008/10/018, 10.1088/1475-7516/2016/04/E01
  [arXiv:0807.4730 [astro-ph]].
  
  
  

 Moskalenko:2005ng 
  I. V. Moskalenko, T. A. Porter and A. W. Strong,
  
  Astrophys. J.  640, L155 (2006)
  doi:10.1086/503524
  [astro-ph/0511149].
  
  

Song:2019nrx 
  D. Song, O. Macias and S. Horiuchi,
  
  arXiv:1901.07025 [astro-ph.HE].
  
  

Boschini:2017fxq 
  M. J. Boschini et al.,
  
  Astrophys. J.  840, no. 2, 115 (2017)
  doi:10.3847/1538-4357/aa6e4f
  [arXiv:1704.06337 [astro-ph.HE]].
  
  

Ishiwata:2008qy 
  K. Ishiwata, S. Matsumoto and T. Moroi,
  
  Phys. Rev. D 79, 043527 (2009)
  doi:10.1103/PhysRevD.79.043527
  [arXiv:0811.4492 [astro-ph]].
  
  
  
Crocker:2010xc 
  R. M. Crocker, D. Jones, F. Melia, J. Ott and R. J. Protheroe,
  
  Nature 468, 65 (2010)
  doi:10.1038/nature08635
  [arXiv:1001.1275 [astro-ph.GA]].
  
  

Cirelli:2016mrc 
  M. Cirelli and M. Taoso,
  
  JCAP 1607, 041 (2016)
  doi:10.1088/1475-7516/2016/07/041
  [arXiv:1604.06267 [hep-ph]].
  
  

  
Batista:2016yrx 
  R. Alves Batista et al.,
  
  JCAP 1605, no. 05, 038 (2016)
  doi:10.1088/1475-7516/2016/05/038
  [arXiv:1603.07142 [astro-ph.IM]].
  
  

Heiter:2017cev 
  C. Heiter, D. Kuempel, D. Walz and M. Erdmann,
  
  Astropart. Phys.  102, 39 (2018)
  doi:10.1016/j.astropartphys.2018.05.003
  [arXiv:1710.11406 [astro-ph.IM]].
  
  




  

Mucke:1999yb 
  A. Mucke, R. Engel, J. P. Rachen, R. J. Protheroe and T. Stanev,
  
  Comput. Phys. Commun.  124, 290 (2000)
  doi:10.1016/S0010-4655(99)00446-4
  [astro-ph/9903478].
  
  

Lee:1996fp 
  S. Lee,
  
  Phys. Rev. D 58, 043004 (1998)
  doi:10.1103/PhysRevD.58.043004
  [astro-ph/9604098].
  
  

Kneiske:2003tx 
  T. M. Kneiske, T. Bretz, K. Mannheim and D. H. Hartmann,
  
  Astron. Astrophys.  413, 807 (2004)
  doi:10.1051/0004-6361:20031542
  [astro-ph/0309141].
  
  
  
Stanev:2000fb 
  T. Stanev, R. Engel, A. Mucke, R. J. Protheroe and J. P. Rachen,
  
  Phys. Rev. D 62, 093005 (2000)
  doi:10.1103/PhysRevD.62.093005
  [astro-ph/0003484].
  
  

Greisen:1966jv 
  K. Greisen,
  
  Phys. Rev. Lett.  16, 748 (1966).
  doi:10.1103/PhysRevLett.16.748
  
  

Zatsepin:1966jv 
  G. T. Zatsepin and V. A. Kuzmin,
  
  JETP Lett.  4, 78 (1966)
  [Pisma Zh. Eksp. Teor. Fiz.  4, 114 (1966)].
  
  

Stecker:2005qs 
  F. W. Stecker, M. A. Malkan and S. T. Scully,
  
  Astrophys. J.  648, 774 (2006)
  doi:10.1086/506188
  [astro-ph/0510449].
  
  
  
Protheroe:1996si 
  R. J. Protheroe and P. L. Biermann,
  
  Astropart. Phys.  6, 45 (1996)
  Erratum: [Astropart. Phys.  7, 181 (1997)]
  doi:10.1016/S0927-6505(96)00041-2
  [astro-ph/9605119].
  
  

Finke:2015ona 
  J. D. Finke, L. C. Reyes, M. Georganopoulos, K. Reynolds, M. Ajello, S. J. Fegan and K. McCann,
  
  Astrophys. J.  814, no. 1, 20 (2015)
  doi:10.1088/0004-637X/814/1/20
  [arXiv:1510.02485 [astro-ph.HE]].
  
  

  
Gondolo:1991rn 
  P. Gondolo, G. Gelmini and S. Sarkar,
  
  Nucl. Phys. B 392, 111 (1993)
  doi:10.1016/0550-3213(93)90199-Y
  [hep-ph/9209236].
  
  

  









  









Ando:2015qda 
  S. Ando and K. Ishiwata,
  
  JCAP 1505, no. 05, 024 (2015)
  doi:10.1088/1475-7516/2015/05/024
  [arXiv:1502.02007 [astro-ph.CO]].
  
  
  
Ando:2013xwa 
  S. Ando, A. Benoit-Lévy and E. Komatsu,
  
  Phys. Rev. D 90, no. 2, 023514 (2014)
  doi:10.1103/PhysRevD.90.023514
  [arXiv:1312.4403 [astro-ph.CO]].
  
  

Ando:2014aoa 
  S. Ando,
  
  JCAP 1410, no. 10, 061 (2014)
  doi:10.1088/1475-7516/2014/10/061
  [arXiv:1407.8502 [astro-ph.CO]].
  
  

Fornengo:2013rga 
  N. Fornengo and M. Regis,
  
  Front. Physics 2, 6 (2014)
  doi:10.3389/fphy.2014.00006
  [arXiv:1312.4835 [astro-ph.CO]].
  
  

Ando:2016ang 
  S. Ando and K. Ishiwata,
  
  JCAP 1606, no. 06, 045 (2016)
  doi:10.1088/1475-7516/2016/06/045
  [arXiv:1604.02263 [hep-ph]].
  
  

Shirasaki:2018dkz 
  M. Shirasaki, O. Macias, S. Horiuchi, N. Yoshida, C. H. Lee and A. J. Nishizawa,
  
  Phys. Rev. D 97, no. 12, 123015 (2018)
  doi:10.1103/PhysRevD.97.123015
  [arXiv:1802.10257 [astro-ph.CO]].
  
  

 
Hashimoto:2018ztv 
  D. Hashimoto, A. J. Nishizawa, M. Shirasaki, O. Macias, S. Horiuchi, H. Tashiro and M. Oguri,
  
  doi:10.1093/mnras/stz321
  arXiv:1805.08139 [astro-ph.CO].
  
  





