Phenomenology of superheavy decaying dark matter from string theory

We study the phenomenology of superheavy decaying dark matter with mass around $10^{10}$ GeV which can arise in the low-energy limit of string compactifications. Generic features of string theory setups (such as high scale supersymmetry breaking and epochs of early matter domination driven by string moduli) can accommodate superheavy dark matter with the correct relic abundance. In addition, stringy instantons induce tiny $R$-parity violating couplings which make dark matter unstable with a lifetime well above the age of the Universe. Adopting a model-independent approach, we compute the flux and spectrum of high-energy gamma rays and neutrinos from three-body decays of superheavy dark matter and constrain its mass-lifetime plane with current observations and future experiments. We show that these bounds have only a mild dependence on the exact nature of neutralino dark matter and its decay channels. Applying these constraints to an explicit string model sets an upper bound of ${\cal O}(0.1)$ on the string coupling, ensuring that the effective field theory is in the perturbative regime.


Introduction
Despite various lines of evidence for the existence of dark matter (DM) in the Universe [1], its nature remains an important unsolved problem at the interface of cosmology and particle physics.Weakly interacting massive particles (WIMPs) have been a promising candidate and have driven theoretical and experimental research in this area for a long time.The lack of any positive signals in direct, indirect, and collider searches has however led to the consideration of DM candidates much lighter or much heavier than the typical mass range for WIMP-like particles.The well-known bounds on the annihilation cross section imply that freeze-out in a standard thermal history results in DM overproduction for masses larger than 100 TeV [2], thereby requiring alternative scenarios for obtaining the correct relic abundance of superheavy DM with m DM ≫ 100 TeV.
Superheavy DM can easily evade direct detection experiments and collider searches.Standard indirect detection signals from DM annihilation are also highly suppressed in this case.This is because, in addition to small annihilation cross sections, such a signal scales like n 2 χ,0 (n χ,0 being the current number density of DM particles).The situation will improve if DM is not absolutely stable as signals from decaying DM scale like n χ,0 .This scenario is conceivable because symmetries that are typically invoked to make DM stable may be broken by high-energy effects (from, for example, string theory or quantum gravity).While this could provide an observational window to probe superheavy DM, a viable scenario must include lifetimes that are orders of magnitude longer than the age of the Universe [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18].This leads to another question on how the very tiny couplings required for such long lifetimes arise.
In this paper, we address both questions regarding the origin of the abundance and indirect signals of superheavy DM in the context of an explicit construction within string theory.Successful inflation in this model results in a high scale of supersymmetry (SUSY) breaking (in the absence of fine tuning).The lightest SUSY particle (LSP) is stable by virtue of R-parity conservation, and is hence the DM candidate.The model gives rise to a non-standard cosmological history that involves one or more epochs of early matter domination (EMD) driven by string moduli.This can yield the observed relic abundance for a DM mass in the 10 10 − 10 11 GeV range [19].As we will see, this model can also accommodate decaying DM through tiny R-parity violating (RPV) interactions whose smallness is due to exponentially suppressed non-perturbative terms that are induced by instantons.
We find that, for setups which realise the Minimal SUSY Standard Model (MSSM), the resulting mass spectrum implies a Bino-like LSP.Small RPV effects lead to three-body decays of DM to Standard Model (SM) particles.We calculate the spectrum of high-energy gamma rays and neutrinos from these decays and constrain the DM mass-lifetime plane by imposing constraints from current experiments as well as sensitivity limits of future observations.These limits have a rather mild dependence on the exact RPV channel and do not change much even for a Higgsino-like or Wino-like LSP.The tightest bounds on the DM lifetime happen to be in the 10 9 − 10 11 GeV mass range.As we will see, the resulting limits on the size of RPV couplings in the model imply an upper bound of O(0.1) on the string coupling, which ensures the effective field theory remains in the perturbative regime.
This paper is organized as follows.In Sec. 2 we will describe a string construction that can naturally realise a SUSY model with superheavy decaying DM.We then discuss how various RPV terms make the LSP unstable in Sec.2.3.In Sec. 3 we discuss the relevant current and future experiments that can constrain this type of DM by means of highenergy gamma rays and neutrinos.We describe the methodology used in the numerical computations in Sec. 4. In Sec. 5 we show for the first time the model-independent neutrino and gamma-ray constraints on superheavy DM decaying into three-body final states, as well as the sensitivity of current and future experiments to the string model presented here.We conclude the paper and give further outlook in Sec. 6. Technical details on the renormalization group equations (RGEs) are given in App.A, while the case of superheavy sneutrino DM with two-body decays is described in App.B.
2 Superheavy dark matter from string theory

General features of 4D string models
Low-energy 4D string models have been explored in detail in recent decades, focusing in particular on the type IIB framework where moduli stabilisation is better under control.
These studies revealed the following generic features: • High SUSY breaking scale: Statistical studies of the distribution of the SUSY breaking scale in the string flux landscape have shown a (power-law or logarithmic) preference for higher scales of SUSY breaking [20,21].1 • High inflationary scale: Realising inflation and low-energy SUSY in string models is a difficult task [23].This is indeed the case in most of the known string inflationary models [24].The main obstacle comes from the requirement of matching the observed amplitude of the primordial density perturbations which fixes the inflationary energy scale at relatively high energies which, in turn, fixes also the value of the gravitino mass during inflation.In principle, the gravitino mass could evolve from the end of inflation to today [25,26], or SUSY particles could be much lighter than the gravitino due to sequestering in the extra dimensions [27,28], but these two features do not seem to be generic.On the other hand, high-scale inflation seems to imply also a high SUSY spectrum in general, regardless of the statistics in the string landscape.
• Moduli domination: 4D string compactifications come along with a plethora of moduli which are new scalar fields with gravitational couplings to matter.During inflation, these fields are naturally displaced from their minimum due to the inflationary energy density.Due to this initial misalignment, which is generically of order the Planck scale [29], the moduli store energy when they start oscillating around their minima which happens when the Hubble scale H becomes comparable to their mass.Being redshifted more slowly than radiation (since they behave as matter), the moduli quickly come to dominate the energy density of the Universe, giving rise to EMD epochs which can alter the standard thermal history from the end of inflation to Big-Bang Nucleosynthesis.In particular, the decay of the moduli can dilute any DM production which took place prior to their decay [19,[30][31][32][33][34][35].
• MSSM with RPV: MSSM-like models can be realised with stacks of intersecting magnetised D7-branes wrapped around internal 4-cycles whose size determines the values of the SM gauge couplings [36].In these constructions U (1) B−L is a gauge symmetry at the fundamental level but it arises as an effective global symmetry below the Stückelberg mass of the corresponding gauge boson which is around the string scale.This effectively global U (1) B−L symmetry is exact at perturbative level but it gets broken by non-perturbative effects like stringy instantons.Some of these instantons break the continuous U (1) B−L symmetry down to a discrete Z 2 symmetry which can be identified with R-parity [37], while other instanton contributions (with a different zero mode structure) would induce tiny R-parity violating couplings proportional to e −S inst ≪ 1, where S inst is the instanton action [38].

A string model for superheavy dark matter
The generic features of 4D string compactifications described in Sec.2.1 might point towards an interesting scenario of superheavy decaying DM.In fact, Ref. [19] has shown that high scale SUSY breaking and dilution from moduli domination can give rise to the right relic abundance of a superheavy WIMP χ with m χ ∼ 10 10 GeV.The model is built within type IIB Large Volume Scenarios (LVS) where the complex structure moduli and the dilaton are stabilised by 3-form fluxes, while the Kähler moduli are fixed by combinations of perturbative (in both α ′ and g s ) and non-perturbative effects [39,40].The minimal model involves 3 Kähler moduli (focusing on the canonically normalised modes): the overall volume modulus ϕ and two blow-up modes, σ and ρ.One of them, which we will identify with σ, drives inflation as in Kähler moduli inflation [41,42] and it is wrapped by a hidden sector stack of D7-branes, while the other, ρ, is a spectator field during inflation and it is wrapped by the MSSM stack of D7-branes.During inflation, the volume modulus ϕ is shifted from its late-time minimum, which we set at ϕ = 0, by an amount of order Y ϕ = ϕ 0 /M p (where M p ≃ 10 18 GeV is the Planck scale), and so drives an EMD epoch after the end of inflation [29].
The moduli break SUSY and generate a non-zero gravitino mass m 3/2 ≃ M p /V ≪ M p , where V is the extra-dimensional volume in units of the string length.Their gravitational interaction with the MSSM induces soft SUSY breaking mass terms for squarks, sleptons and gauginos.At the Kaluza-Klein scale M KK ≃ M p /V 2/3 , all SUSY scalars have a common mass given by m 0 , while all gaugino masses are given by M 1/2 .The mass spectrum at M KK looks like: where to work out the exact coefficient of proportionality between the soft terms and m 3/2 we actually considered a slightly more involved situation where the MSSM is built with two stacks of intersecting and magnetised D7-branes wrapped around two intersecting blow-up cycles; one blow-up modulus is stabilised by setting a Fayet-Iliopoulos term to zero, while the other is fixed by an ED3-instanton.The requirement of reproducing the right amount of density perturbations during inflation fixes the gravitino mass around 10 10 GeV, corresponding to volumes of order V ∼ 10 6 .This corresponds also to the typical mass scale of all SUSY particles and the two moduli σ and ρ.On the other hand, the volume modulus is much lighter since m ϕ ∼ 10 7 GeV.
At the end of inflation, an initial reheating epoch is due to the decay of the inflaton σ which decays mainly into hidden sector degrees of freedom since its coupling to MSSM fields is geometrically suppressed by the distance in the extra-dimensions between the inflaton 4-cycle and the 4-cycle supporting the MSSM D7-stack.Hence the production rate of MSSM particles, and in particular of the LSP, from the inflaton decay is very small since [19]: In fact, the DM number density n χ at inflaton decay reads (with Γ σ = Γ σ→hid +Γ σ→MSSM ≃ Γ σ→hid ): where the inflaton number density is n σ ≃ 3Γ 2 σ→hid M 2 p /m σ , while Br odd is the branching ratio for the inflaton decay into R-parity odd particles which then decay into DM.
The decay of ρ can be safely ignored since this modulus would decay when it is not dominating the energy density.On the other hand, the volume mode ϕ has time to come to dominate the energy density (when H = H dom ≃ Y 4 ϕ Γ σ ) before decaying.Its decay cannot reproduce any DM particles since ϕ is lighter than any superpartner, but it dilutes the χ-particles produced at inflaton decay.To get the final DM relic abundance we need therefore to redshift (2.3) during the radiation-dominated (RD) epoch from σ-decay to the onset of ϕ-domination, which amounts to a multiplicative factor (H dom /Γ σ ) 3/2 , and during the matter-dominated epoch from ϕ-domination to ϕ-decay, which yields another multiplicative factor of order (Γ ϕ /H dom ) 2 .Putting all these factors together, we end up with [19]: where we have normalised by the entropy density s = 2π 2 g * T 3 rh /45 at the reheating temperature from ϕ-decay T rh ∼ Γ ϕ M p ∼ m ϕ m ϕ /M p .The relic density (2.4) has to match the observed value [43]: Reference [19] has performed a detailed analysis, showing that (2.4) and (2.5) indeed match for m χ ∼ 10 10 GeV.This can be easily seen by noticing that, for such a large mass scale, (n χ /s) obs ∼ 10 −20 , which is reproduced by (2.4) for Y ϕ ∼ 0.1 (as computed in [29]), Γ σ→MSSM /Γ σ ≃ 10 −12 (from (2.2) for V ∼ 10 6 ) and T rh /m σ ∼ (m ϕ /m σ ) m ϕ /M p ∼ 10 −7 (from (2.1) and V ∼ 10 6 ).Finally, we perform the RGE evolution of the gaugino masses and SM gauge couplings at one-loop order to estimate the mass spectrum of SUSY particles at lower energies in the model described above.We begin the evolution at the Kaluza-Klein scale with the following boundary conditions: , with the gravitino mass set to m 3/2 = 10 10 GeV and the gauge couplings given by their approximate unification value g a = 0.7.From there we run down to the mass scale of the volume modulus m ϕ ∼ 10 7 GeV to obtain the gaugino and scalar masses at low energies (see App.A for details).In Fig. 1 we depict the resultant RGE running with the scalars labeled in the legend and the Bino, Wino, and gluino shown by M 1 , M 2 , and M 3 respectively.The squarks receive strong contributions from the gluino mass, pushing their masses higher than their initial value, while the slepton masses run more mildly, remaining closer to their initial value.We see that in this example the gluino is the heaviest sparticle, while the Bino is the LSP.

Decaying dark matter and three-body decay signatures
As explained in Sec.2.1, an important generic feature of MSSM model-building in string theory is the presence of exponentially small RPV couplings generated by stringy instantons [37,38].These couplings can therefore make the LSP unstable, leading to a decaying DM scenario if the lifetime of the χ-particles is not smaller than the age of the Universe.
To take into account RPV, we shall include the following trilinear contributions to the superpotential [44]: where the first two terms violate lepton number, while the third one violates baryon number.
Refs. [45,46] for details.In particular, these vertices make a sfermion decay into two Feynman diagrams illustrating the decay of the LSP χ 0 into three SM fermions, mediated by a sfermion f .The RPV vertex is highlighted by the white dots.
fermions, allowing for instance the sneutrino or the neutralino, χ 0 , to decay.Indeed in the MSSM, the neutralino-sfermion-fermion vertex is given by: where k specifies the sfermion mass eigenstate, while P L and P R stand for left and right projectors.The couplings are proportional to the sfermion mixing matrix and to all the neutralino eigenstates (Bino, Wino and Higgsino).The sfermion in Eq. (2.10) is then allowed to decay via one of the RPV terms, in Eq. (2.7), and gives as a result a three-body final state, as illustrated in Fig. 2. The LLE term produces a decay into three leptons, one of which should be a neutrino, while the other two terms produce predominantly quarks.
In what follows, we will consider the SUSY neutralino, χ 0 , as the DM candidate, while in App.B we will discuss the case of a sneutrino DM candidate.The decay of the neutralino will take place regardless of its composition.Note, however, that for the Higgsino case, since the coupling is proportional to the mass of the fermion, there will be a strong suppression due to the Yukawa coupling.From now on we will work in the pure Bino approximation (see Fig. 1) as discussed below.The χ 0 decay due to the RPV terms is given by [45][46][47]: where α is the fine-structure constant at the scale of the DM mass.All decay widths have been computed in the limit of degenerate heavy sfermion masses.The LQD and U DD vertices that depend on coloured particles scale as N c and N c ! respectively, with N c = 3 being the colour factor.The neutralino lifetime τ χ 0 = Γ −1 χ 0 should be larger than τ χ 0 ≳ 10 26 s to be compatible with current bounds on decaying DM, see [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18], and to be larger than the age of the Universe.For m χ 0 ≃ m f ≃ 10 10 GeV, this would require very small RPV couplings of order λ ijk ≲ 10 −26 , in the case of the LLE RPV term.Interestingly, these values are in the right ballpark expected in 4D string models.In fact, RPV couplings are generated by ED3-instantons wrapped around the 4-cycle τ SM supporting the MSSM D7-stack, and so λ ijk ∼ e −S inst where the instanton action is given by S inst = 2πτ SM .In LVS models τ SM is a blow-up modulus which is typically stabilised at values of order the inverse of the string coupling g s , τ SM ≃ g −1 s .Therefore we expect λ ijk ∼ e −2π/gs which, for values of the string coupling g s ≲ 0.1 that keep the effective field theory in the perturbative regime (which can be guaranteed by an appropriate choice of 3-form flux quanta), would give indeed λ ijk ∼ e −20π ≲ 5 × 10 −28 .

Astrophysical signatures of superheavy decaying dark matter
Superheavy DM particles χ distributed in galactic and extra-galactic structures can decay into ordinary matter by producing Ultra-High-Energy (UHE) particles.Here we focus on UHE neutrino and gamma-ray fluxes produced by DM decays.UHE charged cosmic rays are also an interesting probe for DM decay.However, their understanding is subject to larger astrophysical uncertainties, because of their propagation in the galaxy, with respect to neutrinos and gamma rays which trace the source.Our analysis on neutrino and gammaray data is based on Refs.[3,4] which we refer the reader to for further details.

Theoretical astrophysical fluxes
The neutrino and gamma-ray fluxes produced by DM decay originate in principle from both galactic and extra-galactic sources.However, the gamma-ray flux coming from extragalactic production is expected to be substantially attenuated, due to the absorption by γγ scattering happening at high energies and high redshift, and can be subsequently neglected.Additionally, the secondary photon emission from Inverse Compton scatterings of electrons and positrons contributes to gamma-ray observations at energies E γ ≲ 10 5 TeV, and therefore can be safely neglected in our analysis based on UHE gamma rays.Hence, the dominant component at E γ ≳ 10 5 TeV is the prompt galactic emission which gives the following expected flux: where τ γγ represents the total optical depth of photons to γγ collisions and ρ χ = ρ s [r/r r (1+ r/r s ) 2 ] −1 with r s = 24 kpc and ρ s = 0.18 GeV cm −3 is the Navarro-Frenk-White galactic DM density profile [48,49].The optical depth takes into account the photon pair production from the interaction of the gamma rays with CMB radiation, starlight and infrared light (see Refs. [4,5] for details on its computation).While the CMB photon density is homogeneous and thermally defined by the CMB temperature (≃ 2.348 × 10 4 eV), the photon density of starlight and infrared light strongly depends on the position in the Milky Way.For the latter, we take the photon density from the GALPROP code [50].The density profile is expressed as a function of the distance r from the center of our galaxy.The integration is performed along the line of sight s from the Earth, so that r can be written as a function of the colongitude ℓ and colatitude b, according to the following expression: where we identify r ⊙ = 8.5 kpc as the distance between the Earth and the galactic center.Current UHE gamma-ray observatories place bounds on the diffuse gamma-ray flux, which can be linked to the integral gamma-ray flux, averaged over the field of view of each experiment Ω exp (see supplementary Tab.I in Ref. [6]), which takes the following expression: Contrarily to gamma rays, neutrinos travel unimpeded and the neutrino flux from DM decays has both a galactic and extra-galactic component.The first one is given by: where α is the flavour of the (anti)neutrino considered.The important feature of extragalactic neutrinos is that their absorption in the inter-galactic medium is negligible.Hence they produce a sizeable flux that can be written as: where H is the Hubble parameter, ρ c is the critical density and Ω χ is the DM relic density, according to the latest Planck data [51].This flux depends on the cosmological redshift z.
Both galactic and extra-galactic neutrinos travel a long distance before reaching Earth and they oscillate very rapidly with respect to the propagation scale, so that the detected neutrinos have a composition of flavour eigenstates given by averaged fast oscillations.This feature is taken into account by defining a total flux given by the sum of each flavour flux.Moreover, experiments are not sensitive to the neutrino arrival direction, but only to its average.We keep that into account by performing an integration over the solid angle.As a result, the total combined flux is: ) depend on the energy spectrum, dN ν,γ /dE ν,γ , that denotes the energy distribution of neutrinos and photons deriving from the three-body or two-body DM decays.This is a central key feature of this paper that needs to be carefully evaluated to have proper estimates of the sensitivity of current and future experiments.We will discuss how we properly compute this quantity in Sec. 4, after having introduced the experimental neutrino and gamma-ray probes, as well as the corresponding statistical approach we adopt to compute the limits on the DM lifetime τ χ for different DM masses m χ .

Future neutrino telescopes
When neutrinos reach the Earth, they can be detected by the hadronic or electromagnetic showers that they create when interacting in a medium.Many new experiments will be available in the near future to collect these signals for UHE neutrinos.In general, the showers originating from UHE neutrinos are characterized by emissions in the radio spectrum, a phenomenon known as the Askaryan effect.However, since the neutrino fluxes at these high energies are low, the experiments should be equipped with a large interaction volume, so that the main mediums that have been considered so far are air, water and ice.In the following we give a brief description of the main experiments at play in the next decades: RNO-G [52] the Radio Neutrino Observatory in Greenland is a neutrino telescope exploiting an ice medium for neutrino interaction.UHE neutrinos can produce hadronic or electromagnetic showers inside the ice which can then be collected by the antennas equipped to each one of the 35 stations the experiment is made of.
IceCube-Gen2 [53,54] the next generation of the IceCube experiment is located in Antarctica, exploiting the natural ice for neutrino interaction.The detection of the showers relies on the radio array made of 200 stations.
GRAND [55] the Giant Radio Array for Neutrino Detection is a planned observatory sensitive to cosmic rays, gamma rays and neutrinos.Its main strategy for neutrino detection focuses on Earth-skimming underground τ neutrinos: they are astrophysical ν τ crossing the Earth medium at a very small angle with the surface, or passing through mountains.In this environment, neutrinos can interact with the medium by charged-current interaction and can produce τ leptons which can eventually decay outside of the rock, producing an electromagnetic shower which can be detected by the large number of antennas foreseen in the design of GRAND.Here, we consider two different GRAND configurations with 10k and 200k radio antennas, respectively.
Superheavy decaying DM offers a natural source of neutrinos that can propagate through the galaxy and reach the Earth.In particular, given a DM candidate of mass m χ and lifetime τ χ , we can compute the expected number of neutrino events in a given telescope, assuming an observation period of T obs = 3 yr: where E max = m χ /2 is the maximum neutrino energy achieved in the DM decay.The neutrino flux from DM decay has been previously defined in Eq. (3.6), while we define the effective area of the experiment with A eff .This can be defined through the flux sensitivity S(E ν ), which is the quantity provided by the experimental collaborations, and the observation time: where 2.44 is the number of neutrino events for each energy decade.Following Ref. [3], we will consider two main contributions (taken from Ref. [55]) for the UHE neutrino flux coming from astrophysical sources: • Cosmogenic neutrinos: They originate from collisions between high-energy cosmic rays and the CMB.Their contribution to the flux of UHE neutrinos is affected by important theoretical uncertainties and, to be conservative, we will consider the most optimistic estimate on this contribution in order to obtain the weakest limits on DM; • Pulsar neutrinos: The extreme astrophysical environment of newborn pulsars might give rise to one of the highest UHE neutrino fluxes.
Such astrophysical models predict the observation of a certain number of astrophysical neutrino events, N astro ν , which can be computed from the neutrino flux through an equation analogous to Eq. (3.7).Hence, we forecast conservative limits on the DM lifetime by assuming the observation of N obs ν neutrino events (with N obs being a Poissonian stochastic variable with mean N astro ν ) and by taking the test statistic: where the likelihood L is assumed to be a Poisson distribution.According to this test statistic [56], we can place future bounds on the DM lifetime τ χ at 95% CL by taking TS(m χ , τ χ ) = 2.71.

Current gamma-ray experiments
UHE photons and cosmic rays reaching the Earth's atmosphere will interact producing extensive air showers (EAS), generating cascades of charged particles that can leave traces in the detectors placed on the Earth's surface.Gamma-ray experiments usually employ one or more scintillation arrays or huge Čherenkov stations.The latter aim at detecting the Čherenkov light coming from the particles in the EAS travelling through the atmosphere, that is collected by photomultiplier arrays.One of the main obstacles for the measurement of gamma rays with ground-based telescopes is the high level of residual charged cosmic ray background.EAS can not be distinguished by background charged cosmic rays based on directional information.
In the present study, we focus on the following current and decommissioned experiments which have placed several bounds on the UHE gamma-ray diffuse flux above 10 5 GeV (see Fig. 1 in Ref. [4]).[57] the Chicago Air Shower Array-MIchigan Muon Array experiment was a surface array of scintillation counters located in Utah (USA), able to infer the properties of EAS and it was active until 1998.The main detection strategy relied on discriminating gamma-ray primary showers from hadronic (cosmic ray)-generated ones, by identifying muon-less or muon-poor events.This principle is based on the fact that muons originate mainly from charged pion and kaon decays which are abundant in a hadronic shower, while the muon content in an electromagnetic shower is mainly coming from photoproduction events whose interaction rate is suppressed with respect to nucleus-nucleus events.

KASCADE and KASCADE-GRANDE [58] the KArlsruhe Shower Core and Array
DEtector was located at the Karlsruhe Institute of Technology (Germany) and operated until the end of 2012.It was made of three detector systems: a large field array, a muon tracking detector and a central detector.The experiment's aim was to measure the electromagnetic, muonic and hadronic components of EAS.
TA [59] the Telescope Array experiment operates in Utah (USA) since 2008 and it is made of a surface detector with plastic scintillators and three fluorescence detectors.It obtained limits on diffuse photon flux with energies greater than EeV.
PAO [60] the Pierre Auger Observatory is a surface detector array of Čherenkov stations, sensitive to the electromagnetic, muonic and hadronic components.An additional fluorence detector completes the setup to observe the longitudinal development of the shower.It operates in Argentina and it performed a search for UHE photons with energies above EeV.
Given the upper bounds on the angle-averaged integral gamma-ray flux from each experiment, we compute the limits on DM by simply performing a χ 2 analysis using the test statistic: where the index i runs over the experimental energy bins, Φ data γ = 0 (no observation of photons) and σ i is the corresponding standard deviation at 68% CL of each data point.The limits on the DM lifetime at 95% CL can be obtained by taking TS = 2.76, a threshold value which derives from the constraint that the gamma-ray measurements cannot be negative [4].

Methodology and energy spectra
The DM decay rates and final states have been computed using the UFO model that is publicly available in the FeynRules [61,62] database. 3It encodes the MSSM with trilinear vertices for RPV interactions and it is based on [45,63].It allows for the choice of the SUSY mass spectrum, in such a way that either the neutralino or sneutrino can be set as the LSP.
The SUSY mass spectrum has been chosen in such a way that the gauginos and sfermions are all at the same mass scale, which is given by the gravitino mass, as explained in Sec. 2. In order to optimise the generation, we have discarded irrelevant and time-consuming contributions in the decay matrix element computation.In particular, we have set to a much heavier mass scale all SUSY particles that do not participate in the decay.
The UFO model files have been loaded in MadGraph5 aMC@NLO [64] to obtain the decay rates, including the three-body decays of the neutralino.The kinematics of the final states is fixed for a two-body decay, while it is not for a three-body final state.In this case the distribution in energy of the final state particles of the hard process needs to be simulated accurately.Tab. 1 shows all possible final states for the neutralino decay.We have turned on one RPV vertex at a time, in order to study separately the three channels LLE, LQD and U DD.
Considering the case of neutralino DM, the LLE vertex gives rise to three leptons of different flavour in the final state, one of them always being a neutrino, such as χ 0 → µ + e − ν τ , χ 0 → e + e − ν µ and so on.There are 36 different final states.We run Mad-Graph5 aMC@NLO to simulate the hard process, asking for 10 6 events per decay, which ensures a statistic large enough to sample accurately each single final state.We proceed in a similar manner for the LQD and U DD vertices, requiring 10 6 events per decay.The former gives rise to final states with one lepton and two quarks, such as χ 0 → e − u d, χ 0 → ν e s d, for a total of 108 final states.The third RPV vertex gives rise to 18 final states with different combinations of quarks, e.g.χ 0 → u s b, χ 0 → b t s, ....The hard process final states of the three-body neutralino decay are highly energetic, because we are considering DM particles with masses m χ ≃ 10 10 GeV.Every final state will generate a shower, producing a cascade of new particles.Given the high energy of the process, it is important to properly include electroweak corrections that affect the final energy spectra of neutrinos and gamma rays (for discussion, see [65,66]).Notice that the electroweak evolution is collinear, and therefore local.Hence, given the mother particle and its energy, it is enough to perform the evolution.This can be done consistently using the electroweak corrections obtained with the tool HDMSpectra [66].
HDMSpectra provides electroweak-corrected spectra for γ, ν's and protons originating from the two-body annihilation (or decay) of a generic DM particle, for energies up to 10 19 GeV: for which the hard spectrum is highly peaked.However, in our case we need to compute the spectra for a three-body decay process, distributed accordingly to a continuum of energies, hence a convolution operation between the hard spectra and HDMSpectra results is necessary.We proceed as follows.For each RPV vertex, we sum over all the neutralino three-body decay processes, getting P different particles overall, e.g.: The convolution operation requires us to consider each p i , compute the electroweak shower spectra with HDMSpectra and sum over every particle.To avoid doing that for each event, we bin the particles in the hard spectra and do a convolution.This is possible assuming that the particles belonging to each bin in the hard spectra have the same energy, and in particular equal to the mid value of the bin.This is true if the binning scheme is fine enough to have very small bins, which is realised by taking 1000 bins in our case.The convolution is made by considering that each binned particle is coming from a fake two-body annihilation with a proper center-of-mass energy to match the energy of the bin.The result will be the electroweak-corrected spectra coming from the three-body decay of the neutralino.More precisely, the hard process spectrum generated by MadGraph5 aMC@NLO is stored in an event file in the LHE format [67], containing the different generated processes and kinematics for each event.This file can be parsed to extract the hard spectra of the final states.Each hard spectrum of the particle p i , with i ∈ [1, P ], is a one-dimensional histogram with M bins in the variable x = E/m χ 0 , where E is the energy of the particle p i .Each bin j, with j ∈ [1, M ], of the hard spectrum of the particle p i will contain n ij particles and we have: where D is the total number of generated events for the decays in (4.2).Furthermore, the histogram is normalised according to D, so that: w ij is the probability that the three-body decay of χ will yield the particle p i with energy fraction x j (belonging to the j-th bin).
Then we need to compute the shower spectrum for each of the three-body decays.Every particle p i coming from the decay will produce a shower, generating several different particles: where H is the number of different particles types generated (γ, ν's, . . .), each one being a state s k , with k ∈ [1, H].We remark that the process written in (4.5) happens for each particle in the histogram of p i , where every bin j is associated to the weight w ij and contains particles (p i ) j with energy E j = x j m χ 0 .For each bin j, we run HDMSpectra simulating the following process: where χ j is a fake DM with a mass consistent with the energy of the particles (p i ) j : Because HDMSpectra samples its output on the basis of m χ , with an array binned in y j = E/m χ j , we substitute: to have results in terms x.HDMSpectra returns the values of (dN k /dy j ) i , which represents the electroweak-corrected spectra of the final state s k , as coming from the shower of the hard spectrum particle p i with energy E j .This spectrum should be rescaled appropriately: for every bin j of the particle p i hard spectrum and then for any particle p i of the decay processes.The total spectrum of the particle s k is eventually given by summing all of the contributions for each bin j and for each particle p i , weighting each result on the basis of the probability w ij discussed before: In Fig. 3 we show a comparison between the shower spectra in photons and neutrinos, coming from the two-and three-body decay spectra for the LLE vertex in the case of a sneutrino and a neutralino DM candidate, respectively, fixing m DM = 10 10 GeV.As expected, the three-body decay spectra are shallower and peaked towards lower x with respect to the case of the two-body decay spectra.

Neutrino and gamma-ray bounds for three-body decays
After describing the statistical procedure in Sec. 3 and the computation of the energy spectra in Sec. 4, we can now obtain the lower bounds on the DM lifetime τ χ for DM masses m χ in the range from 10 7 to 10 15 GeV.Here, we focus on the three-body decay channels reported in Tab. 1, which occur in the case of neutralino DM.The case of two-body decay channels (sneutrino DM) is discussed in App.B. We consider that the decay of DM occurs through only one of the couplings LLE, LQD and U DD (see Tab. 1) at a time.Moreover, we assume that such couplings are independent from the particle flavours.This results in 36, 108, 18 different final states (taking into account different flavour combinations) with equal branching ratio for the LLE, LQD, U DD scenarios, respectively.Hence, we first discuss the model-independent constraints for each coupling, and then report the constraints projected in the parameter space of the string model featuring superheavy neutralino DM.
In Figs. 4 and 5 we report the bounds obtained with current gamma-ray data and with forecast neutrino data, respectively, in the case of the three different decay channels.In the former figure, different lines correspond to the bounds placed by different telescopes.In the latter, the left (right) panels show the bounds for the future neutrino telescopes RNO-G and IceCube-Gen2 (GRAND10k and GRAND200k).Moreover, the neutrino forecast bounds are shown as bands which correspond to the uncertainty on the different scenarios for the expected neutrino observations after 3 years of data-taking according to the astrophysical pulsar and cosmogenic neutrinos, as well as to zero neutrino events (null observation of neutrino events).The latter provides the strongest constraints on the DM properties from neutrino data.The common behaviour of all the constraints is due to the fact that each experiment is sensitive to a specific energy range.This translates into a specific interval of the DM mass that can be probed.In Fig. 6 we make a direct comparison between the constraints placed by all the gammaray experiments (red lines) and the ones which will be placed by IceCube-Gen2 assuming zero neutrino events (blue lines).Different line styles correspond to different DM decay channels.As can been seen, we generally expect future neutrino data to provide stronger bounds on the DM lifetime for DM masses below 10 11 GeV.This behaviour depends on the decay channels.In particular, for the pure leptophilic LLE channel, which involves only leptons (solid lines), we expect an improvement on the constraints from the neutrino data for a larger range of DM masses.On the other hand, the impact of neutrino data is expected to be very limited in the case of the pure hadronic U DD channel (dot-dashed lines).This is simply due to the different energy injection into UHE photons and neutrinos for the   different decay channels.We emphasise that our forecast analysis is not strongly sensitive to spectral features and to angular DM anisotropy since it is based on energy-integrated and angle-averaged quantities.However, the neutrino constraints could be further improved with a more refined analysis once actual data will be taken by the future telescopes.
All of the previous bounds hold for a generic DM candidate with mass m χ and threebody decay channels summarised in Tab. 1.For the specific model discussed in Sec. 2, we can recast the constraints on the DM lifetime τ χ into model-dependent bounds on the RPV coupling λ by means of Eqs.(2.11), (2.12) and (2.13) and, remarkably, on the string coupling g s through the relation λ ∼ e −2π/gs .To obtain such bounds, indeed one needs to also define the mass of the mediator (see Fig. 2) as well as the RGE running of the finestructure constant which has to be computed at an energy scale Q ∼ 10 10 GeV.As shown in Fig. 1, at this scale the mass of the neutralino DM is m χ 0 = 7.3 × 10 9 GeV, while the averaged masses m f of the mediators are 8.0 × 10 9 GeV, 9.1 × 10 9 GeV and 9.9 × 10 9 GeV, for the three different couplings LLE, LQD and U DD, respectively.Moreover, the finestructure constant is α(m 2 χ 0 ) = 0.01172.In Fig. 7 we report the constraints on the three different couplings as a function of the neutralino DM mass m χ 0 .The red dashed lines correspond to the gamma-ray constraints, while the blue solid ones to the strongest bounds that can be placed with future neutrino observations assuming zero neutrino events (no astrophysical backgrounds).The vertical thin black lines mark the prediction of the model for the DM mass, while the horizonal lines show a benchmark value for the string coupling g s .Hence, these bounds can be then translated into constraints on g s as done in Fig. 8 where we fix the DM mass to be m χ 0 = 7.3 × 10 9 GeV and vary the mediator masses.Here, the vertical lines correspond to the mediator masses obtained from the RGE running.Interestingly, for the DM mass predicted by the model, the allowed values for g s are g s ≲ 0.1.This implies that the phenomenological constraints from UHE gamma rays and neutrino observations coincide with the theoretical requirement on the string coupling to be in the regime where the low-energy effective field theory is under perturbative control.Let us also mention that in string compactifications the value of g s is determined dynamically as the vacuum expectation value of the dilaton field which is stabilised by background 3-  form fluxes.Given that different possible choices of these flux quanta are compatible with tadpole cancellation, a landscape of string vacua arises with an approximately uniform distribution of the string coupling, i.e.N vac (g s ) ≃ g s [21,22,68,69].One might therefore naively think that the phenomenological constraint λ ≲ 10 −30 from Fig. 7 could leave only a very sparse region of the landscape.This is however not true since the RPV coupling λ depends exponentially on g s , λ ∼ e −2π/gs , and so the distribution of λ turns out to be milder than logarithmic: Thus the increase of the number of vacua for larger values of λ is milder than logarithmic, implying that the region with λ ≲ 10 −30 contains still several string vacua which are in agreement with current data and can be probed with future DM observations.

Conclusions
In this paper we have studied indirect detection signals from decaying superheavy DM.Generic features of string constructions (namely a high scale of SUSY breaking and epochs of EMD driven by string moduli) provide a suitable framework to accommodate superheavy LSP DM and obtain the correct DM relic abundance [19].Moreover, tiny RPV couplings induced by stringy instantons naturally lead to a DM instability with decay lifetimes longer than the age of the Universe.
As a case study, we considered an explicit string model with high scale SUSY breaking that gives rise to successful inflation and yields the correct DM relic abundance in the 10 10 − 10 11 GeV mass range.We derived the mass spectrum of SUSY particles for a benchmark point that realises the MSSM at low energies and found that the LSP DM is a Bino-like neutralino.The RPV couplings lead to three-body decays of DM to quarks and leptons.
By adopting a model-independent approach, we computed the resulting gamma-ray and neutrino fluxes as a function of the neutralino mass and lifetime.We obtained constraints in the mass-lifetime plane by imposing the limits from current experiments and using the sensitivities of future neutrino observations.The most stringent limits on the lifetime, τ χ ≳ 10 30 s happen for m χ ∼ 10 9 − 10 11 GeV.We found that these bounds have a rather mild dependence on the specific decay channel (leptonic/hadronic/mixed as well as flavour combination of final states).Furthermore, they do not change significantly for Higgsino-like or Wino-like neutralino DM.It is interesting to note that for the specific string model considered, these limits result in an upper bound of O(0.1) on the string coupling which is compatible with the effective field theory being in the controlled regime where perturbation theory does not break down.
Our work provides an explicit example for indirect detection of superheavy DM which would otherwise evade (in)direct detection and collider searches, in the context of a UVcomplete model from type IIb string compactifications.A natural extension of this model could include RH neutrinos (and their superpartners) that are usually introduced to address neutrino mass and mixing.This results in new final states for DM decay that include the RH neutrinos.This is worth exploring as such high energy RH neutrinos might provide an explanation for the recent anomalous events reported by IceCube and ANITA experiments (for example, see [70]).

B Sneutrino dark matter and two-body decay spectra
In the MSSM typically the sneutrino cannot be the LSP, as one can see from Fig. 1, unless the GUT unification condition of scalar and gaugino masses is uplifted.A possibility to have a sneutrino LSP, and hence a good DM candidate, is to extend the MSSM content with a right-handed neutrino superfield.The right-handed sneutrino can then be DM (with or without mixing with the left-handed sneutrino).In this way it is possible to tackle at the same time the neutrino mass issue as well (for example, see [73][74][75][76][77][78][79] for further discussions on how to achieve sneutrino DM).Sneutrino DM is made directly unstable by the RPV terms in (2.7).However, note that only the LLE and LQD terms can contribute to ν decay, because of the superfield L. The sneutrino decays into two fermions and the main decay diagrams are illustrated in Fig. 9.The expressions for the sneutrino decay into fermions are given by [80]: where the fermion masses have been neglected.The possible final states are given in Tab. 2. It is interesting to remark that it is not possible to have monochromatic neutrino lines from the two-body decay, because of the structure of the RPV terms.This is in contrast with the case of annihilating sneutrino DM (for example, see [81]).
From the numerical point of view, the sneutrino DM decay is simpler to handle than the one of χ 0 , to obtain the energy spectra.The sneutrino decays via LLE into two charged leptons of different flavour, and there are 12 possible contributions, such as ν → e + µ − , ν → τ + µ − , ....The LQD vertex gives rise to only 6 final states featuring one lefthanded quark and one right-handed down quark, because it contains only one L, such as ν → d d, χ 0 → s b, .... To simulate the sneutrino decay we used the tool HDMSpectra, which allows one to obtain the spectra for two-body DM decay.
In Fig. 10 we show the lifetime constraints on the basis of gamma-ray data and future neutrino predictions, for a sneutrino DM candidate, in the case of the two RPV decay channels.
The exclusion bounds from gamma rays and the sensitivity reach of neutrino telescopes are rather similar to the case of neutralino DM, as the information coming from the energy spectra is integrated over the full interval of sensitivity of the experiment.This has the effect of making less prominent the differences between two-body and three-body decay spectra.Figure 10: Current gamma-ray limits (top panels) and future neutrino limits (center and bottom panels) on the DM lifetime as a function of the DM mass in case of sneutrino DM ν which can decay into two particles as ν → ℓ + ℓ − (left panels) and ν → dd (right panels).A sum over different flavour combinations of the final states is performed.The neutrino bands cover the uncertainty on the astrophysical UHE neutrino background for DM searches.

Figure 1 :
Figure 1: One-loop RGE running of the gaugino and scalar masses as a function of energy scale Q, with boundary conditions as described in the text.M 1 , M 2 , and M 3 depict the Bino, Wino and gluino, respectively, while the scalars appear as shown.The vertical dashed line marks the modulus mass m ϕ ∼ 10 7 GeV.

. 6 )
Considering the angle-averaged flux means to neglect neutrino flux anisotropies which are particularly relevant for galactic neutrinos in the direction of the galactic center.This implies that the combined flux expressed by Eq. (3.6) allows us to understand only the flux normalization, staying conservative on the angular distribution.The neutrino and gamma-ray fluxes in Eqs.(3.1), (3.4) and (3.5

Figure 3 :
Figure3: Comparison between the shower spectra in γ (red) and ν (blue, averaged over the three flavours) from sneutrino (dashed lines) and neutralino (solid lines) DM candidates, both with m DM = 10 10 GeV and decaying via the LLE vertex according to two and threebody processes, respectively.

Figure 4 :
Figure 4: Current gamma-ray limits on the DM lifetime as a function of the DM mass for different experiments (different lines), assuming at a time one of the three different threebody decay channels arising from the RPV couplings (see Tab. 1): LLE (top left), LQD (top right), and U DD (bottom center).The horizontal dotted line acts as a reference.

Figure 5 :
Figure 5: Forecast neutrino limits on the DM lifetime as a function of the DM mass for different future neutrino telescopes, assuming at a time one of the three different threebody decay channels arising from the RPV couplings (see Tab. 1): LLE (top panels), LQD (center panels), and U DD (bottom panels).The bands cover the uncertainty on the astrophysical neutrino flux which leads to a different expected number of detected neutrino events after 3 years of observation.The horizontal dotted line acts as a reference.

10 7 Figure 6 :
Figure6: Comparison of the constraints on the DM lifetime placed with all current gamma-ray data (red lines) and with no neutrino events in IceCube-Gen2 after 3 years of data-taking (blue lines).The solid, dashed and dot-dashed lines refer to the LLE, LQD, U DD channels, respectively (see Tab. 1).

Figure 7 :
Figure7: Upper bounds on the LLE (top left), (top right) and U DD (bottom center) couplings from current gamma-ray (red dashed lines) and future neutrino (blue solid lines) data, as a function of the neutralino DM mass m χ 0 .For the three cases, the mediator mass is fixed according to the RGE running (see the main text).The thin black lines display a benchmark value for the string coupling g s = 0.09 and the model prediction for m χ 0 .

Figure 8 :
Figure 8: Upper bounds on the string coupling g s as a function of the mass ratio m f /m χ 0 with m χ 0 = 7.3 × 10 9 GeV.The different panels correspond to the different RPV couplings of the model: LLE (top left), LQD (top right) and U DD (bottom center).Red dashed lines (blue solid lines) refer to the current gamma-ray (future neutrino) limits.The vertical thin black lines display the model prediction for the mass ratio.

Figure 9 :
Figure 9: Feynman diagrams illustrating the decay of the LSP ν ℓ into two SM fermions.The RPV vertex is highlighted by the white dots.

Table 1 :
Three-body neutralino final states induced by the RPV vertices in Eqs.(2.7).The indices i, j, k label the fermion generations.

Table 2 :
Two-body sneutrino final states induced by the RPV vertices in Eqs.(2.7).As usual, i, j, k denote the generation of fermions.