Probing heavy dark matter decays with multi-messenger astrophysical data

TL;DR

This paper establishes new constraints on decaying heavy dark matter particles across a wide mass range by analyzing multi-messenger astrophysical data, significantly limiting their possible lifetimes based on cosmic-ray and gamma-ray observations.

Contribution

The study introduces a comprehensive multi-messenger approach using advanced simulation and propagation tools to set conservative limits on dark matter decay lifetimes over a broad mass spectrum.

Findings

Dark matter lifetimes shorter than 10^{28} seconds are excluded.

Constraints reach up to 10^{30} seconds for very heavy dark matter.

The analysis covers dark matter masses from 10^4 to 10^{16} GeV.

Abstract

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 $\bar{b}b$ channel, we simulate the spectra of dark matter yields by using the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi equations and the $\texttt{Pythia}$ package. We then propagate the CRs of dark matter origin till Earth by using the state-of-the-art numerical frameworks $\texttt{CRPropa}$, $\texttt{GALPROP}$ and $\texttt{HelMod}$ 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.

Tables

Table 1. Table 1 : Observations of cosmic-ray particles which are used in the analysis. The fourth column shows whether each experiment detected the corresponding CRs. Otherwise, the last column shows the confidence level (CL) of the upper limits quoted in the references.

CRs	Observations	Energy [GeV]	Detected	CL upper limits
Gamma ($\gamma$)	Fermi-LAT [30]	${10}^{-1}$ – ${10}^3$	✓
	CASA-MIA [36]	${10}^5$ – ${10}^7$		90%
	KASCADE [35]	${10}^5$ – ${10}^7$		90%
	KASCADE-Grande [35]	${10}^7$ – ${10}^8$		90%
	PAO [40, 41]	${10}^9$ – ${10}^{10}$		95%
	TA [44]	${10}^9$ – ${10}^{11}$		95%
Proton ($p$)	PAO [47]	${10}^9$ – ${10}^{11}$	✓	84%
Anti-proton ($\overline{p}$)	PAO [47]	${10}^9$ – ${10}^{11}$	✓	84%
	AMS-02 [31]	${10}^{-1}$ – ${10}^2$	✓
Positron ($e^{+}$)	AMS-02 [32]	${10}^{-1}$ – ${10}^3$	✓
Neutrino ($\nu$)	IceCube [45]	${10}^5$ – ${10}^8$	✓	90%
	IceCube [46]	${10}^6$ – ${10}^{11}$		90%
	PAO [47]	${10}^8$ – ${10}^{11}$		90%
	ANITA [48]	${10}^9$ – ${10}^{12}$		90%

Table 2. Table 2 : Main GALPROP propagation parameter setup considered in this study. Our baseline fore/background model corresponds to the best-fit propagation parameter setup obtained in Ref.  [ 81 ] . Other propagation parameters (that impact relatively less the results) are taken from Tab. 2 and 3 of that reference.

$D_0$	$z_h$	$V_{\mathrm{Alf}}$	$\delta$	$V_{\mathrm{convec}}$	$d{V}_{\mathrm{convec}}/dz$
(1028 cm2 s-1)	(kpc)	(km s-1)		(km s-1)	(km s-1)
4.3	4.0	28.6	0.395	12.4	10.2

Table 3. Table 3 : HelMod propagation parameters considered in this study.

${\rho }_i$	$K_0$	$g_{\mathrm{low}}$
	(AU2 GV-1 s-1)
0.065	$3\times {10}^{-5}$	0.4

Equations

$$\frac{d{\Phi }_X(E_X,\psi )}{d{E}_X}=\left(\frac{1}{4\pi {\tau }_{\mathrm{dm}}{m}_{\mathrm{dm}}}\frac{d{N}_X}{d{E}_X}\right)\left(\frac{1}{\Delta \Omega }{\int }_{\Delta \Omega }d\Omega {\int }_{l.o.s}dl{\rho }_{\mathrm{dm}}(r(l,\psi ))\right),$$

$$\rho (r)=\frac{{\rho }_s}{(r/r_s)(r/r_s+1{)}^2},$$

$$\frac{d{N}_I}{d{E}_I}=\frac{2}{m_{\mathrm{dm}}}\frac{d{N}_I}{dz},$$

$$\frac{d{N}_I}{dz}=2\sum _h{\int }_z^1\frac{dy}{y}{D}_b^h(y,m_{\mathrm{dm}}^2)f_h^I(z/y).$$

$$K_{||}=\frac{\beta }{3}{K}_0\left[\frac{P}{1\text{GV}}+g_{\mathrm{low}}\right]\left(1+\frac{r}{\text{1 AU}}\right),$$

$$F(l_{\mathrm{fixed}},N-l_{\mathrm{Null}})=\left(\frac{{\chi }_{\mathrm{Alternate}}^2}{{\chi }_{\mathrm{Null}}^2}-1\right)\frac{N-l_{\mathrm{Null}}}{l_{\mathrm{fixed}}},$$

Full text

1 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., [1]), primordial abundances of heavy isotopes created by Big Bang nucleosynthesis (e.g., [2]), 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 [3]. 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 $\sim 10^{4}$ and $10^{16}\,{\rm GeV}$ assuming a finite DM lifetime. Such (ultra)heavy DM was proposed in the literature [4, 5, 6, 7, 8, 9]. 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. [10, 11, 12, 13, 14, 15, 16, 17, 18, 19] for earlier works), and recently Ref. [20] have extended the study for heavier gravitino whose mass is around EeV. When the DM mass is much larger than $\sim 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$, $\bar{p}$, $\gamma$, $e^{\pm}$, $\nu$ and $\bar{\nu}$ 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}\,{\rm 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. [21, 22, 23, 24, 25, 27, 26, 28, 29]). 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}\,{\rm 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$, $\bar{p}$, $\gamma$, $e^{\pm}$, $\nu$ and $\bar{\nu}$, in the Galactic and extragalactic regions. Various types of CRs have been observed in a wide energy range; MeV–TeV $\gamma$, $\bar{p}$ and $e^{+}$ with Fermi-LAT [30] and AMS-02 [31, 32], respectively; in the PeV energy range, photons are observed or constrained with, e.g., KASCADE [33], KASCADE-Grande [34, 35], CASA-MIA [36, 37], CASA-BLANCA [38], and DICE [39]. Furthermore, for energies in the EeV range, photon flux upper limits have been obtained by (PAO) [40, 41] and Telescope Array (TA) [42, 43, 44]. Astrophysical $\nu$ have been observed/constrained by IceCube [45, 46], Pierre Auger Observatory (PAO) [47], and ANITA [48]. 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 1 along with the corresponding references that we use to extract the data.

This paper is organized as follows. Section 2 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. 3 we show the predicted CR spectra after propagation and the resulting limits on the DM lifetime. Finally, we conclude in Sec. 4.

2 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

[TABLE]

where $m_{\rm dm}$ is the DM mass, $\tau_{\rm dm}$ the DM lifetime, $r(l,\psi)$ is the Galactocentric distance, $l$ and $\psi$ 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$, $\bar{p}$, $e^{\pm}$, $\gamma$, and $\nu$, $\bar{\nu}$.

For local DM energy density $\rho_{\rm dm}$, we adopt the spherically symmetric Navarro-Frenk-White (NFW) profile:

[TABLE]

where we select $r_{s}=11$ kpc, $\rho_{\odot}=0.43$ GeV/cm3 and $R_{\odot}=8.34$ kpc for the scale radius and local DM density. We extract these parameters by inspection of Fig. 6 in Ref. [49].111We have checked that the gamma-ray intensity does not change significantly if other halo profile is adopted. For example, we find a $\sim 10$ % difference if a Burkert profile [50] is used.

For definiteness, here we consider a scenario in which DM decays into $b\bar{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$, $\bar{p}$, $e^{\pm}$, $\gamma$ and $\nu$, $\bar{\nu}$ 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. 1.

2.1 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 Pythia 8.2 [51]. 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 $\sim 10$ PeV (in the case of $b\bar{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_{\rm dm}\geq 100$ PeV, otherwise we use the Pythia package. The procedure to solve the DGLAP equations consists of two parts; calculation of fragmentation functions (FFs) of hadrons $h=\pi^{\pm}$, $\pi^{0}$, $K^{\pm}$, $K^{0}$, $\bar{K}^{0}$, $n$, $\bar{n}$ and $p$, $\bar{p}$ (using DGLAP equations); and calculation of the energy spectra of stable particles resulting from the decay of unstable hadrons (using Pythia). Similar attempts have been made in earlier studies using HERWIG [52, 53], the QCD event generator [54], HERWIG and FFs in (supersymmetric) QCD  [55], Monte Carlo simulation and FFs in (supersymmetric) QCD [57, 56] and FFs in the (supersymmetric) standard model [58, 59].

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 $\overline{\rm MS}$ scheme are available in [60, 61, 62], and the uncertainties of parton distribution functions are provided by [63, 64, 65, 66, 67, 68, 69, 70, 71]. In our study we use the code made available in Refs. [72, 73].

The energy spectra $f^{I}_{h}(x)$ of stable particles ($I=p$, $e$, $\gamma$, $\nu$) from unstable hadrons $h$ with energy fraction $x$ are calculated with Pythia. 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. [74], for pion decay products.

Consequently, the energy spectra of stable particles from DM decays (in the $b\bar{b}$ channel) are given by

[TABLE]

where $z=2E_{I}/m_{\rm dm}$ and

[TABLE]

The factor of 2 included in the right-hand side of Eq. (2.4) results from taking into account contributions from both $b$ and $\bar{b}$ final states. Figure 2 shows the predicted spectra from DM decaying into the $b\bar{b}$. For clarity, all panels display the quantity $dN_{I}/dx=2dN_{I}/dz|_{z=2x}$ (where $x=E_{I}/m_{\rm dm}$). As it can be seen, the spectra present asymptotic behavior as the $m_{\rm 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_{\rm dm}\leq 10~{}$TeV agrees well with those produced by the PPPC4 [75] package with or without electroweak corrections.222It has been pointed out that electroweak corrections become important for DM masses $\gtrsim 10$ TeV [76]. This effect is essential in the simulation of stable particles. Specially for leptophilic decaying DM because $p$, $\bar{p}$, $e^{\pm}$ and $\gamma$ 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 $\gamma$ can be compared with, e.g., Refs. [25, 27]. We have found that the $\gamma$ spectra shown in the published version of Fig. S12 in Ref. [25] are quantitatively different from our results.333In 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 $\gamma$ spectra shown in Fig. 1 of Ref. [27] 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 Pythia-only versus those using DGLAP+Pythia for $m_{\rm dm}<100$ PeV and found that the results obtained with Pythia-only gave softer spectra in the $x\lesssim 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.

2.2 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 GALPROP [77] or DRAGON [78]. However, in the very high energy regime, gamma rays of DM origin can be attenuated via pair production. As pointed out in Ref. [79], 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 $\gamma\gamma\to e^{+}e^{-}$. Since one of the primary science goals of a propagation code such as GALPROP 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. [79] which we implement outside the GALPROP framework, but using the ISRF data that comes with that package. This is also outlined in the flowchart of Fig. 1.

For CR particles of energies $\lesssim 10^{8}$ GeV, our method consists of using the propagation packages GALPROP v54444For this part of the analysis we use a customized GALPROP version explained in Ref. [80]. and HelMod v4.0 for the solution of the transport equation in the interstellar medium and the heliosphere, respectively. At its core, GALPROP 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), GALPROP makes detailed predictions of relevant observables for all CR species. The processes accounted for by GALPROP 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 [79]. 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}_{4}$Be, ${}^{26}_{13}$Al, ${}^{36}_{17}$Cl and ${}^{54}_{25}$Mn [79].

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 GALPROP 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 ($\lesssim 30$–50 GeV).555We 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 HelMod 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. HelMod 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, HelMod uses Monte Carlo methods to solve the two-dimensional Parker equation for CR transport through the heliosphere [81]. For rigidities greater than 1 GV, it assumes a parallel component to the magnetic field of the diffusion tensor given by:

[TABLE]

where $\beta=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_{\rm 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_{\perp,i}/K_{||}=\rho_{i}$ and $i$ refers to Cartesian coordinate index.

In order to propagate energetic Galactic CRs to the Earth, we first use GALPROP to obtain the local interstellar spectra (LIS) and its output is subsequently fed into HelMod which allows us to calculate modulated CR spectra for the particular time periods in which the AMS-02 observations were carried out. In Ref. [81] 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 2 displays the best-fit main GALPROP propagation parameters obtained in that reference which we adopt as our baseline propagation model setup.666We use the default setting for the magnetic fields. It produces synchrotron emissions that can be another signal of DM, e.g., [82, 83, 84]. Ref. [84] shows the conservative (progressive) bounds for DM decaying to $b\bar{b}$, $\tau_{\rm dm}\gtrsim 10^{24}\,(10^{26})\,{\rm s}$ for $m_{\rm dm}\geq 10\,{\rm 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 $\delta$, Alfvén velocity $V_{\rm Alf}$, convection velocity $V_{\rm conv}$ and convection velocity gradient $dV_{\rm conv}/dz$. Other GALPROP propagation parameters used in our analysis are as given in Tab. 2 and 3 of Ref. [81]. In turn, the HelMod propagation parameter setup assumed in our simulations is displayed in Tab. 3.

Of particular relevance to this study is the production and propagation of CR $\bar{p}$ and $e^{\pm}$. GALPROP classifies the $\bar{p}$ produced by our Galaxy as “secondary” and “tertiary” depending on their origin. Namely, “secondary $\bar{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 $\bar{p}$” result from inelastic scattering of $\bar{p}$ at propagation. In the case of $e^{\pm}$, GALPROP also considers primary $e^{-}$ which are accelerated in CR sources (e.g., supernova remnants) as well as secondary $e^{\pm}$ from the collisions of nuclei with the interstellar media. We use the same parameter setup shown in Tab. 2 for the computation of the fore/background $\bar{p}$ as well as those of DM origin. Also, Ref. [81] did not include $e^{\pm}$ data in their MCMC scans. In this sense, we do not expect to obtain a suitable astrophysical background model for $e^{\pm}$ using Tab. 2. In light of this we have opted for using the same propagation setup for $\bar{p}$ as for $e^{\pm}$ when propagating $e^{\pm}$ of DM origin but only simulated an astrophysical background model for $\bar{p}$ particles. It will be detailed in a later section that this is a conservative assumption as in the $e^{\pm}$ case we will compute DM constraints by imposing that our DM predicted fluxes do not saturate $e^{\pm}$ measurements.

We note that $p$ and $\bar{p}$ of energies $\gtrsim 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. [26].

2.3 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 CRPropa 3.1 [85, 86] for the simulation of such processes. Within the CRPropa framework, SOPHIA [87] and DINT [88] 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. CRPropa is specially suitable to study the propagation of CR nuclei, photons and electrons/positrons.

In the case of $p$ and $\bar{p}$, two photo-hadronic processes are relevant (see also ‘CRPropa/photo-hadronic’ in Fig. 1):

• •

Photo-pion production: $p+\gamma_{\rm bg}\to p+\pi$,

• •

Pair production (Bethe-Heitler): $p+\gamma_{\rm bg}\to p+e^{+}+e^{-}$.

Here $\gamma_{\rm 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 [89]). Through the two processes mentioned above, $e^{\pm}$, $\gamma$ and $\nu$, $\bar{\nu}$ are produced as secondary CRs. The threshold energies for photo-pion production and pair production are estimated as $6.8\times 10^{10}({\rm meV}/E_{\gamma_{\rm bg}})~{}{\rm GeV}$ and $4.8\times 10^{8}({\rm meV}/E_{\gamma_{\rm bg}})~{}{\rm GeV}$, respectively [86], with $E_{\rm bg}$ being the energy of the background photons. For $p$ with energies above $\sim 10^{11}\,{\rm GeV}$, photo-pion productions becomes the dominant dissipation process with mean energy-loss length of 10 Mpc [90]. This is the main process of the Greisen–Zatsepin–Kuzmin (GZK) effect [91, 92]. 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^{\pm}$ and $\gamma$ case, four different electromagnetic cascading effects need to be taken into account (see also ‘CRPropa/EM cascades’ in Fig. 1),

• •

Inverse Compton scattering (ICS): $e^{\pm}+\gamma_{\rm bg}\to e^{\pm}+\gamma_{\rm bg}$,

• •

Triplet pair production (TPP): $e^{\pm}+\gamma_{\rm bg}\to e^{\pm}+e^{+}+e^{-}$,

• •

Pair production (PP): $\gamma+\gamma_{\rm bg}\to e^{+}+e^{-}$,

• •

Double pair production (DPP): $\gamma+\gamma_{\rm bg}\to e^{+}+e^{-}+e^{+}+e^{-}$.

For the photon background fields, we assume the default setting in DINT: CMB, EBL (Stecker 2006 model [93]), and radio background (Protheroe 1996 model [94]).777Sometimes the EBL and radio background are called IRB and URB in CRPropa, respectively. The impact of each process can be seen in Fig. 5 of Ref. [86]. Regarding $e^{\pm}$-$\gamma_{\rm bg}$ scattering, ICS (TPP) with the CMB is dominant for energy of $e^{\pm}$ smaller (larger) than $10^{8}~{}{\rm GeV}$. As for $\gamma$-$\gamma_{\rm bg}$ scattering, DPP is subdominant compared to PP. The later is most relevant in the energy range of $10^{5}\,{\rm GeV}\lesssim E\lesssim 10^{11}\,{\rm 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. [95], that is round $10^{-19}\,{\rm G}$. On the other hand, it is shown in Ref. [86] 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, $\nu$ and $\bar{\nu}$ are produced via the photo-hadronic interactions in addition to the prompt DM decay. Such high-energy neutrinos may suffer from resonant absorption processes [96]. 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.

3 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 $\gamma$-rays (Fermi-LAT) and $\bar{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 $\chi^{2}_{\rm Null}$ and $\chi^{2}_{\rm Alternate}$ cannot be explained by the loss of a degree of freedom within a 95% confidence. In particular,

[TABLE]

where $l_{\rm Null}$ is the number of parameters of the null model, $N$ is the number of data points and $l_{\rm fixed}$ is the difference of number of parameters between the null and alternate hypotheses. We solve Eq. (3.1) with the non-linear least-squares minimization package lmfit.888https://lmfit.github.io/lmfit-py/intro.html

3.1 Cosmic ray fluxes

Using Eq. (2.1) 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$, $\bar{p}$, $\gamma$-rays, and $\nu$. All the CR spectra shown in this section assume a DM lifetime of $10^{27}$ s.

Figure 3 displays the $p+\bar{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_{\rm dm}\lesssim 10^{11}\,{\rm 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_{\rm dm}\gtrsim 10^{11}\,{\rm GeV}$, its overall intensity remains unchanged up to $m_{\rm dm}\sim 10^{11}\,{\rm GeV}$. This behavior is a result of the GZK effect. Namely, $p$ ($\bar{p}$) lose their energies due to photo-pion production process which is relevant for $p$ energy over $10^{11}\,{\rm GeV}$. Then part of that lost energy is converted into pions, whose decay products emit a given amount of $\gamma$, $e^{\pm}$ and $\nu$, $\bar{\nu}$. Although their fluxes are suppressed for $E\gtrsim 10^{11}\,{\rm GeV}$, these are nonetheless comparable to the observed CR fluxes at Earth. Thus, models of new physics predicting DM particles with $\tau_{\rm dm}\lesssim 10^{27}\,{\rm s}$ and $m_{\rm dm}\gtrsim 10^{10}\,{\rm GeV}$ are expected to be constrained by observations.

We show the $\bar{p}$ spectra for $m_{\rm dm}=10^{3}$ and $10^{4}$ GeV in Fig. 4. 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. [81].999The antiproton background computed with the GALPROP-HelMod method in Ref. [81] slightly overpredicts the AMS-02 measurements at  10 GeV. However, no such discrepancy is observed when the predictions are compared to PAMELA data [81]. 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 $\bar{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_{\rm dm}\sim 1\,{\rm 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 $\tau_{\rm dm}=10^{27}\,{\rm 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 5 shows $\gamma$ fluxes for the same mass values assumed in Fig. 3. The spectral bump seen in the high energy regime corresponds to the contribution from the Galactic component. We find that the $\gamma$ 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 $\gamma$ and $e^{\pm}$ coming from photo-hadronic processes. In all the panels we observe an energy range ($10^{5}\,{\rm GeV}\lesssim E\lesssim 10^{10}\,{\rm GeV}$) where the emission of $\gamma$ is suppressed. This is because the PP process is so effective that photons with these energies lose most of their energy producing lower energy $\gamma$ and $e^{\pm}$. 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.101010We found that there is a few factor uncertainties in the $\gamma$-ray flux in the Fermi-LAT energy range using DINT, which was also stated in Ref. [88]. 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. 5, $\gamma$ with energies larger than $10^{11}\,{\rm GeV}$ also survive. Consequently, the CR fluxes observed by PAO and TA can be used to constrain such $\gamma$ fluxes.

Figure 6 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_{\rm dm}\gtrsim 10^{9}\,{\rm GeV}$ and by TA and PAO for $m_{\rm dm}\gtrsim 10^{12}\,{\rm GeV}$.

Finally $\nu+\bar{\nu}$ fluxes are displayed in Fig. 7. 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}\,{\rm GeV}\lesssim m_{\rm dm}\lesssim 10^{12}\,{\rm GeV}$. As such, this observations (upper limits) can be used to constrain the DM lifetime in this mass range.

3.2 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 8 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 1. Fig. 9 shows a combination of extragalactic and Galactic limits together and include a comparison with previous results in the literature [25, 27, 28].

The left panel of Fig. 8 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 $\tau_{\rm dm}$ due to Galactic $\gamma$ rays. The PAO data bounds the DM lifetime $\tau_{\rm dm}\gtrsim 10^{30}\,$s for $10^{10}\,{\rm GeV}\lesssim m_{\rm dm}\lesssim 10^{15}\,{\rm GeV}$. In the mass range $10^{8}\,{\rm GeV}\lesssim m_{\rm dm}\lesssim 10^{10}\,{\rm GeV}$, KASCADE-Grande gives the most stringent constraints, i.e., $\tau_{\rm dm}\gtrsim 10^{29}\,$s. We note that these results are consistent with Ref. [28]. Finally, Fermi-LAT constrains the DM lifetime at roughly $\tau_{\rm dm}\gtrsim 10^{28}\,$s for $10^{3}\,{\rm GeV}\lesssim m_{\rm dm}\lesssim 10^{6}\,{\rm GeV}$. However, constrains obtained using $p+\bar{p}$ and $\bar{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. 8. 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}\,{\rm GeV}\lesssim m_{\rm dm}\lesssim 10^{8}\,{\rm GeV}$ where we find that the constrains on the neutrino flux using IceCube observations are the most stringent, i.e., $\tau_{\rm dm}\gtrsim 10^{28}\,$s. This is consistent with limits reported in Ref. [25]. 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. [97] in $m_{\rm dm}\leq 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 CRPropa instead of analytic methods described in that reference. Reference [27] reports a qualitatively similar result for $10^{7}\,{\rm GeV}\lesssim m_{\rm dm}\lesssim 10^{12}\,{\rm GeV}$, except that their bound is a factor of a few weaker. In addition, we find that for $10^{13}\,{\rm GeV}\lesssim m_{\rm dm}\lesssim 10^{16}\,{\rm 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 [98, 99, 100, 101, 102, 103].

4 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 104 and 1016 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. 9. Dark matter less massive than 108 GeV is most stringently constrained by unresolved diffuse gamma-ray emission measured by Fermi-LAT. For dark matter with much heavier masses above $\sim$1010 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 1030 s. For masses between 108 and 1010 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 1028 s throughout the wide mass range between 104 and 1016 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 $b\bar{b}$ 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.

Acknowledgments

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 CRPropa, 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. 2, 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.