Gamma ray tests of Minimal Dark Matter

We reconsider the model of Minimal Dark Matter (a fermionic, hypercharge-less quintuplet of the EW interactions) and compute its gamma ray signatures. We compare them with a number of gamma ray probes: the galactic halo diffuse measurements, the galactic center line searches and recent dwarf galaxies observations. We find that the original minimal model, whose mass is fixed at 9.4 TeV by the relic abundance requirement, is constrained by the line searches from the Galactic Center: it is ruled out if the Milky Way possesses a cuspy profile such as NFW but it is still allowed if it has a cored one. Observations of dwarf spheroidal galaxies are also relevant (in particular searches for lines), and ongoing astrophysical progresses on these systems have the potential to eventually rule out the model. We also explore a wider mass range, which applies to the case in which the relic abundance requirement is relaxed. Most of our results can be safely extended to the larger class of multi-TeV WIMP DM annihilating into massive gauge bosons.


Introduction
The Large Hadron Collider (LHC) has completed the Standard Model (SM), by unveiling the existence of the higgs boson with a mass of 125.1 GeV. So far, however, the LHC has found no convincing evidence whatsoever for the long-sought-after New Physics at the TeV-scale, expected to be responsible for keeping the mass of the higgs boson so light. Hence, doubts are cast on the very relevance of the naturalness idea and in particular on the special role of such a scale in modern particle physics.
On the other hand, thanks to the argument that goes under the name of WIMP miracle (Weak Interacting Massive Particle), the TeV-scale does remain appealing in Dark Matter (DM) terms, independently of its possible role or not for the naturalness problem. Indeed, a particle in that mass range and with Weak SM interactions can naturally provide the required DM abundance observed in cosmology via the thermal freeze-out mechanism, which is not a small feat. While this may eventually turn out to be a mirage, it still constitutes a very well motivated driving principle in the quest for the nature of the new particle that has to exist to constitute the observed DM abundance.
Indeed, betting on the possibility of a pure WIMP nature of Dark Matter and of the SM being the ultimate theory all the way up to the GUT or Planck scale (which is admittedly a bold set of assumptions), the model of Minimal Dark Matter (MDM) was proposed back in 2005 [1,2]. In a nutshell, the construction proposes to add to the SM the minimal amount of new physics (just one extra EW multiplet χ) and search for the minimal assignments of its quantum numbers (spin, isospin and hypercharge) that make it a good DM candidate without ruining the positive features of the SM. No ad hoc extra features are introduced: the stability of the successful candidates is guaranteed by the SM gauge symmetry and by -1 -JCAP10(2015)026 renormalizability. By following these principles of consistency to their end and applying the most evident phenomenological constraints (e.g. direct detection bounds from scattering on nuclei), the theory selects a fermionic SU(2) L 5-plet with null hypercharge as the only one which provides, in terms of its electrically neutral component, a viable DM candidate. 1 The annihilation cross sections of such particles can be fully computed in EW theory, including the Sommerfeld enhancement. Thus, by solving the Boltzmann equation for the relic abundance and requiring such abundance to match the measured cosmological value, the mass of the DM candidate can be univocally determined, turning out to be around 9.4 TeV (we will return in detail on these aspects). The components of the multiplet are then split by 1-loop electroweak corrections, which produce small differences. 2 This implies that the charged components χ ± and χ ±± disappear from the thermal bath by decaying, with lifetime O(ns), into the lightest one χ 0 , which constitutes the cosmological DM. Due to its minimality, the theory is remarkably predictive: no free parameters are present and therefore the phenomenological signatures can be univocally calculated.
Almost ten years later, and in the context described above, it makes sense to test whether the MDM construction still stands the comparison with data. Broadly speaking: its large predicted mass makes its production cross section at the LHC extremely suppressed [6], and even prospects for future colliders are bleak (see [7,8] for related studies); scattering on nuclei for direct detection is much suppressed with respect to initial estimates, as pointed out in a series of recent works (see [9] and references therein); indirect detection (ID), instead, initially considered in [10,11], remains promising. Among the ID messengers, we focus in this work on gamma rays, and we will comment in the end on other possible channels.
Beyond the motivations specific to the MDM construction, the scale for DM is being pushed towards the TeV and multi-TeV range by null results at the LHC and in direct and indirect detection experiments, as we mentioned above. It seems therefore natural to explore that regime. We will indeed work in a broad range of masses, spanning from 100 GeV to 30 TeV. More precisely, we are concerned with three levels of generality in this analysis: The specific original MDM 5plet candidate: it has fully determined mass (9.4 TeV) and annihilation cross sections. We will refer to this as 'the MDM 5plet' and point to it with a vertical band in our constraint plots (figure 4, 6 and 7). 1 For such a DM candidate no gauge invariant effective operator, made of a DM field and SM fields, hence mediating its decay, can be written with dimension less than 6. A dimension 6 operator induced by the exchange of a GUT scale (or higher scale) heavy particle leads to a large enough DM lifetime, i.e. not ruled out by the constraints which hold on the fluxes of cosmic rays that such a decay scenario would induce. Another originally considered candidate, the scalar 7-plet, is now ruled out by higher dimensional operators (initially overlooked) that mediate its fast decay [3,4]. In addition, this candidate develops a Landau pole at a low scale [5]. 2 The explicit lagrangian of the model, for reference, reads: where χ is the fermion 5-plet Minimal Dark Matter candidate, g is the SU(2) gauge coupling, and sw and cw are the sine and cosine of the Weinberg angle. Mχ is the degenerate mass of the multiplet and M χ 0 , M χ ± = M χ 0 + ∆M and M χ ±± = M χ 0 + 4∆M , with ∆M = 166 MeV, are those of the individual components after splitting.

JCAP10(2015)026
A 'generic' 5plet DM candidate: its annihilation cross sections are fixed by EW theory but it can have an arbitrary mass, having relaxed the relic abundance requirement. Namely, for masses smaller (larger) than 9.4 TeV the candidate will be thermally underproduced (overproduced) in cosmology so that alternative mechanisms have to be devised. 3 All our plots apply to this case.
A generic DM candidate which annihilates into gauge boson channels, such as W + W − , ZZ, Zγ and γγ: most of the bounds that we show can apply to this case.
This paper is organized as follows. In section 2 we review the computation of the Sommerfeld enhancement and in section 3 that of the relic abundance. In section 4 we discuss the constraints from galactic diffuse gamma rays (section 4.1), dwarf galaxies (section 4.2) and line searches (section 4.3). In section 5 we draw a summary and conclude.
Detailed phenomenological analyses of the ID signals of pure SU(2) L multiplets as DM candidates have been performed in the past, mostly focussed on the case of the triplet, i.e. the pure Wino (recent examples include [13][14][15][16][17]). On the other hand, the case of the 5plet has attracted significant attention in many different contexts recently [3,6,16,[18][19][20][21][22][23][24], further motivating the analysis that we are embarking on. Other recent works with points of contact to ours include [25,26].

Sommerfeld corrections
At non-relativistic velocities, when the mass of the EW multiplet is much larger than the gauge boson masses, certain radiative corrections to the annihilation cross-section can become large, and dominate the lowest order result. In this regime, the gauge bosons mediate longrange attractive (or repulsive) interactions and the wave functions of the particles involved in the annihilation process become distorted with respect to the plane-waves, enhancing (reducing) the annihilation cross-section. These effects, the so-called Sommerfeld corrections, can be described as a series of ladder diagrams, where the gauge bosons are exchanged before the annihilation process. In practice, the usual pertubation expansion fails, and this class of ladder diagrams has to be resummed. The Sommerfeld corrections have received a lot of attention during recent years in the context of DM theories, see e.g. [27][28][29][30][31][32][33][34][35][36][37][38][39]. In the following we review how to compute these effects.

Sommerfeld factors
The basic step is the factorization of the short range contribution from the long-range effects, which produce the Sommerfeld corrections. The latter ones can be accounted for by computing the wave-function of the system involved in the annihilation, in non-relativistic quantum mechanics.
We classify the initial χχ state with the conserved quantum numbers Q, the electric charge, and S, the total spin. We do so because the potential V relevant for the computation depends both on Q and on S, so that there are in principle different Sommerfeld effects for different values of Q and S. For definiteness, as well as to make the discussion easier to follow, we often refer to the Q = 0, S = 0 initial state: this is the one relevant for indirect detection signals, because the χ 0 χ 0 initial state indeed has Q = 0 and S = 0 (the latter due JCAP10(2015)026 to its Majorana nature and Fermi statistics). All the other initial states are important for the relic abundance computation, because temperatures in the early universe make all the χ components on shell, and coannihilations become relevant (as we will discuss in detail in section 3).
An intuitive picture of the computation can be obtained as follows. We are interested in the cross section for the process a → SM SM, where e.g. a = χ 0 χ 0 . We factorize the computation in a → i and i → SM SM, where i runs over all the states that can be mixed by the interactions, i.e. those with the same Q and S of the initial state a. In the Q = 0, S = 0 case, this amounts to say that i = 1, 2, 3 = χ ++ χ −− , χ + χ − , χ 0 χ 0 . The first factor of the computation has to do with the Sommerfeld effect that we describe in this section, the second factor with the short-distance annihilation process.
We define ψ ia as the two body non relativistic wave-function in the center of mass frame, where the indexes are defined as above. The coupled Schrödinger equation governing ψ ia in the center of mass frame is: where M r i is the reduced mass of the system and E is defined with respect to the χ 0 χ 0 state: Here k 3 = 1/2M χ 0 χ 0 v with v the relative velocity of the χ 0 χ 0 particles. The potential V ij takes into account both the mass splittings and the interactions induced by the gauge bosons. We will give explicit expressions for V ij in section 2.2.
The scattering states at large distances have the asymptotic behaviour: 2) The first term describes the incoming plane-wave wave and δ ia selects the initial state, for instance for the case of χ 0 χ 0 annihilation a = 3. The second part corresponds to the out-going spherical scattering waves. The wave-numbers k i are: with ∆ i the mass splitting, therefore ∆ 1 = 4∆M , ∆ 2 = ∆M and ∆ 3 = 0. The 3D Schrödinger equation can be solved using standard techniques for the scattering problem in quantum mechanics. We move to a spherical coordinate system and we decompose the wave-function in partial-waves l as: where P l (cos θ) are Legendre Polynomials. The 1D Schrödinger equation for the radial function u(r) is: Imposing the asymptotics of eq. (2.2), the radial function should satisfy:

JCAP10(2015)026
where S ia is closely related to the scattering matrix. The Schrödinger equation admits two sets of solutions, which are either regular around the origin, or singular ∼r −l as r → 0. We can select the regular solutions imposing the asymptotics (l = 0 for the s-wave case under study): Summarizing, we should solve eq. (2.5) and impose suitable boundary conditions at the origin and at infinity in order to enforce the asymptotics of eqs. (2.6), (2.7). 4 Finally, the Sommerfeld factors s ia , which should be convolved with the amplitudes describing the short-range interaction (see later), are simply computed for s-wave annihilations taking the value of the wave-function at the origin, see e.g. [27,28]. In terms of the reduced functions u ia , they read: The discussion above shows in simple steps how to solve the scattering problem and extract the Sommerfeld factors. In practice, the differential equation should be solved numerically and some modifications of this method are more appropriate from the point of view of numerical stability. One possibility is to combine the two sets of solutions of the Schrödinger equation, the regular and the singular solutions. The details are discussed for instance in [29,30]. Here we report only the final recipe: We solve the coupled eq. (2.5) for each value of a (therefore 3 times for the case of the Q = 0, S = 0 state) with boundary conditions: where the prime stands for the derivative with respect to r.
This condition originates from imposing the asymptotic behaviour u ia (r → ∞) = R ia e ik i r , where the matrix R ia is determined once eq. (2.5) (of course with the boundary conditions i) and ii)) is solved.
The Sommerfeld factors are s ia = R T ia , where T denotes the transpose.
We work with the adimensional variable y = k 3 r. The numerical origin is set to y min = 10 −7 and we have checked that smaller values do not affect our results. The numerical infinity is chosen in such a way that a stable solution of the differential equation is achieved. A different method to solve numerically coupled Schrödinger equations have been proposed in [41], and it has been adopted to compute the Sommerfeld corrections in [31]. We find this method quite efficient for our purposes, and it nicely agrees with the other technique discussed above.

Sommerfeld-improved annihilation cross-sections
Now we come back to the potential V in eq. (2.5). It can be derived constructing a nonrelativistic field theory for the fermions and integrating out the high energy modes of the spinors and the gauge bosons (see [32]). In addition to the gauge boson interactions, the potential in eq. (2.5) includes also the mass differences between the different components of JCAP10(2015)026 the quintuplet. For the Q = 0, S = 0 states, in the basis χ ++ χ −− , χ + χ − and χ 0 χ 0 , one obtains: . One should pay attention to the difference in normalization between the identical and non-identical two particle states, the former ones get in fact an extra √ 2 factor. The factor √ 2 in the off diagonal 23 entry comes from this mismatch (see e.g. [32]).
Finally, the annihilation cross-section including the Sommerfeld corrections is: with s the Sommerfeld matrix introduced before, and where the index a is not summed over. The factor c a is equal to 2 for the χ 0 χ 0 state (c 3 = 2) and 1 otherwise. For the case of χ 0 χ 0 annihilations we can simply set a = 3. The matrix Γ encodes the short-range annihilation cross-sections. The diagonal entries are simply the σv for the a state into the final state under consideration. 6 Non diagonal entries instead correspond to the quantum interference among amplitudes involving different initial states. They are imaginary parts of the two-point function of a → b. The formalism we have described so far can be applied to compute the annihilation cross-sections of all the components of the quintuplet, in addition to the χ 0 χ 0 initial state that we have mentioned so far. All these processes (annihilations and co-annihilations) are relevant for the calculation of the thermal relic abundance. The computation can be organised as follow. We label the system with the conserved quantum numbers Q and S and we adopt the basis: The potentials, in agreement with [33], are: Using them we can compute the Sommerfeld factors s ia as described in the previous section. Then, to obtain all the relevant annihilation cross-sections σv from eq. (2.10), one JCAP10(2015)026 needs the following tree-level s-wave annihilation cross sections [33]: where in (2.13) the first matrix refers to the cross section into W W , and the second one to the cross sections into all the other gauge bosons. The cross sections into ZZ, γγ and Zγ can be obtained from that second matrix via multiplication by c 4 w , s 4 w and 1 − s 4 w − c 4 w respectively. We conclude this section showing, in figure 1, the χ 0 χ 0 annihilation cross-sections into several final states. For the galactic environment, we fix the typical DM relative velocity to the value v = 10 −3 c. The series of peaks and dips in the figure are indeed the manifestation of the Sommerfeld effect. 7 These results update those obtained in [33]. We find a good overall agreement, although the precise details of the peaks at large mass (which are important for our purposes) are somewhat modified.
We first notice that the Sommerfeld enhancement can be very big, enhancing the value of the cross sections by several orders of magnitude. Also, while the tree-level short-distance computation would yield σv = 0 for the process χ 0 χ 0 → γγ, one observes here a rather large value of this cross section. The presence of such a process is due to the fact that the Sommerfeld effect mixes the χ 0 χ 0 initial state with χ + χ − and χ ++ χ −− , which couple to JCAP10(2015)026 photons. The fact that the related cross section is rather large is due to the initial state consisting in two doubly charged particles, which results in a factor of 16 enhancement in the cross section, with respect to the analogous process with χ + χ − instead of χ ++ χ −− . We also show the behaviour of the cross sections as a function of the relative velocity, for two given fixed values of the M DM . This makes evident that, below a certain threshold which depends on the specific mass, the cross sections stop growing and reach a constant value.

Relic abundance
In this section we review the computation of the thermal relic abundance of the lightest component of the 5plet (the Dark Matter particle) as a function of its mass M DM . By demanding that it makes all of the measured DM in the Universe (Ω DM h 2 = 0.1188 ± 0.0010, as determined by the latest Planck results [42]) we can univocally determine its mass.
In principle, the evolution of the number density n α of each α-th component of the 5plet χ α ≡ (χ 0 , χ + , χ − , χ ++ , χ −− ), with internal degrees of freedom g α ≡ 2 and mass M α ≡ M DM + ∆M α , has to be computed by solving a system of five coupled Boltzmann equations. The densities are affected, in addition to the process of χ 0 χ 0 annihilation into any SM state, by all the possible χ α χ β co-annihilations. Indeed, since the states are almost degenerate in mass, the χ α χ β co-annihilation cross sections play a major role in setting today's relic abundance. However, in practice, it turns out that it is sufficient to follow the evolution of the total number density n = α n α , i.e. it is sufficient to just solve one simple Boltzmann equation in terms of the total thermally averaged annihilation cross section which includes all the possible χ α χ β co-annihilation channels. This is because n α is much smaller than the SM thermal bath number density n SM 8 and, at the same time, the scattering cross sections σ α SM → β SM are of the same order of the annihilation cross sections σ αβ → SM SM . Hence the scattering rates off SM particles are much faster than the χ α χ β annihilation rates and the χ α particles are kept in equilibrium with the thermal bath. As a consequence, the only relevant processes are those that modify the total density n. This greatly simplifies the problem.
As it is customary, we actually define the total comoving density Y = n/s, where s is the total entropy density of the Universe, and we follow the evolution of Y in terms of the adimensional parameter x = M DM /T , with T the temperature of the thermal bath. In terms of x, the entropy density reads s(x) 2π 2 /45 g * s (x)M 3 DM x −3 and the Hubble rate in the radiation domination era reads Here g * s (x) and g * (x) are the effective relativistic degrees of freedom as given in ref. [43].
We compute the cross sections as described in section 2, i.e. combining the Sommerfeld effect with the s-wave short distance contributions of eqs. (2.13), (2.14), (2.15). This procedure has an uncertainty coming from the neglection of the p-wave contributions to eqs. (2.13), (2.14), (2.15) [2], as well as from having neglected the relativistic corrections to the Schrödinger equation (2.1). We account for these effects by assigning a rough v 2 /c 2 uncertainty to the thermal relic mass, which is of the order of 5% at freeze-out. We have accounted for the dominant thermal corrections (as in [33]), which however have a very small impact.

JCAP10(2015)026
We can now write the Boltzmann equation for Y as a function of x. Following refs. [43,44], it explicitly writes is the equilibrium comoving density, K 2 is the BesselK function of the second order and g tot (x) = α g α (1 + ξ α ) 3/2 e −x ξα , with ξ α = ∆M α /M DM , are the total effective degrees of freedom of the 5plet at a given temperature corresponding to x. Notice that, in the limit of degenerate states (ξ α ≡ 0 for all i), g tot (x) = α g α = 10. The initial condition of Y in eq. (3.1) reads as usual Y (x 0 ) = Y eq (x 0 ), where we fix x 0 = 4 in order to follow the whole subsequent evolution.
In eq. (3.1), σ eff v is the total thermally averaged effective annihilation cross section in the cosmic comoving frame, which is the most important quantity to determine the today's relic abundance of the MDM 5plet. Under the reasonable assumption of thermal equilibrium of the χ α state with the SM thermal bath, as discussed above, it writes [45] where K 1 is the BesselK function of the first order and the σ αβ are all the possible annihilation cross sections into any SM state (see section 2). 9 Since, for velocity dependent cross sections (e.g. the ones which include the Sommerfeld corrections) and large x, the integral in eq. (3.3) is numerically hard to solve, for x > 100 we introduce the variable = 1/x and we expand the function K 1 √ s/( M DM ) /K 2 2 (1/ ) for small . We have checked that this approximation is very accurate by comparing the exact integral in eq. (3.3) with the approximated one, for velocity independent cross sections. Considering all the possible CP combinations of the initial states, eq. (3.2) in the limit of degenerate masses explicitly writes The right panel of figure 2 shows the total thermally averaged annihilation cross section (with (solid lines) and without Sommerfeld (dashed lines)) for two indicative values of the DM mass (1 TeV (gray), 9.4 TeV (black)). As one can see, for x → 0 (relativistic regime), σ eff v is not approaching a constant value because the thermal averages taken in the cosmic comoving frame and in the center-of-mass frame do not coincide in the relativistic regime [43]. As an aside the bump of the Sommerfeld enhanced cross sections, before the plateau at large x, is basically the convolution of the many peculiar resonant peaks of the Sommerfeld. On the other hand, the little increase at x M DM /(166 MeV) of the non-enhanced cross sections is due to the decoupling of the χ + , χ − , χ ++ , χ −− states.
Having at our disposal the total thermally averaged annihilation cross section, we can finally integrate numerically eq. (3.1) in order to determine the asymptotical value of the JCAP10(2015)026 1.05 × 10 −5 GeV/cm 3 are the today's entropy and critical energy densities of the Universe respectively [42]. The left panel of figure 2 shows Ω DM h 2 as a function of the DM mass. The dashed and solid lines refer respectively to the computations with and without the Sommerfeld effect, while the horizontal strips individuates the measured DM in the Universe at 1σ (2σ) CL by Planck 2015. The solid line crosses the 2σ CL band when M DM = 9.4 ± 0.094 TeV. We also show, as a vertical yellow band, the O(5%) correction due to theory uncertainty on the determination of the cross sections and, in turn, on the value of M DM . In this case we get then M DM = 9.4 ± 0.47 TeV. In the rest of the paper, since the experimental error in the determination of M DM is smaller than the theoretical one, we will always consider the latter.

Gamma ray constraints
In this section we compare the different components of the gamma ray spectrum, as predicted by MDM, to the experimental data. We will consider in turn the measurement of the galactic diffuse emission by Fermi, the results on dwarf galaxies observations by, again, Fermi, Hess and Magic and the gamma ray line searches performed by the same experiments.
First, however, let us shortly remind of the characteristics of the DM distribution in the targets that we consider. For the Milky Way the range of possible profiles, as schematized e.g. in [46], spans from the Burkert and Isothermal profiles, featuring a constant density core in the inner kiloparsecs of the Galaxy, to the peaked Navarro-Frenk-White (NFW) one, which is formally divergent as r approaches the GC. The Einasto profiles are not divergent but still peaked at the GC. All these profiles need to be normalized by assuming a value for a scale radius and a scale density: we follow the procedure described in [46] which amounts to fix the density at the location of the solar system to ρ = 0.3 GeV/cm 3 and to impose that the total DM content of the Galaxy (within 60 kpc) agrees with the recent estimates. The Einasto profiles depend on an additional parameter (α) which controls their steepness in the GC region: as in [46], we will consider in some cases the standard value α = 0.17 and also a steeper version featuring α = 0.11 (dubbed EinastoB). For the case of the dwarf galaxies, the uncertainty on the DM distribution reflects essentially in the uncertainty on the so-calledJ factor. We will discuss this in some detail in section 4.2.

Galactic gamma ray diffuse emission measurement by Fermi
In this section we derive the bounds from the whole sky measurements provided by Fermi. We will derive two kinds of constraints: Conservative constraints: we suppose a vanishing gamma ray background from astrophysics and we just impose that the DM signal does not exceed the measured flux.
Constraints obtained with a modeling of the background: this of course reduces the room for a DM signal and therefore leads to bounds that are more stringent. They are however also less robust, as they depend on the reliability of such modeling.

Analysis setup
Choice of the regions of interest and data extraction. We divide the whole galactic sky observed by Fermi in 35 non-overlapping 'regions of interest', as depicted in figure 3, masking out the 2 • around the galactic plane (but retaining a 2 • × 2 • region around the GC). The regions are designed to be smaller near the GC and wider at high latitude and longitude.

JCAP10(2015)026
Some of them correspond to areas already used by previous analyses (e.g. our RoIs number 12 and 24 coincide with the 'mid-latitude strip' used by the Fermi collaboration itself in [47]), which allows us to make quantitative comparisons with previous results. For reasons that will be clear later, we distinguish between the Inner Galaxy (|b| < 15 • , | | < 80 • , RoI's from 12 to 24) and the Outer Galaxy (the corresponding complement).
We extract Fermi data using the Fermi Science Tools v9r32p5. We use 5 years of data within the event class CLEAN. We perform the following selection cuts: (DATA QUAL==1) && (LAT CONFIG==1) && (ABS(ROCK ANGLE)<52) && (IN SAA!=T). Events with zenith angles larger than 100 • are excluded. The exposure is computed using the Fermi-LAT response function P7REP CLEAN V15. The data are binned in 31 energy bins equally spaced in log scale between 300 MeV and 500 GeV. Then, in order to increase the statistics at high energies, we have grouped the last four bins into two wider bins. The counts and the exposure maps have been produced using the HEALPix pixelization scheme [48], with a resolution n side = 256, corresponding to a pixel size of ∼ 0.23 • . The error bars on the differential photon flux are obtained summing in quadrature statistical and systematic (from [49]) uncertainties.
Background modeling. Anything other than the gamma rays produced by DM annihilation is astrophysical background for our purpose. We need therefore to have a reliable modeling of it, in order to be able to gauge the DM contribution still allowed by data. However, designing such a modeling is a challenging task by itself. We discuss in the following the procedure that we follow and how we assess its reliability.
We consider several background components: (I) a template for the diffuse galactic emission produced by charged CR, via (Ia) interactions on the interstellar gas and via (Ib) the Inverse Compton process; (II) a template for point sources, as derived in [50,51]; (III) a template for the so-called 'Fermi bubbles', as provided in [52]; (IV) the isotropic gamma ray background as measured in [53], including their estimate for the irreducible charged CR contamination.
While the components (II), (III) and (IV) are rather straightforwardly implemented, (I) deserves a dedicated discussion. The Fermi collaboration provides the template for (I) (in the supplementary material of [53]), which we adopt. They provide three different versions (Model A, B and C, with A being the benchmark one, B adding extra CR sources and C considering variations of the charged CR diffusion coefficient): we adopt Model A for definiteness and we have checked that our procedure is only very marginally affected if choosing another one. The collaboration then corrects the template with renormalizing coefficients, energy bin per energy bin, in order to best-adapt it to the data in the wide region |b| > 20 • . We follow therefore the same procedure, by inferring the renormalization coefficients from [53]. This is however subject to two caveats. First, we have found that the inferred coefficients are good up to a ∼ 20% error, since in [53] some of the subdominant foregrounds accounted for in the analysis have not been shown explicitly [54]. Second, we are applying the renormalizing coefficients deduced from the wide |b| > 20 • region to our RoI's which are in general different (smaller and covering also lower latitudes). However, a posteriori we will find that the agreement with the data is good within the limits that we impose (see below). Following the procedure of the collaboration (see section 3 of [53]), we infer a renormalizing coefficient for

JCAP10(2015)026
each energy bin up to 13 GeV and a single coefficient for all data points above that threshold. This is done separately for the (Ia) and (Ib) components. The procedure discussed above is as accurate as possible for the scope of our analysis. We test it against the data in the different RoI's and we indeed find that it provides a good description of the background for photon energies within the following windows: for the RoI's in the 'inner galaxy' 1.5 GeV E γ 500 GeV; for the 'outer galaxy' 1.5 GeV E γ 100 GeV. Indeed, the reliability of our template for point sources worsens rapidly beyond 100 GeV in the outer regions. Also, in these 'low contrast' areas a more rigorous treatment of the contamination from charged cosmic rays would probably be needed.
To summarize: we borrow the background modeling used by the Fermi collaboration for their galactic diffuse study, but we adapt it to our spatial RoI's and we limit our analysis to the energy ranges in which we find that it provides a reasonable description. In doing so we will derive constraints that are necessarily less stringent, i.e. more conservative than they could in principle be (if we had kept the whole dataset and if we had re-optimized the background to all the regions).
Dark Matter contribution. 5plet Dark Matter annihilations contribute to the galactic halo continuum gamma ray spectrum measured by Fermi mainly via the DM DM → W + W − channel but also via the channels DM DM → ZZ and DM DM → Zγ. To a much lesser extent, also the DM DM → γγ contributes: EW radiation can degrade the monochromatic final state photons and generate a continuum lower-energy flux. The latter three channels give a contribution which is as relevant as the W + W − one, thanks to the Sommerfeld enhancement of their cross sections, as discussed in section 2. We therefore sum the gamma ray yields from all these processes, with the relative cross sections precisely computed in section 2.
In addition, besides the gamma rays promptly emitted in the annihilation, we also include in the computation of the full spectrum the Inverse Compton Scattering secondary gamma rays (which originate when electrons and positrons from the DM annihilation scatter against the ambient galactic light). For this purpose we use the tools provided in [46,55], to which we refer for all details.
It is worth noticing that, while we are here interested in considering the specific case of the 5plet, the bounds that we obtain can be safely extended to a more general class of models, in which multi-TeV DM generically annihilates into gauge bosons. This is because the shapes of the continuum gamma rays from any V V channel (with V a gauge boson) are essentially indistinguishable from one another, as explicitly shown e.g. in figure 3 of [46], apart from the features at the endpoint of the spectrum which are present for the Zγ and γγ channels (from the monochromatic γ) and for the W + W − one (from the final state radiation of a hard photon). Hence, provided that one considers DM heavy enough that these features fall beyond the sensitivity range of Fermi (i.e. for M DM 500 GeV), the constraints that we derive apply to the spectral shape corresponding to any gauge boson channel. One just needs to rescales the bounds with the different γ-ray multiplicities from the different channels.

Background-free conservative constraints
In figure 4 (upper panels) we present the constraints that we obtain by imposing that the DM signal does not exceed by more than 3σ any of the data points in any given RoI. We consider two representative DM profiles (NFW and Burkert) and we span a large range of DM masses, from 100 GeV to 30 TeV.  We see that there is no single RoI which dominates over all the range of masses. However, as expected, regions close to the GC (e.g. RoI 21, 15 and 18) impose the strongest constraints for the peaked NFW profile while wide regions in the outer galaxy (e.g. 3,4,27) are the most relevant for the cored Burkert profile. To better illustrate this point, in figure 5 we choose a definite value of the DM mass (9.4 TeV, as obtained in section 3) and we report on the galaxy map the strength of the different constraints.   On the other hand, considering the envelop of the most constraining curves on the whole range of M DM , we see that the NFW hypothesis provides stronger bounds than those from Burkert, but the difference never goes beyond a factor of a few.
In figure 4 we superimpose the total cross section in the channels mentioned above. 10 The mass intervals for which the line enters in the shaded regions are excluded. We see that, while significant constraints can be imposed especially below ∼ 2 TeV, the model is still allowed for large intervals of masses. For the specific case of the MDM 5plet (M DM 9.4 TeV), the bound lies almost 2 orders of magnitude above the predicted cross section, thanks to the fact that the Sommerfeld enhancement displays a trough at that value in mass.

Constraints including background
In figure 4 (lower panels) we present instead the constraints that we obtain by including an astrophysical background. Here the bounds are derived by looking for the best fit (background + DM) configuration and then requiring that the addition of more DM does not worsen the best fit χ 2 by more than ∆χ 2 = 9. The situation changes in two respects. First, of course the bounds are much stronger, as less room is left for Dark Matter. Second, the relative importance of the different RoI's in setting the bounds changes, as a consequence of the fact that our background modeling may describe accurately or not the measured flux in 10 We denote the bounds as being on the annihilation cross section DM DM → V V , with V a gauge boson W , Z or γ. Namely, here we have computed the bounds considering the case of a pure W + W − annihilation, but they can be recast for different channels. We refer to the discussion in the previous paragraph for more details. We see that now the measurements rule out essentially all the region below M DM 7 TeV. For larger masses, small islands are excluded up to about 25 TeV. The MDM 5plet is again spared both for the NFW and the Burkert cases. The constraints from NFW and from Burkert are now even more similar.

Dwarf galaxies observations of gamma ray continuum
Dwarf satellite galaxies are believed to be some of the cleanest possible laboratories to search for DM in gamma rays, thanks to their presumed high DM content and relatively reduced stellar emission foreground. On the other hand, the scarcity of stellar tracers makes it difficult to precisely reconstruct how much DM they actually host and how it is distributed, leading to large uncertainties. For the case at hand, a related uncertainty has to do with the Sommerfeld enhancement: in order to compute it precisely, one would of course need to know the DM velocity dispersion, which is in principle different in each galaxy and in different radial positions within each galaxy. As customary, we assume that the DM velocity is the same as that of the stellar tracers (which is plausible if the systems have reached relaxation) and we adopt a common value of 10 km/s = 3 × 10 −5 c, in line with the measurements that typically span 3 to 15 km/s (see e.g. [56] for a recent compilation).
The uncertainties on the DM content and distribution for each galaxy are the subject of a long ongoing debate in the literature (see e.g. [57][58][59][60][61][62]). They are commonly expressed as uncertainties in the determination of the so-called averageJ factor, where ∆Ω is angular area subtended by the dwarf spheroidal galaxy or by the instrument and θ is the angle between the axis connecting the Earth and to the dwarf galaxy and the line of sight. The Fermi collaboration quotes a rather small uncertainty at the level of 10% to 40% at most [63] for the stacked sample of galaxies. For the individual dwarf galaxies (in the 'classical' class), a recent study [64] finds values from 17% to 124%. The very recent dedicated study in [56] typically finds larger uncertainties, which can however differ case by case. Finally, in case of somewhat more exotic scenarios in which dwarf galaxies host intermediate mass black holes, the flux can be modified by a factor from a few up to 10 6 [65].
In view of this situation we are prompted to consider, alongside our constraints, alternative ones that are somewhat relaxed (as we will detail below). Future observational work on these systems has clearly the potential of reducing these uncertainties significantly.
We base our analysis mainly on the most recent observation of 15 dwarf galaxies by Fermi [63]. We also consider the constraints imposed by Hess [66], which has observed a subset of four of these dwarf galaxies, plus Sagittarius (not included in the Fermi analysis), and Magic [67], which has intensively observed Segue1. The Hess constraints complement those of Fermi in the window M DM = 10 → 20 TeV while the Magic ones remain subdominant to those of Fermi, since they do not extend beyond 10 TeV.
In figure 6 we report the bounds on the annihilation in the W + W − annihilation channel obtained by the experimental collaborations and we compare them with the predicted 5plet DM cross section into W + W − + ZZ + Zγ/2, which we denote as V V like before. Indeed, the spectral shapes of continuum γ-rays from these channels are very similar (as discussed above) and we neglect here the largely subdominant continuum γ-rays from the γγ channel.  The fact that the W + W − channel features a sharp peak at E ∼ M DM , from the emission of a hard final state photon, is not expected to modify these results significantly, also in light of the much larger uncertainties connected to the DM contents determination.
We find that the constraints from Fermi, taken at face value, exclude all the range M DM 7 TeV. We also consider bounds relaxed by an order of magnitude (labelled 'Fermi -'), intended as an average conservative assessment of the uncertainties related tō J-factors discussed above. This slightly reduces the excluded interval to M DM 6.3 TeV. The constraints from Hess rule out very small intervals between 10 and 20 TeV. The bounds that do not include the most stringent dwarf (Sagittarius) are slightly less constraining.
We stress that all these numbers are strongly dependent on the detailed shape and position of the Sommerfeld peaks, which in turn suffer from a 'theory uncertainty' that we have quantified before to be of 5% to 10%. However, the global results are robust. The mass for the MDM 5plet is again not affected by these bounds, unless very aggressive assumptions on theJ-factor determinations, such as e.g. to lower significantly the Fermi bound, are made.
To conclude this section, we point out that significant progress could be made by current and future experiments. Improving the Fermi, Magic or Hess bounds by a factor of just a few at the largest masses would allow to probe almost the entire parameter space of the model. In this respect, choosing one of the dwarf galaxies with the most promisingJ-factors (such as Coma or Ursa Major II, according to [56]) can perhaps allow to reach such a goal.

Gamma ray lines searches
As already mentioned above, γ-ray lines (or sharp spectral features which are degenerate with lines for experimental purposes) arise at the endpoint of the spectrum from DM annihilation. Searching for these features has been regarded since a long time as a very promising strategy, and they have been often referred to as the proverbial 'smoking gun' for Dark Matter [68][69][70][71][72][73][74][75][76][77][78][79][80][81].

JCAP10(2015)026
The Hess telescope has imposed line constraints in [82], the relevance of which, for pure wino DM, was initially recognized by [13,14]. The region observed by Hess, the Central Galactic Halo (CGH), is a promising one due to its relative proximity and large predicted DM concentration. The search region is defined as a circle of 1 • radius centered on the GC, where the Galactic plane is excluded, by requiring |b| > 0.3 • . 11 Magic has also published γ-ray line constraints [67] from the observation of the Segue 1 dwarf galaxy. Finally, we will also quote the bounds from the Fermi collaboration, presented in [83], which are relevant at small DM masses.
We consider the 5plet DM annihilation channels γγ and γZ. Some recent works have advocated the need to go to higher orders in the cross sections for heavy EW multiplets, and have provided refined fixed order [36] and resummed [85,86] calculations, in particular for the exclusive cross section into γγ. The effect of these resummations is a reduction of the annihilation cross section by a factor of ∼ 3 for EW triplets of mass of a few TeV. Recently, however, the issue has been revisited in [30], where it has been pointed out that a less exclusive cross section would be more appropriate for this kind of channels, given the limited resolution (hundreds of GeV) of the current experiments for high energy γ rays. Indeed the relevant cross section is the one for γγ + X, where X is anything soft enough not to affect the line, once the experimental resolution is taken into account. 12 Ref. [30] presents a result of this procedure for a pure Wino Dark Matter candidate and shows that the resummed cross section is barely distinguishable from the lowest order one, in the interesting mass range. In light of the above discussion, as well as for simplicity, we stick to the lowest order cross section in γγ + γZ/2, to which we add the effect of Sommerfeld enhancement.
In figure 7 we show the constraints superimposed to the DM annihilation cross section as computed in section 2. Notice that Hess quotes in [82] a bound on an Einasto profile with α = 0.17, scale radius r s = 20 kpc and scale density ρ s = 2.8 × 10 6 M /kpc 3 , 13 which we have to rescale to our standard profiles. We see that the Fermi bounds (which we report for two different profiles and regions, marginally different) rule out the low mass portion. The Hess constraints are very relevant in the 500 GeV → 20 TeV window: if the profile is assumed to be peaked (such as NFW or Einasto) the whole region is essentially ruled out. If the profile is cored (such as Burkert), some intervals are reopened.
The MDM 5plet falls in one of these intervals: if the profile is peaked, it is ruled out; if the profile is cored, it is again spared.
The constraints from Segue1 by Magic are equally relevant and have the power to rule out essentially the entire window up to M DM = 7 TeV. However, as we commented in section 4.2, they are subject to a large indetermination. In particular, the Segue1 dwarf spheroidal is found in [56] to have a largeJ-factor uncertainty. If, for illustration, we rescale 11 The Hess search is performed adding to the background a Gaussian-shaped line, in different bins of energy, and thus is effective for flat DM profiles. This is opposite to the Hess searches for a γ ray signal in the inner region of the galaxy [84], whose constraints do not apply for profiles with a large core, since there the background is estimated with the On/Off procedure from a control region close to the center. 12 We are consistently neglecting the radiation of a very hard γ from W bosons. Such an emission is Sudakov enhanced, and it gives a feature at Eγ MDM which resembles a line -except for a shift of order mW , which serves as a cut-off of the soft divergence. Given that the experimental resolution is larger than mW , this contribution is observationally indistinguishable from γγ and γZ lines. Moreover, to be consistent in including this effect, one would need to compute the same processes, at the same order in the EW couplings, including splittings like γ → W W . For inclusive enough observables (as it is the case, given the experimental resolutions), these effects compensate with the extra photons radiated by the W bosons. 13 It is important to realize that these parameters imply a total DM content of the Milky Way and a local DM density much higher than the ones that we adopt (see the discussion at the beginning of section 4).  the bound by the '1σ' of the value determined in [56] for theJ-factor of Segue1, the constraints relax significantly (see the line labelled 'Magic -' in figure 7).
We conclude this section by briefly mentioning the prospects for improvements. In figure 7 (left) we report the sensitivity reach for the Cta observatory as quoted in [86] (on the basis of [87]), rescaled to the case of the Burkert profile. We see that significant additional portions of the parameter space would be covered, but still large gaps would remain at high masses, due to the peculiar peak structure of the Sommerfeld enhanced cross section. In particular, the MDM 5plet would still not be tested.
Concerning the line searches in dwarf spheroidals, we stress again that improving the existing bounds by a relatively small factor, e.g. focussing on target galaxies with potentially highJ-factors, would allow to complete the coverage in mass, including the value of the MDM 5plet.

Conclusions
Motivated by the Minimal Dark Matter model, and more generally by the model-class of multi-TeV Dark Matter with EW interactions (multi-TeV WIMPs), we have explored the constraints that come from several gamma ray probes. A crucial ingredient for these kinds of models, as recognized since quite some time, is the Sommerfeld enhancement which arises from the exchange of EW bosons among the heavy DM particles: it modifies significantly the annihilation cross section, both at DM thermal freeze out and in the current universe, and gives rise to a peculiar structure in peaks, that we recomputed in detail (section 2).
In figure 8 we compile the bounds in a summary chart: we shade away the intervals in mass which are excluded by each one of the gamma ray probes that we have considered and we contour with a dashed line the regions explorable with near future improvements.
The MDM 5plet, which has a mass of 9.4 TeV as determined by a careful computation of its relic abundance (section 3), is severely tested by the Galactic Center line constraints from  Hess (section 4.3): if the DM galactic profile is peaked like NFW, the model is ruled out; if the profile is cored like Burkert (or Isothermal), the model is still allowed. In this latter case, future Cta line searches may enlarge the explored window but still are not expected to test that precise value in mass. In addition, future Cta observations of the inner Galaxy are expected to have the power to explore the whole mass range 100 GeV M DM 30 TeV [88][89][90] (applicable however only to non-flat profiles).
More generally, our results show that multi-TeV WIMPs can be significantly probed with Indirect Detection via gamma rays. Also, independently on the specific case of the 5plet with its peculiar roller-coaster cross section, the constraints in figure 4 (and figure 6) actually apply to any DM with M DM 500 GeV that annihilates into W + W − , ZZ, Zγ or γγ, as we have discussed above.
Other channels of Indirect Detection for multi-TeV WIMP DM are possible, but are expected to have a somewhat shorter reach than the one we have considered here. The case of antiprotons is particularly interesting and has been explored intensively recently. Using the results in [91][92][93][94] and, in particular, the most recent assessment based on Ams-02 data [95], we can obtain constraints on the W + W − cross section of section 2. We show the outcome in the summary chart as a shaded orange area: despite their relevance for smaller masses, these constraints do not change the global picture sketched by the gamma ray ones.
In conclusion: Minimal Dark Matter, and more generally multi-TeV WIMP DM candidates, are arguably even more motivated than before, in the current context of absence of New Physics from the LHC. Gamma rays are a powerful probe for this class of models. The MDM 5plet candidate is ruled out or still allowed depending on the DM profile at the Galactic Center. Significant future progress is possible and may notably come from the observation of dwarf spheroidal galaxies.

JCAP10(2015)026
Note added. During the preparation of this work, we became aware of another group independently investigating gamma ray signals of Minimal Dark Matter scenarios [97].