Millicharge or Decay: A Critical Take on Minimal Dark Matter

Minimal Dark Matter (MDM) is a theoretical framework highly appreciated for its minimality and yet its predictivity. Of the two only viable candidates singled out in the original analysis, the scalar eptaplet has been found to decay too quickly to be around today, while the fermionic quintuplet is now being probed by indirect Dark Matter (DM) searches. It is therefore timely to critically review the MDM paradigm, possibly pointing out generalizations of this framework. We propose and explore two distinct directions. One is to abandon the assumption of DM electric neutrality in favor of absolutely stable, millicharged DM candidates which are part of $SU(2)_{\text{L}}$ multiplets with integer isospin. Another possibility is to lower the cutoff of the model, which was originally fixed at the Planck scale, to allow for DM decays. We find new viable MDM candidates and study their phenomenology in detail.

The presence of Dark Matter (DM) in the universe is a clear evidence for new physics beyond the Standard Model (SM). Despite lacking a unique description of DM in terms of elementary particles, a number of general requirements have been identified for DM candidates to fit observations. One of these is stability on cosmological scales. Stability may be explained in terms of symmetries. One may impose a symmetry on a DM model by hand to force stability of the DM candidate, hoping this symmetry can be later justified naturally in ultraviolet completions of the model. Another elegant way to ensure stability is instead through accidental symmetries, the same mechanism that makes the proton stable in the SM. In fact, if one considers only local symmetries as fundamental, other exact or approximate global symmetries can arise as accidental "gifts" given by the specific matter content of the model, which are preserved up to a certain dimension in an effective theory description. This is the main idea behind the "Minimal Dark Matter" (MDM) setup first presented in ref. [1], where the SM is augmented with a new generic multiplet X with mass M and quantum numbers (c, n, Y ) under the SM gauge group SU(3) c × SU(2) L × U(1) Y , without introducing new symmetries. The requirement that the multiplet contains a suitable DM candidate with the correct relic abundance, which is stable on cosmological time scales and is not excluded by present observations, is then used to constrain X 's quantum numbers. For example, the MDM multiplet must be color neutral to avoid the stringent constraints on colored particles [2,3] which seem to exclude the parameter space of thermal relics. Moreover, n must be odd or X 's components would all have sizable tree-level interactions with the Z boson, which are excluded by direct DM searches. The authors of ref. [1] then go on and assume DM electric neutrality, which implies Y = 0. In order to avoid Yukawa couplings with SM fields, as well as dimension-5 effective operators that would cause the DM to decay quickly, it must be n 5 for spin- 1 2 multiplets and n 7 for scalars. Finally, the consistency condition that the theory does not produce a Landau pole below the assumed cut-off at the Planck scale is used to set an upper bound on n, n 5 for Majorana fermions and n 7 for real scalars 1 (these bounds are conservative with respect to those for Dirac fermions and complex scalars). As a result, the authors of ref. [1] single out a fermionic SU(2) L quintuplet and a scalar septuplet as the only viable MDM multiplets.
The eptaplet candidate was recently excluded by the presence of an overlooked dimension-5 operator trilinear in X [4] which makes the DM candidate decay too quickly, while the fermionic quintuplet is seriously constrained by gamma-ray line searches in the Galactic Center [6,7]. In the light of these recent results and of the good sensitivity of present searches to MDM candidates, a critical reanalysis of the MDM framework is timely. Despite the extended literature on the subject and variations thereof (see e.g. refs. [8][9][10][11][12][13][14][15][16][17][18][19][20][21] for some very recent works), some assumptions and basic aspects of the MDM setup remain, that could be more thoroughly examined, and others that could be easily generalized or naturally extended. These are, for instance, the assumption that the cutoff of the theory is at the Planck scale; or the choice of taking small X -Higgs quartic couplings; or even the seemingly natural assumption that the DM is electrically neutral. The aim of this work is to exam-JCAP04(2016)048 ine these aspects in detail, proposing generalizations and studying their phenomenological consequences, in the spirit of Gell-Mann's Totalitarian Principle "everything not forbidden is compulsory". In extending the MDM framework we find a new DM candidate, a Dirac SU(2) L triplet, with a larger degree of compatibility with present bounds with respect to the standard MDM quintuplet. However, the constraints are not necessarily relaxed for the other candidates we propose, meaning that most of the scenarios we study explicitly can be ideally probed by present or near-future experiments.
After critically reviewing the MDM setup and its assumptions in section 2, we explore two main directions in section 3 and 4. In section 3 we abandon the assumption of DM electric neutrality, thus allowing multiplets to get non-zero (although phenomenologically small) hypercharge. These candidates feature a millicharged DM particle which is absolutely stable due to electric charge conservation. This removes the need of a very high cutoff and opens the possibility of large SU(2) L representations without having to worry about Landau poles. We discuss the phenomenology of these candidates and compute the mass needed to achieve thermal production for few of them. Interestingly, the millicharged SU(2) L fermionic triplet is found not (yet) to suffer from the stringent gamma-ray line constraints afflicting the standard fermionic quintuplet. In section 4 we explore the consequences of lowering the cutoff from the Planck scale, so that the MDM fermionic quintuplet decays with observable consequences in the gamma-ray sky. We compute in detail the photon flux (both continuum and line-like features) from DM decays and constrain the cutoff Λ using Fermi data on the diffuse isotropic flux and H.E.S.S. data on gamma-ray lines. Were a clear photon line from this candidate's annihilations to be soon detected, gamma-ray data could also be used to gain insight on the scale of new physics Λ above the DM mass. We conclude in section 5. One interesting technical aspect of our work concerns the presence (or absence) in the Lagrangian of operators of the form X 3 or X 3 H 2 , which cause the DM to decay quickly. We study the issue with the method of Hilbert series in appendix A, while in appendix B we give analytic expressions for the total and differential decay rates of the fermionic MDM quintuplet at dimension 6 in the effective Lagrangian. Finally, we give an analytic treatment of the phase space for 4-body decays into massless particles for our case of interest in appendix C.

Minimal Dark Matter, a critical review
As explained in the Introduction, the MDM setup features the addition of an extra multiplet X to the SM, with quantum numbers (c, n, Y ) under the SM gauge group SU(3) c × SU(2) L × U(1) Y . Further requirements characterizing the MDM framework [1], reported and individually commented below, allow to reduce the number of suitable candidates to a few. For the sake of restricting our discussion to the phenomenologically viable candidates, let us anticipate here some consequences of the requirement that the DM candidate still be allowed by DM searches, see point 4 below.
Given the stringent constraints on colored particles [2,3] which seem to exclude the parameter space of thermal relic DM [1], we restrict ourselves to color-neutral multiplets, c = 1. n must be odd for X to contain a viable DM candidate, with no sizable tree-level interactions with the Z boson and the photon which are excluded by direct DM searches. We do not enforce electric neutrality for the DM at this point, so that Y is allowed to take non-zero (but nevertheless very small, see section 3) values. Fermion multiplets are taken to be vector-like so that they can be given an invariant mass term in the Lagrangian and to cancel anomalies. Notice that the (1, n, Y ) representation with odd n is real for Y = 0, while it is complex for Y = 0.

JCAP04(2016)048
Barring Yukawa interactions of X with two SM fields, which are explicitly forbidden (by gauge invariance) in the MDM setup to avoid DM decay, the renormalizable Lagrangian of the model is for Y = 0, and for Y = 0, where C is the charge conjugation matrix and V (X , H) denotes X 's potential plus possible X -H interaction operators. The lightest state contained in the X multiplet (the DM candidate) is stable under a Z 2 symmetry transforming X → −X for X in a real representation (Y = 0), or a U(1) symmetry transforming X → e iθ X for X in a complex representation (Y = 0). n and Y are chosen so that no renormalizable and dimension-5 interactions exist, that spoil this symmetry thus inducing fast DM decays. This dictates the absence of Yukawa interactions and restricts the operators entering V (X , H).
The MDM setup is characterized by the following requirements, that allow to further reduce the number of suitable candidates [1]. After stating each requirement (written in boldface below), we critically review its implications and comment upon possible generalizations.
1. "The lightest component is automatically stable on cosmological time-scales".
The easiest way to satisfy this condition is probably assuming a very light X , so that the DM particle cannot decay to anything. However, such a multiplet would have been already discovered if charged under strong or weak interactions. The only possibility of having a light DM particle seems therefore for X to have quantum numbers (1, 1, 0) or (1, 1, ) with very small (but positive 2 ) . The former, for a real scalar X , is the well studied scalar singlet DM, deemed to be one of the simplest DM models (see e.g. refs. [22,23] and references therein). UV-complete models of fermionic singlet DM usually also feature a scalar messenger connecting the dark and visible sectors, see e.g. ref. [24], since there exist no renormalizable interactions of a fermionic singlet alone with SM fields (the lowest-order interaction of this DM candidate is through the Higgsportal dimension-5 operator X X H † H, see e.g. refs. [23,25] for recent references). The (1, 1, ) candidate is the so-called millicharged DM. We defer a further examination of this candidate to section 3.
For a heavy multiplet, the stability condition implies that the only acceptable quantum numbers are those for which an accidental symmetry exists, protecting the DM from decay. Assuming the cutoff Λ of the theory to be the Planck scale, as in the original MDM setup, this symmetry must be respected both by renormalizable interactions and by dimension-5 effective operators, in order to avoid fast decay of the DM candidate. Decays induced by higher-dimensional operators violating the accidental symmetry have a negligible impact on the DM phenomenology, but can become important if a lower cutoff is assumed. 2 Here and in the following we take to be a positive number. The case of negative hypercharge is trivially related to that of positive hypercharge.

JCAP04(2016)048
While it is natural to assume a Planck-scale cutoff, new physics is required to explain experimentally observed phenomena like neutrino oscillations or the matter-antimatter asymmetry in the universe, that may occur at a much lower scale. With no further assumption added to the MDM paradigm, one cannot prevent the new physics responsible for these phenomena to break the accidental symmetry stabilizing the DM. It is then natural to study the effect of higher-order symmetry breaking operators and the phenomenology of decaying MDM as a function of the cutoff scale Λ. We will discuss all this for the MDM fermionic quintuplet in section 4. The scalar eptaplet decays too quickly, as discussed in the following, and it is therefore not a good MDM candidate.
When considering DM decays, the most obvious operators to consider are those that are linear in the DM field, which break the accidental symmetry with the minimum number of X fields. However, one should worry about all symmetry-breaking operators, including those with a larger number of X fields. In particular, we show in appendix A that for any SU(2) L n-plet X with integer weak isospin I = (n − 1)/2, three X can be uniquely combined into an SU(2) L singlet for even I (i.e. n = 1, 5, 7, . . . ) or into a triplet for odd I (i.e. n = 3, 9, . . . ). Therefore all scalar (1, n, 0) multiplets with odd n allow for dimension-5 symmetry-breaking operators of the form X 3 H 2 , with H 2 either a SU(2) L singlet or triplet (notice that the Z 2 symmetry protecting scalar DM from decay is already broken at dimension 3 by the operator X 3 for multiplets with even I). Upon closing two of the three X legs in a loop (see e.g. figure 9 of ref. [4]), these operators induce fast DM decays even assuming a Planck-scale cutoff.
The argument just presented is very general and can be also applied outside the MDM framework. In fact, it concerns any model featuring a color and hypercharge-neutral scalar multiplet containing a DM candidate, unless ad-hoc symmetries are introduced to prevent DM from decaying. A possibility to bypass this drawback of scalar DM could be to assign the scalar multiplet a tiny hypercharge. This generalization of the MDM paradigm will be the subject of the next section.
A similar argument as above can be applied to fermion (1, n, 0) multiplets with odd n, for which there exists always a dimension-7 symmetry-breaking operator of the type X 3 LH. Cosmological bounds on the DM lifetime, τ DM > 150 Gyr ≈ 5 × 10 18 s [26,27], then fix a minimum cutoff scale that can be estimated using naive dimensional analysis: implying Λ 10 11 GeV for M ≈ 10 TeV. This bound in turn allows to fix an upper limit on n with the requirement that the model has no Landau poles below this minimal cutoff, n 5 (see footnote 1). As for the scalars, the existence of symmetry-breaking operators relies on the multiplet having zero hypercharge as assumed in the original MDM setup; in the next section we will relax this assumption and show that MDM with Y = 0 can be phenomenologically viable. In line with Gell-Mann's Totalitarian Principle, this condition cannot be satisfied for Lorentz scalars (as already noticed in ref. [1]). In fact, the dimension-4 operators -5 -

JCAP04(2016)048
X † X H † H and X † t a X X H † t a H H with t a the SU(2) L gauge group generators in the proper representation, cannot be forbidden by any choice of X 's quantum numbers. If taken into account, these couplings can affect the annihilation cross section, which is relevant for the computation of the relic abundance and thus for determining M . The operator X † t a X X H † t a H H also affects the mass splitting between the different components of the multiplet. The requirement that the splitting is determined only by loop corrections (see next point) constrains its coefficient to be smaller than O(M/100 TeV) [1]. Moreover, it was found in ref. [5] that the renormalization group evolution of quartic couplings in V (X , H) generates a Landau pole below the Planck scale, and below the Landau pole due to running of the gauge couplings [4], even if the coupling constants are set to zero at the scale M . Although the presence of these operators spoils the minimality of the model by introducing extra free parameters, their effect must be included in any truly generic analysis of scalar MDM candidates.
3. "Quantum corrections generate a mass splitting ∆M such that the lightest component of X is neutral. We compute the value of M for which the thermal relic abundance equals the measured DM abundance". If condition 2 is met, then the splitting can only be radiative, and its size is fixed by X 's quantum numbers. In this case, as shown in ref. [1], the lightest component of X is electrically neutral as long as Y = 0. Even letting the hypercharge take non-zero values, although small enough to be allowed by DM searches as required by condition 4 below, the lightest state is a viable DM candidate albeit electrically charged. Therefore, a lightest state with small |Y | is automatically obtained when conditions 2 and 4 are enforced simultaneously. 4. "The DM candidate is still allowed by DM searches". As anticipated above, thermal-relic colored DM seems to be excluded by the stringent constraints on strongly-interacting DM [1][2][3]. Moreover, constraints from direct DM searches imply that interactions with the photon and the Z boson must be suppressed.
This only leaves open the possibility of (1, n, Y ) MDM with odd n and either Y = 0 or Y = with very small but positive .
Summarizing, we confirmed that, of the two MDM candidates which were so far considered to be viable, the scalar eptaplet is actually ruled out [4]. This singles out the fermionic quintuplet as the only viable MDM candidate. We also proposed to extend the MDM framework by separately dropping two of the original assumptions, namely the assumption of a Planck-scale cutoff and the assumption of DM electric charge neutrality. The assumption of a Planck-scale cutoff can be relaxed in favor of a generic cutoff, which then enters the model as a new free parameter that can be probed by studying cosmic-ray signatures of DM decays. We will do that in section 4. Finally, by lifting the hypothesis of electric neutrality of the DM we established the existence a new class of MDM models featuring millicharged DM. We explore this possibility in section 3. In this section we explore the possibility of MDM candidates with small hypercharge, (1, n, ) with = 0. As in the standard MDM scenario, n must be odd to avoid tree-level interactions of the DM particle with the Z boson. Notice that DM-Higgs interactions can induce -6 -JCAP04(2016)048 a mass splitting in the DM components that makes the DM-nucleus scattering at direct detection experiments inelastic, thus drastically reducing the scattering rate. In this case the stringent bounds from direct DM searches become ineffective and relatively large hypercharge assignments are possible, see refs. [28,29]. We do not pursue this direction here, and stick to small .
An important feature of these candidates is that the DM has electric charge equal to (in units of e), and this makes it absolutely stable. In fact, its stability is protected to all orders in the effective field theory expansion by electric charge conservation. What is usually an unwanted feature in a DM candidate, i.e. electric charge, is here what stabilizes the multiplet making it a potentially successful candidate! Since the DM is stable to all orders, one does not need to worry any more about cutoffs. In the original MDM setup, large multiplets were discarded because the presence of Landau poles in the running of the electroweak gauge couplings could indicate new physics that may spoil the accidental symmetry stabilizing the DM. Millicharged DM being absolutely stable now allows to consider, in principle, even large n's. In this case, a criterion for setting an upper bound on n could be computability. For example we may require that the 1-loop amplitude does not exceed the tree-level result: roughly speaking, (α 2 /4π)G < 1 with α 2 the SU(2) L fine structure constant and G a n-dependent group factor.
Although it may seem odd to consider a field with such a small hypercharge, there is no a priori reason to exclude this possibility: in fact, this choice is allowed by gauge symmetry, and gauge anomaly cancellation is unaffected as long as fermion DM candidates are vector-like.
From a GUT standpoint, one may object that small values of hypercharge are difficult to accommodate in models of grand unification. While this is definitely true, we note that the whole MDM framework is not particularly GUT friendly, since its large multiplets badly modify the running of the SM gauge couplings and moreover they supposedly require a large GUT representation to embed the X field, 3 thus generating a severe doublet-triplet splittinglike problem.
There is also a more theoretical advantage that is worth commenting. According to the no-hair theorem [30,31], gravitational effects break global but not local symmetries. As we observed above, stability of millicharged DM is guaranteed by a local symmetry, the unbroken U(1) EM . Remarkably, this is the only symmetry that could be used to completely stabilize the DM without extending the gauge group of the SM model. This being said, for phenomenological purposes it is enough for a global symmetry stabilizing the DM to be accidentally preserved up to dimension 5 in an effective theory expansion: in fact, the effects of breaking that symmetry at the Planck scale are small enough to guarantee the stability of the DM on cosmological timescales.
In the following we first review the most stringent constraints on the DM electric charge , and then discuss the possible millicharged MDM candidates (1, n, ) and compute the mass of few of them.

Constraints
Constraints on heavy millicharged particles are inferred from cosmological and astrophysical observations as well as direct laboratory tests [32][33][34]. The most stringent upper bounds on , summarized in the following, are shown in the right panel of figure 1 below. A conceivable lower bound could be obtained by considering the weak gravity conjecture [35], which requires > M/M Pl .

Bounds from CMB
Millicharged DM particles scatter off electrons and protons at the recombination epoch via Rutherford-like interactions. It was shown that if millicharged particles couple tightly to the baryon-photon plasma during the recombination epoch, they behave like baryons thus affecting the CMB power spectrum in several ways. This kind of bounds were derived by different groups [32,33]. In particular, ref. [33] found that in order to avoid the tight-coupling condition the DM millicharge must be for a DM particle much heavier than the proton.

Direct searches
Millicharged DM scatters off nuclei via Rutherford-like interactions. In the non-relativistic limit the differential cross section for DM scattering off a nuclear target T with mass m T and electric charge eZ T is given by [36,37] Here E R is the nuclear recoil energy, related to the momentum transfer q by q 2 = 2m T E R , and α is the electromagnetic fine structure constant. F T (q 2 ) is the nuclear Helm form factor [38,39], which takes into account the loss of coherence of the interaction at large q. Since the interaction is spin-independent, the most stringent bound to date is set by the LUX experiment [40]. We use the tools in ref. [41] to infer a 90% CL bound on from LUX. For M 100 GeV, only values 7.6 × 10 −10 M 1 TeV are allowed by LUX with 90% confidence. Notice this bound does not apply in the range because for these values millicharged particles have been evacuated from the galactic disk by supernova explosion shock waves, and galactic magnetic fields prevent them from entering back [32,42]. For respecting eq. (3.1), we do not expect DM self-scattering to sufficiently randomize the direction of motion of DM particles before they gyrate out of the disk [43]. These constraints, depicted in the right panel of figure 1, allow for relatively large values of , which may give rise to interesting phenomenology of millicharged DM candidates. However, for values below the LUX bound this parameter does not contribute to the DM phenomenology and can be safely ignored, the only relevant effect being the doubling of the number of X 's degrees of freedom due to passing from a real to a complex representation of the gauge group.
In the following we discuss the possible millicharged candidates and their phenomenology. We first consider (1, 1, ) candidates, which do not have weak interactions, and then we focus on (1, n, ) with n 3.

(1, 1, ) Dirac fermion
This candidate has only electromagnetic interactions at the renormalizable level. Ref. [32] showed that the parameter space where the DM can be produced thermally with the correct relic abundance is ruled out by observations, most notably by the CMB bound commented above. Therefore, without introducing non-renormalizable interaction, this candidate must be non-thermally produced. Since the details of this production mechanism are highly model dependent, we do not explore this possibility further. Another possibility is to assume a low enough cutoff so that production of this candidate can occur through the dimension-5 Higgs-portal operator X X H † H. In this case, the assumption of thermal production fixes the cutoff as a function of the DM mass, which remains as a free parameter. This candidate has been widely studied in the literature (see e.g. refs. [23,25] for recent references), therefore we do not dwell further on this possibility.

(1, 1, ) complex scalar
With no interactions other than electromagnetic, thermal production of this candidate is ruled out on the same ground as the fermionic (1, 1, ) candidate. However, as already mentioned, scalar DM can interact with the Higgs through the Higgs portal X † X H † H, which opens a new window for thermal production. Given the strong bounds on , this candidate behaves basically as a complex scalar field which is completely neutral under the SM. In the assumption we can neglect the quartic self-coupling (X † X ) 2 , the real and complex components of X do not interact with each other (see below) and can be therefore treated independently as two degenerate real scalar DM particles. Real scalar DM is considered one of the simplest models of DM, and has been widely studied in the literature (see e.g. refs. [22,23] and references therein). The most stringent constraints on thermal relic DM come from the Higgs's invisible decay width [44], which excludes DM masses below ∼ 50 GeV, and the LUX bound [40], which excludes masses from about 10 GeV to roughly 200 GeV [45], except for a very narrow (few GeV wide) interval around half the Higgs boson mass (M ≈ 60 GeV) where the annihilation cross section is resonantly enhanced. As shown in ref. [22], DM masses above 200 GeV will be probed in the near future by both direct and indirect detection experiments, most notably XENON1T [46,47] and CTA [48].

(1, n, ) with n 3
For satisfying the bounds on millicharged DM presented above, DM particles interact mainly with massive gauge bosons. Therefore, the phenomenology of a (1, n, ) multiplet is basically identical to that of (1, n, 0). The only difference, for odd n, is that (1, n, 0) is a real representation of the gauge group while (1, n, ) is complex. This implies a doubling of the number of degrees of freedom with respect to the real case, which affects the computation of the relic density and therefore the DM mass M . In the following we show that, under some conditions, the relic density for (1, n, ) can be obtained by simply scaling results for (1, n, 0) appeared in the literature.
We start by expressing a complex scalar multiplet in terms of its real components, and a Dirac fermion as two degenerate Majorana states with opposite parity under charge conjugation: for fermion X . In this new basis, and in the absence of quartic couplings in V (X , H) for scalar MDM, we get two separate Lagrangians for X 1 and X 2 that are bilinear in these fields, and in fact we can define a global Z 2 symmetry acting on (X 1 , X 2 ) as If we consider now the main annihilation mode of this candidate, i.e. 2 → 2 DM annihilations into SM vectors V , this symmetry tells us that the only possibile annihilation channels are has initial and final states with different parity. This means that, at tree level, the two sectors completely decouple. This is not the whole story, however. As well known, Sommerfeld enhancement provides an important non-perturbative correction relevant in the non-relativistic regime, and thus must be included. Since all relevant diagrams are of ladder type, we have again that if a process is initiated e.g. by X i X i then no X j =i particles appear in the diagram. Therefore Sommerfeld enhancement respects the complete factorization of X 1 and X 2 . As a consequence, the computation of the relic density of X 1 is completely independent from the computation of the relic density of X 2 . Moreover, these two states having same mass and same gauge interactions, they must have the same relic density and therefore the relic density of X is twice that of a single X i . Figure 1 shows the Sommerfeld-corrected DM relic density as a function of the DM mass for a complex scalar triplet and eptaplet and a Dirac triplet and quintuplet (solid lines). These functions are taken to be twice the value of the relic density of a real scalar triplet [49], a real scalar eptaplet in the approximate SU(2) L -symmetric limit [11], and a Majorana triplet [49] and quintuplet [6], respectively (the relic density of the Majorana quintuplet is also shown as a dashed line). Since the real scalar quintuplet was found in ref. [49] to have the same mass of the Majorana quintuplet, we assume the same holds also for the complex case for both the quintuplet and the eptaplet. While candidates with larger n can be perfectly viable, we only consider here SU(2) L triplets, quintuplets and eptaplets as careful computations of the relic abundance are available in the literature only for these candidates. The horizontal red strips in figure 1 show Planck's measurement of DM density [50], at the 1σ (inner strip) and 2σ (outer strip) CL. The DM mass for each case is determined by the crossing of the relevant solid line and the red strips (notice the relic density line for the Dirac triplet crosses the DM abundance band twice, thus there are two allowed values for its mass). This interval in DM masses is indicated with a vertical band whose width is given by the 2σ uncertainty of Planck's result. A larger uncertainty, of order 5% of the total result, comes however from the theoretical determination of the cross sections [6]. This translates into an uncertainty on the determination of the DM mass shown as a lighter vertical band for each case. Considering the latter uncertainty the DM has mass 1.55 ± 0.08 TeV complex scalar triplet, Also shown in figure 1 is the mass for a Majorana quintuplet (which is determined by its relic density, shown as a dashed line). This information will be useful for our study of this candidate in the next section.  For scalar X , quartic couplings such as X † t a X X H † t a H H can break the symmetry in eq. (3.6) and thus affect the above scaling argument. In this case, the annihilation cross section will be in general larger than that, discussed above, due solely to DM couplings to gauge bosons. Therefore, in order to fit the observed DM abundance, the DM mass must be larger with respect to the case of DM with only gauge interactions. The values of M given above and shown in figure 1 can thus be thought of lower bounds on the true value of the DM mass in presence of quartic couplings. See e.g. ref. [51] for a dedicated analysis on the effect of these couplings.
Once the mass of the (1, n, ) multiplet is known, the phenomenology of these candidates is univocally determined (up to free terms in V (X , H) for scalar multiplets). In particular, the most stringent constraints on electroweak multiplets come from indirect DM searches. The bounds from gamma-ray line searches are particularly relevant due to the Sommerfeldenhanced annihilation cross section into gauge bosons.
The phenomenological advantage of millicharged MDM candidates is that, since the DM particle and its antiparticle are distinct, the annihilation probability is half that of a selfconjugated DM candidate with the same quantum numbers and mass. Therefore, all bounds on the annihilation cross section are a factor of 2 less stringent. However, the DM mass for these candidates is in general lower than for their self-conjugated version, and this may be a -11 -JCAP04(2016)048 drawback for the following reason. Ideally, bounds on annihilation cross sections scale with the inverse of the DM number density squared, and thus with (ρ/M ) −2 with ρ the assumed DM energy density. Therefore, lighter DM candidates are ideally more constrained. However, realistic bounds depend on the experimental resolution, which is particularly relevant for gamma-ray line searches given that the expected signal is a very narrow spectral feature. A finite and energy-dependent resolution leads in general to an uneven sensitivity on the DM mass. Moreover, the theoretical dependence of the annihilation cross section on M may be very irregular, especially in presence of non-perturbative effects (see e.g. figure 7 of ref. [6]). Therefore, it is difficult to predict whether a lighter DM particle is more or less constrained than a somehow heavier particle. For this reason, constraints on millicharged MDM candidates must be checked case by case.
Bounds on some of the candidates considered above can be determined by properly rescaling existing bounds on self-conjugated multiplets with the same quantum numbers. Constraints on a (supersymmetric Wino) Majorana triplet, on the MDM Majorana quintuplet, and on the real scalar eptaplet can be found in refs. [6,7,49,[52][53][54][55][56], and [11], respectively. We do not have enough information on the scalar triplet and fermion eptaplet to determine bounds on these candidates.
Interestingly, the Dirac triplet with M = 2.00 TeV is allowed by gamma-ray searches even with the most aggressive choices of DM profile made in figure 12 of ref. [52]. In the assumption of a cuspy profile, forthcoming experiments like CTA [48] will be able to probe this candidate. The situation of the Dirac triplet with M = 2.45 TeV is closer to (although worse than) that of the Majorana triplet with mass 3.1 TeV [53], which is already excluded by bounds assuming cuspy profiles while allowed when choosing a cored profile. The 6.55 TeV Dirac quintuplet is in the same situation as the Majorana quintuplet, whose mass is given in eq. (4.4), i.e. it is badly excluded with the choice of a cuspy profile, while it is still viable if a cored profile is considered (see e.g. figure 7 of ref. [6]). The complex scalar eptaplet, while excluded for a cuspy Einasto profile, may be either excluded or allowed for a cored Isothermal profile, depending on the precise value of its mass within the 5% uncertainty reported above and shown in figure 1; notice however that our calculation of mass and constraints for this candidate rely on the computations carried out in ref. [11] in the approximate limit of exact SU(2) L symmetry.
We have only discussed here bounds from gamma-ray line searches, which are the most constraining, as mentioned above, when a cuspy profile is chosen. Other bounds, that a rough evaluation reveals not to exclude these candidates at present, may become relevant in the near future (see e.g. refs. [6,52,57]). The most entertaining possibility is to probe MDM with a future 100 TeV proton-proton collider [57], but the hope is to find other evidences for it well before that.

Decaying quintuplet MDM
In this section we study the possibility that the MDM setup has a generic cutoff Λ. We consider here the 'standard' MDM scenario with Y = 0, thus the only viable candidate (as discussed above) is the fermionic SU(2) L quintuplet. The main effect of lowering the cutoff from the Planck scale is that DM stability on cosmological timescales is spoiled. In fact, dimension-6 operators can break the Z 2 symmetry protecting DM from decay [58], so that for a low-enough cutoff we can expect to observe the signature of DM decays in cosmic-ray spectra. We thus perform a thorough analysis of the gamma-ray spectrum produced in DM decays and use the Fermi [59] and H.E.S.S. [60] data to set bounds on Λ.

Relevant Lagrangian and decay modes
We represent the fermionic SU(2) L quintuplet as a Dirac four-spinor X with only righthanded components, so that P R X = X and P L X = 0, where P R and P L are the right and left projectors, respectively. Notice that the quintuplet is a real representation of SU(2) and therefore the neutral component of X is a Majorana fermion. DM decays are induced at dimension 6 by the two operators X LHHH † , X σ µν LW µν H, and their hermitian conjugates. We are therefore interested in the following Lagrangian: where a = e, µ, τ is a lepton flavor index and σ µν ≡ i 2 [γ µ , γ ν ]. We neglect dimension-5 and all other dimension-6 operators, since they do not contribute to DM decays. DM annihilations are of course dominated by X 's renormalizable gauge couplings, as the contribution of non-renormalizable operators is suppressed by powers of M/Λ. To show how the multiplet components of the various fields contract, we represent X as a rank-4 completely symmetric tensor in the anti-fundamental representation of SU(2) L , X ijkl with i, j, k, l = 1, 2 (see e.g. appendix B of ref. [4] for more details). The W -boson multiplet is also written as a symmetric rank-2 tensor, while the lepton doublet L and the Higgs doublet H are represented by rank-1 tensors in the fundamental representation. Indices are raised and lowered with the completely antisymmetric SU(2) L -invariant tensor , with 12 = − 12 = 1. Making the SU(2) L indices explicit we get The fields can be rewritten as where the X t 3 R are Dirac spinors with only right-handed components. The DM candidate is the self-conjugated neutral component X 0 R , which from now on we will denote X 0 for simplicity. As can be seen from figure 1, the DM mass is fixed by its relic abundance to be M = 9.4 ± 0.47 TeV . (4.4) We study the effect of the two dimension-6 operators separately. Detailed analytic formulas for the matrix elements and the differential and total decay rates for the relevant processes can be found in appendix B. Since the DM is much heavier than all SM particles it decays to, we neglect all final state particle masses in the calculations.

JCAP04(2016)048
The first operator, X LHHH † , induces ν L -X 0 mixing and the following 2, 3, and 4-body DM decays: where all gauge bosons are longitudinal and charged leptons and neutrinos are left-handed (the polarization of final state particles is an important ingredient entering the code [61], that we use to compute the gamma-ray flux from DM decays). Notice that the squared amplitude for decays into final states with many Higgs fields are enhanced by powers of (M/v) 2 ≈ 10 3 : in fact, adding a Higgs field to the final state removes a factor of v from the Lagrangian coefficient, that is replaced in the decay amplitude by a factor of M . By the Equivalence Theorem, the same holds for decays into many longitudinal gauge bosons as well. Therefore, despite the phase-space suppression, decays into many particles can be favored over 2-body decays. We check explicitly that this is the case by computing both X 0 → νh, νhh, and νhhh decay rates. The remaining decay rates are computed with the Equivalence Theorem. All our analytic results can be found in appendix B. Our analytic computation of the 4-body phase space (approximating all final states as massless) is described in appendix C. The second operator, X σ µν LW µν H, induces the following 2 and 3-body decays: where one gauge boson is always transverse, while the other, if present, is longitudinal, and the charged lepton or neutrino is left-handed. Contrary to the previous case, 4-body decays do not receive the (M/v) 2 enhancement factor with respect to the 3-body modes, and thus can be neglected. 2-body decays are also suppressed with respect to the 3-body channels, but they deserve special attention since they produce very narrow features in the gamma-ray spectrum. These peaks, appearing at the very end of the produced photon spectrum (i.e. at energies equal to half the DM mass), are due to photons produced by the monochromatic decay products of 2-body decays, and may be visible on top of the continuum produced by 3-body decays. In appendix B we compute explicitly the X 0 → W, νZ, νγ, W h, νZh, νγh decay rates, and apply the Equivalence Theorem to compute the remaining rates. To avoid the shortcomings of the Theorem (see e.g. footnote 7 of [62]) we checked the result by also performing the computation in the Equivalent gauge [62].

Gamma-ray fluxes from DM decay
The strongest limit on models of decaying DM is arguably set by observations of gamma rays. For this reason we focus here on production of secondary gamma rays, and compare the model expectation with Fermi data to obtain a bound on the relevant parameter c a 1,2 /Λ 2 . Moreover, the photon flux does not suffer from the same astrophysical (e.g. diffusive) uncertainties as charged particles, thus our analysis is quite reliable and does not depend much on the modeling of the cosmic environment. Other relevant constraints may be set for instance by looking at the sum of electron and positron fluxes up to 1 TeV measured by

Production rate of stable SM particles in DM decays
Upon decay of the DM particle, X 0 → f 1 . . . f n , the primary decay products f undergo a series of processes such as decay and radiative processes (hadronization, showering, . . .) which generate a set of stable SM particles α = e ± ,p, γ, . . . . The production rate of each stable state α at the source per single DM decay is where dN f α /dE α is the spectrum arising solely from the primary f with energy E f , and dΓ/dE f is the DM decay rate into f summed over all decay channels which include f in the final state.
While propagating away from the source, these stable particles can interact with the cosmic environment thus modifying their spectrum in a position-dependent way. For instance, the photon flux at Earth gets a contribution from the prompt emission in eq. (4.7) with α = γ, and a contribution from low-energy background photons (e.g. from the CMB or the interstellar photon field) being up-scattered by e ± from DM decays; in the latter case, the α = e ± rate at the source in eq. (4.7) must be convolved with the probability of undergoing inverse-Compton (IC) processes with the inhomogeneous photon field (see e.g. refs. [64,65]). The so-modified rate, which we call dR α /dE α , depends e.g. on the distance r from the Galactic Center (GC) for decays within our galaxy, or from the redshift z for extragalactic decays. We compute dR α /dE α from the production rate at the source dR s α /dE α following ref. [61]. We adapt the spectra per single primary dN f α /dE α from ref. [61], with the following cautions. The primary spectra given in ref. [61] are meant for DM decays to particle-antiparticle pairs X 0 → f f , so that the primary energies E f are not parameters that can be varied but are instead fixed to half the input DM mass, call it M PPPC . The latter is a parameter whose value is possible to vary, and therefore in using the primary spectra from ref. [61] we adopt the prescription M PPPC → 2E f . One has also to take into account the fact that the primary spectra given in ref. [61] include the spectrum generated by the primary antiparticle f , besides that due to f . However, in the assumption of CP invariance, the rate for the decay X 0 → f 1 . . . f n equals the rate for X 0 → f 1 . . . f n (notice that X 0 is a Majorana fermion). Therefore, when summing the two rates (as part of the sum over all decay channels), we will have in eq. (4.7) dN f α /dE α + dN f α /dE α . We use the spectra of ref. [61] in place of this sum. Consequently, the only channels that remain to be summed are those that are non mutually conjugated. In practice eq. (4.7) can be operatively written as

Continuum photon emission
The residual isotropic gamma-ray flux observed by Fermi [59] extends from 100 MeV up to 820 GeV. The origin of this isotropic emissions is not well understood and can be due to different phenomena such as unresolved sources or truly diffusive processes. DM decays can contribute to this isotropic flux, with two components: i) galactic (Gal) residual flux due to DM decays within the Milky Way halo, and ii) extragalactic (ExGal) flux due to cosmological DM decays integrated over redshift. The latter is of course isotropic, while the former is not, however its minimum constitutes an irreducible contribution to the isotropic flux. Therefore, the isotropic diffuse gamma-ray flux as measured by Fermi is Here we make the reasonable approximation that the minimum of the angular flux in the galaxy is found at the anti-GC, as in refs. [66][67][68]. This approximation is well justified because, for the decay channels relevant in our analysis, the prompt flux (which is the dominant contribution) follows the angular distribution of DM density, which is of course minimum at the anti-GC. For reasonable sizes of the diffusive halo, moreover, we expect also the IC contribution to approximately follow the angular DM distribution. The flux from the galactic halo, observed from a given direction and within a solid angle dΩ, is in general given by with parameters ρ s = 0.184 GeV/cm 3 and r s = 24.42 kpc. In any case, given the linear dependence of the photon flux on the DM density for decaying DM (as opposed to the quadratic dependence for annihilating DM), and the fact that we are mainly interested in the anti-GC where all profiles are similar, the final result will bear little dependence on this choice of profile and parameters. The galactic gamma-ray flux has two main components: i) the prompt gamma rays originating from the fragmentation of the primary products of decay, whose spectrum can be obtained by taking α = γ in eq. (4.7) for all decay channels given above; and ii) IC gamma rays, produced by the up-scattering of low-energy photons of CMB, infrared light, and starlight, by energetic e ± produced by DM decays. To obtain the IC spectrum we integrate the α = e ± rate in eq. (4.7) with the IC halo functions given in refs. [61,70].
The extragalactic flux, integrated over the redshift at emission z, is given by -16 -

JCAP04(2016)048
where H 0 is today's Hubble expansion rate and Ω M , Ω DM , Ω Λ are respectively the matter, DM, and cosmological constant energy density in terms of today's critical density ρ c,0 . E γ = (1 + z)E γ is the photon energy at emission redshift z so that the same photon is detected on Earth with energy E γ . The factor e −τ (z,Eγ ) accounts for the absorption of DM-produced gamma rays due to scattering off low-energy background photons, which results in production of energetic e ± pairs. The converted energy is in turn redistributed to lower-energy gamma rays via IC scattering off CMB photons. We take into account this effect in our analysis, which is sizable for channels with a pronounced prompt emission (the most relevant cases being X 0 → νγ, eW + , µW + ). We take the optical depth of the Universe τ (z, E γ ) from ref. [61]. As for the galactic flux, the extragalactic spectrum is again the sum of the prompt and IC contributions. The z dependence of the prompt flux is obtained by simply "redshifting" E γ , while the IC contribution due to e ± scattering off the warmer CMB photons at z > 0 scales as [66] dR IC Clearly the DM signal does not agree in shape with the data, which are instead well fitted by a simple power law with an exponential cutoff. We therefore use this functional form to model the background, adopting the best-fit parameter values from table 4 of ref. [59] and the Fermi baseline diffuse galactic emission model (model A of ref. [59]). Each plot reports the value of used to normalize the photon flux from DM decays, which has been determined in each case by performing a minimum-χ 2 analysis of model + background against the Fermi data. As apparent from a comparison of the three figures for each operator, the total DM signal for each operator bears little sensitivity to the lepton flavor a, the most noticeable difference being the size of the broad bump peaking between 10 and 100 GeV. Since this ����� ���� ����� ���� ��� ���� Figure 2. Isotropic gamma-ray flux due to DM decays induced by the operators (Λ a 1 ) −2 X L a HHH † (left) and (Λ a 2 ) −2 X σ µν L a W µν H (right), assuming DM coupling to electrons and electron neutrinos (a = e). Fluxes from 2, 3 and 4-body decays are separately shown in red, green, and blue, respectively, while the total flux is in black. Dashed lines indicate the extra-galactic component of the flux, dotted lines the galactic flux, and solid lines their sum. Fermi data on the diffuse isotropic gamma-ray flux are shown in brown, and the astrophysical background is displayed as a gray line. The thick black line indicates the sum of the total flux from DM decays and the background. The best-fitting value of Λ a 1,2 , adopted here to normalize the fluxes, is reported in the upper part of the plots. bump is due to IC processes that populate the high-energy gamma-ray spectrum at the expense of high-energy e ± , DM decays generating more e ± are expected to make it larger. For this reason, the bump is largest for DM coupling to a = e, somehow smaller for a = µ, and smallest for a = τ since τ 's mainly decay into hadrons. The main difference between the photon fluxes of the two operators is the narrow feature at the very high-energy end of the spectrum for X σ µν LW µν H, which is absent for the other operator. This is due to the prompt photon emission, especially from the 2-body X 0 → νγ decay which generates a gamma-ray line at E γ = M/2. 3-body decays with γ's in their final states also contribute to the peak, although with a broader spectrum. We will analyze this feature in more detail in the next section, where we derive a complementary bound coming from line-like searches with the H.E.S.S. telescope.
In deriving a bound on the maximum allowed DM signal, i.e. on the minimum value of Λ a i allowed by data, we adopt two methods.
• DM signal only. This method, which yields a very conservative limit, consists in demanding that the gamma-ray flux from DM decays alone, i.e. assuming no background, does not exceed any one of the Fermi data points by more than a given significance, which we take to be 3σ. This option is largely conservative for two main reasons. First, it is quite clear from figures 2, 3, 4, that allowing the flux to exceed one data point would result in excesses in nearby data points as well, and therefore the global significance of the exclusion is in principle higher than the chosen significance in one single bin. This is due to the smooth nature of the DM signal (and of the data) in the Fermi energy window. Second, the assumption of a negligible background is clearly physically untenable. In fact, the spectral shape of the signal is so different from the data that background is needed in order to obtain a good fit.
• DM signal + background. A more realistic method consists in demanding that the sum of astrophysical background and DM signal does not exceed a chosen level of global significance, which we take to be 3σ.
The values of Λ a i allowed by Fermi data as computed with both methods are summarized in table 1, together with the respective bounds on the DM lifetime τ DM . The constraints on τ DM can be compared for a reference with the bounds obtained in refs. [67,72] from an -19 -

JCAP04(2016)048
earlier Fermi data release. The enhanced constraining power of the new data set is apparent. One also needs to bear in mind that, compared to the phenomenological analysis carried out in ref. [67], where only one decay mode is present at a time, the DM candidate considered here features several decay modes some of which contributing negligibly to the gamma-ray emission. For this reason, our bounds on τ DM may appear less strong than naively expected.

Gamma-ray lines
As commented above, the gamma-ray flux from DM decays induced by the X σ µν LW µν H operator displays a sharp feature at energies close to M/2 (see right panels of figures 2, 3, and 4). This is due to the presence of decay channels with prompt photon emission, most notably the X 0 → νγ decay which generates a gamma-ray line at E γ = M/2 ≈ 5 TeV. This feature is not constrained by Fermi measurements, which only reach up to 820 GeV. Therefore we compute here a bound from H.E.S.S. gamma-ray line searches [60], which extend up to 25 TeV.
The H.E.S.S. Collaboration performed two separate searches for line-like features in the gamma-ray flux in two sky regions of interest, namely the extragalactic sky and the central galactic halo (CGH) region, the latter defined as a circle of 1 • radius around the GC, where the Galactic plane is excluded by requiring |b| > 0.3 • . We compare our gamma-ray flux with H.E.S.S. limits in both sky regions, thus producing two sets of bounds. To take into account the finite experimental resolution we convolve the photon flux with a Gaussian G(E γ , E) centered around E γ , where E denotes the energy detected by the instrument. We take the Gaussian function to have resolution 15% of E γ [60,73]. We then integrate the signal over each bin in detected energy, 16) and compare the result with the 95% CL limits on the gamma-ray flux in both sky regions shown in figure 2 of ref. [60]. We neglect IC processes as they only contribute to the continuum gamma-ray spectrum, not to line-like features, thus we compute the photon flux only using the position-independent prompt emission in eq. (4.7) with α = γ. Contrary to the extragalactic flux, the flux in the CGH region is sensitive to the assumed DM density profile due to the pronounced differences between cored and cuspy profiles close to the GC. For this reason we use the H.E.S.S. bound in the CGH region to set a profile-independent bound on Λ a i /J 1/4 with ds r ρ halo (r(s, ψ(b, ))) ρ (4.17) the angular-averaged J factor in the sky region of interest. Notice that this bound is truly profile-independent as long as position-dependent processes such as IC can be neglected. For reference, the value ofJ for a cuspy profile like NFW [69] and a cored profile like Burkert [74] computed with the functions in ref. [ The Fermi limits from the continuum gamma-ray flux prove to be stronger than the H.E.S.S. bounds from gamma-ray lines (as also found by ref. [73]), up to a factor of 8.7 for the limit on τ DM almost independently on the DM profile in the galaxy.  Table 1. Gamma-ray bounds on the new-physics scale Λ a 1,2 defined in eq. (4.15) and on the DM lifetime τ DM , separately for the two operators (Λ a 1 ) −2 X L a HHH † and (Λ a 2 ) −2 X σ µν L a W µν H and for each lepton flavor a = e, µ, τ . Both operators are constrained by the Fermi measurement of the isotropic diffuse gamma-ray flux, which is used here to derive a conservative bound considering the DM signal alone, and a realistic bound considering DM signal + background. The dipole operator induces a gamma-ray line-like feature in the photon spectrum, and thus is also constrained by H.E.S.S. searches of gamma-ray lines in the CGH and extragalactic regions. The DM-profile dependence of the bounds from the CGH region is factored in theJ factor, values for which are given in eq. (4.18) for two example density profiles. All other bounds are reasonably independent of the DM profile in the halo. current sensitivity on the partial decay width into channels with prompt photons is at the level of Γ −1 γ ∼ 10 28 s [73], our bounds on τ DM 10 27 s are less stringent due to the fact we include all relevant decay channels. In other words, we constrain the full DM decay width rather than the partial width into channels with prompt photons, by also taking into account the important contribution of other decay modes.
The bounds shown in table 1 were derived separately for the two operators X LHHH † and X σ µν LW µν H, assuming only one was turned on at a time. However, in general, both operators are expected to arise in the effective theory description, and, if they are generated by the same physics at the scale Λ, we also expect their coefficients c a 1 and c a 2 to be somehow related. Since the dipole operator X σ µν LW µν H is certainly generated at loop level, while X LHHH † can conceivably originate at tree level, we can guess that c a 2 ≈ (α 2 /4π)c a 1 with α 2 ≈ 1/25 as expected from the renormalization group evolution of the weak coupling if the new physics in the loop is at the GUT scale. Therefore, the prospects of detecting the gamma-ray line-like feature originating from the dipole operator are much worse than naively expected from the study of the operator alone. Figure 5 shows the gamma-ray flux due to DM decays induced by the operators 1 Λ 2 X LHHH † (red line) and α 2 4πΛ 2 X σ µν LW µν H (green line), and their sum (black line). It is clear that the resulting line-like feature is much less visible against the continuum of photons than in our previous analysis considering just one operator. This result shows that analyses of gamma-ray line-like signatures of specific operators within an effective theory description should be accompanied by an assessment of the contribution to the continuum photon flux of all other operators that are expected in the effective theory.  Figure 5. Isotropic gamma-ray flux due to DM decays induced by the operators (Λ a ) −2 X L a HHH † (red line) and (α 2 /4π)(Λ a ) −2 X σ µν L a W µν H (green line), and their sum (black line). The suppression factor for the dipole coupling is expected from a radiative nature of the operator and causes the gamma-ray line-like signal to be dwarfed by the continuum photon flux. The three plots assume DM coupling to a = e (left), a = µ (center ), and a = τ (right). The best-fitting value of Λ a , adopted here to normalize the fluxes, is reported in the upper part of each plot.

Conclusions
Minimal Dark Matter (MDM) [1] is a theoretical framework highly appreciated for its minimality and yet its predictivity. Contrary to many models where DM stability is imposed by hand through a global symmetry, MDM candidates are made stable on cosmological timescales by accidental symmetries occurring through a careful selection of the DM quantum numbers. When the cutoff of the model is taken to be the Planck scale, internal consistency conditions (the absence of Landau poles below the cutoff scale) and phenomenological constraints single out a fermionic SU(2) L quintuplet and a scalar eptaplet as the only viable MDM candidates.
Recently, the MDM model was endangered by the discovery that the eptaplet decays quickly due to a previously overlooked dimension-5 operator [4], and thus it is not a viable candidate; and by stringent gamma-ray line constraints in the Galactic Center, which do or almost do rule out the quintuplet, depending on the assumed DM density profile in the halo [6,7]. In the light of these recent results, a critical reanalysis of MDM aiming at generalizing and extending this framework was in order. This is the purpose of the present paper. After reviewing the MDM setup and its assumptions, we proposed two possible generalizations and studied their phenomenological implications. First, we found that MDM multiplets with a small enough hypercharge provide viable DM candidates, which possess small electric charges (the so-called millicharged DM) and are therefore absolutely stable. We discussed the case of millicharged singlet DM, and determined the thermal relic of triplets, quintuplets and eptaplets thus obtaining their mass. Interestingly, we found that a Dirac triplet is not constrained by the gamma-ray line searches that, for a cuspy DM halo profile, rule out a Wino (Majorana triplet) and the original MDM quintuplet.
Second, we proposed the possibility of lowering the Planck-scale cutoff for the original model of MDM quintuplet with zero hypercharge. As a consequence, the DM can decay by means of two dimension-6 operators which break the accidental symmetry and we can observe the signature of these decays in the gamma-ray sky. We found the cutoff to be constrained by Fermi data on the diffuse isotropic gamma-ray flux at about the GUT scale. We also discussed the constraints set by H.E.S.S. on the gamma-ray line-like feature produced by the dipole operator, finding that the Fermi data set a stronger bound for a 10 TeV DM. We also found that, when the dipole operator is assumed to be generated by loop processes, this linelike feature is completely dwarfed by the photon continuum induced by the other operator. Were a clear photon line from this candidate's annihilations to be soon detected, gamma-ray data could also be used to gain insight on the scale of new physics above the DM mass.
B MDM quintuplet decay rates B.1 X LHHH † As explained in section 4, this operator induces the DM decay modes listed in eq. (4.5). We give here detailed analytic expressions for these decay rates. We compute explicitly the X 0 → νh, νhh, νhhh decay rates and then derive all other rates by applying the Equivalence Theorem.
Our analytic computation of the 4-body phase space (approximating all final states as massless) is described in appendix C. When computing decay rates into final states with n identical particles, we consider the n! identical diagrams contributing to the scattering amplitude, and the 1/n! phase space reduction factor to prevent double-counting identical configurations.
DM decays into a neutrino plus Higgses are given by with v = 246 GeV. The relevant polarization sum entering the spin-averaged squared matrix element is 1 2 s,r |ū s (p ν )P R u r (p X )| 2 = p ν · p X = M E ν , (B.2) where in the last equality we set ourselves in the rest frame of the decaying particle.

JCAP04(2016)048
These values take into account the appropriate n! factors in the decay rate due to the presence of indistinguishable particles in the final state, both on the left and right-hand side of eq. (B.13). As per assumptions of the Equivalence Theorem, these relations hold in the high-energy limit where the DM particle is much heavier than its decay products, which is true in this case given eq. (4.4). Only longitudinal gauge bosons contribute significantly to the rate in this limit.
B.2 X σ µν LW µν H From the second operator we get two terms inducing X 0 decay: The most relevant decays induced by this operator are X 0 → W, νZ, νγ, W h, νZh, νγh (see section 4). The relevant polarization sum entering the spin-averaged squared matrix element of the processes X 0 → f V (h) (with f a fermion and V a vector boson) is 1 2 q s,r ū s (p f )σ µν P R u r (p X )p µ V ε ν * q (p V ) 2 = 4(p X · p V )(p f · p V ) . (B.22)