The 750 GeV Diphoton excess, Dark Matter and Constraints from the IceCube experiment

Recent LHC data show hints of a new resonance in the diphoton distribution at an invariant mass of 750 GeV. Interestingly, this new particle might be both CP odd and play the role of a portal into the dark matter sector. Under these assumptions and motivated by the fact that the requirement of $SU(2)_L$ invariance automatically implies the coupling of this alleged new resonance to $ZZ$ and $Z\gamma$, we investigate the current and future constraints coming from the indirect searches performed through the neutrino telescope IceCube. We show that these constraints can be stronger than the ones from direct detection experiments if the dark matter mass is larger than a few hundred GeV. Furthermore, in the scenario in which the dark matter is a scalar particle, the IceCube data limit the cross section between the DM and the proton to values close to the predicted ones for natural values of the parameters.


Introduction
The ATLAS and CMS collaborations have recently announced [1,2] (for a recent update see Refs. [3,4]) the observation of an excess of events in the search for two photons in the final state. The shape of the excess suggests the detection of a resonance at an invariant mass of approximately 750 GeV, and preliminary analyses suggest a rather large value for the width of the resonance, around 45 GeV (equivalent to 6% of its mass).
The statistical significance of this observation is still far from conclusive. The ATLAS collaboration, with 3.2 fb −1 of data, claims a statistical significance of 3.9σ (or 2.3σ by taking into account the look-elsewhere effect) with an excess of about 14 events, corresponding to a cross section of about 10 ± 3 fb. The CMS collaboration partially supports this observation, with a weaker statistical significance (local significance of 2.9σ) due to a smaller integrated luminosity of 2.6 fb −1 , and with an estimated cross section of 6 ± 3 fb.
If the excess were confirmed by the data collected in the continuation of the Run-2 of the LHC experiments, this discovery would represent a historical cornerstone in the investigation of fundamental interactions beyond the Standard Model (SM).
It is tempting, although speculative, to try to relate this hypothetical new particle to the fact that about 30% of the energy density of the universe seems to be in the form of Dark Matter (DM) particles. If the indication on the total width Γ P of the resonance persists, then the possibility of a coupling of this new resonance to DM with a sizeable branching ratio (BR) would gain stronger support, because it would be easier to obtain such a large width in a simple, weakly coupled model. In that case, a natural scenario would be that P acts as a portal between the SM and DM [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24]. Given that the mass scale of the resonance is above the electroweak scale, it is reasonable to write an effective theory describing the interaction between P and the photon in an SU (2) L invariant way. In this case, an interaction of P with photons automatically implies the existence of the vertices P ZZ and P Zγ, with fixed couplings (see Ref. [25] for a Run-2 ATLAS analysis analysis of Zγ final states). This very general statement has important implications for phenomenology, first of all at LHC but also in DM searches, if indeed P acts as a portal into the DM sector.
The observed decay channel into two photons already constrains the spin of the resonance (which we denote by P ) by virtue of the Landau-Yang theorem [26,27], forcing its spin to be either 0 or 2. 1 In this work, we focus on the possibility that P is a particle of spin zero, which interacts with a DM particle charged under a Z 2 symmetry that prevents it from decaying and makes it a viable DM candidate. We consider the two cases of a Dirac fermion χ and of a complex scalar φ. The discussion of the DM phenomenology of this model is strongly affected by the assumption on the CP properties of P . At the energy scales involved at the LHC, the production cross section of P is basically independent of whether P is a scalar or a pseudoscalar particle, whereas at low energies the two options bring to different non-relativistic effective operators.
In the scalar case, the effective operator for the interaction between DM and nuclei is spin independent (SI), yielding strong bounds from direct detection (DD) experiments. If however P is a pseudoscalar particle, then the low energy effective interaction between DM and a nucleus is spin dependent (SD). In this case, the exclusion reach of DD is much weaker and has to be complemented with the constraints coming from indirect detection (ID), because the DM annihilation at low velocities occurs through an s-wave for the pseudoscalar case (while in the scalar case it occurs through a p-wave process).
In this paper we assume that the particle P is a pseudoscalar 2 and assess the constraints coming from the IceCube (IC) experiment [33], a neutrino telescope that can be used to study the DM annihilations occurring in the Sun [34,35]. DM particles can get captured in the gravitational well of the Sun if they scatter with atomic nuclei inside it and they lose some energy. The accumulation of DM particles is partly compensated by their annihilation, until the DM density in the Sun reaches an equilibrium level. If this is already achieved today, then one can directly relate the annihilation rate to the capture rate in the Sun and constrain the interactions of the DM with the 750 GeV resonance as well as the interactions of the resonance with the SM particles.
The following features make particularly interesting the study of the 750 GeV resonance 1 Consequences of the Landau-Yang theorem can be evaded in scenarios in which a vectorial resonance decays via a cascade into a final state of three photons, two of which are too collimated to be discriminated in the detector (see e.g. Ref. [28]). For a discussion about the fate of Landau-Yang theorem in non-Abelian gauge theories, see [29]. 2 See Refs. [30][31][32] for a discussion of a pseudoscalar field coupled to B B as a heavy hypercharge axion.
with the IC experiment, if the pseudoscalar P acts like a portal to DM: • The decay of P into ZZ and Zγ, granted by the assumption of SU (2) L invariance, ensures the presence of at least two annihilation channels for DM that yield energetic neutrinos which can be observed by IC.
• Given the assumption of equilibrium inside the Sun, neutrino fluxes at Earth depend only on the BR's of DM annihilating into SM primary products. In the expression of the BR's, interesting simplifications occur and the bound from IC turns out to be independent of the total width Γ P , of the coupling of P to the DM particle, and of the mass m P of the resonance. Furthermore, the expressions of BR's can be rewritten as functions of the partial decay widths of P , which are the quantities that can be directly measured via a resonant production at LHC.
• IC constraints reach the highest exclusion power for DM masses of order 10 −1 to 1 TeV [34]. Indeed, for DM masses higher than a TeV, the DM number density in the Sun is reduced and this, together with the fact that primary neutrinos from DM annihilations hardly escape without interacting within the Sun and losing much of their energy, deteriorates the limit. At lower DM masses, the angular resolution of IC is poorer because of the lower number of Cherenkov photons produced by muons at these energies. Interestingly, this mass window coincides with the order of magnitude m P . Thus, if the DM particle has a mass close to m P , IC might give the strongest bounds on DM.
This paper is structured as follows. Sec. 2 describes the different models we consider, and the benchmarks we choose for the couplings of P to SM particles which are consistent with LHC observations. Sec. 3 illustrates the physics that links the observations of IceCube and DM annihilations, and the importance of electroweak (EW) corrections at the energy scales of interest to IC. Sec. 4 contains the results for the constraints from IceCube, DD and γ-rays observations (Fermi-LAT, HESS) on the benchmark models we consider. Finally, in Sec. 5 we summarise our results.
2 The 750 GeV resonance as a portal to DM: the models The observation of P in final states with two photons forces P to be coupled at least to the photon field strength. According to our assumption, P is a pseudoscalar and the effective vertex for its interactions with photons must be of the form P F µν F µν in order not to introduce a source of CP violation, where F µν = 1 2 ε µνρσ F ρσ is the dual field strength. Since the theory must have a cut-off at least higher than m P = 750 GeV, it is more than reasonable to write down an explicitly SU (2) L -invariant Lagrangian. Given the current lack of excess in W W , ZZ and Zγ searches, the coupling of P to B B is favoured with respect to a coupling to W i W i (where i = 1, 2, 3 is the SU (2) L index). Thus, for minimality, we set to zero the coefficient of the vertex P W i W i . 3 The absence of anomalies in the diphoton searches of Run-1 at LHC strongly favours the hypothesis that P couples also to gluons 4 and/or quarks, because the parton distribution function (PDF) for photons changes only mildly from 8 TeV to 13 TeV, whereas the experimental results of Run 1 and 2 can be consistent if P is produced mainly by gluon fusion or from heavy quarks (see for instance Ref. [5]). We ignore the coupling of P to leptons, since this is irrelevant for diphoton production and has a minor impact for IC and DD. We consider therefore the following Lagrangian density where Λ is the dimensionful scale of the effective theory, a = 1, . . . , 8 is the gluon SU (3) c index, i is the family index, H is the Higgs doublet and H c = iσ 2 H * is its conjugate, Q i L is the quark weak doublet, and c gg , c BB , y qP , y χP are real coefficients.
As for the DM Lagrangian we envisage two possible (and mutually excluding) cases: either the DM is a Dirac fermion χ (the results being analogous in case of a self-conjugate DM particle) or it is a scalar particle φ with CP-conserving Lagrangian where A P has to be thought of as a spurion field which changes sign upon CP (for example it might be the vacuum expectation value of some heavy parity-odd field). Notice that in both cases the DM could be identified with one of those particles in the fermion or scalar multiplets coupled to the resonance and which, upon integrating them out, give rise to the effective interaction between P and photons and gluons. If so, the DM must be of course a singlet under the SM gauge group.
The couplings c P P φφ and c HHφφ of Eq. (2.3) are not currently constrained by the observation of the decay of P into two photons, thus they cannot be linked to the quantities measured by the experimental collaborations. Furthermore, the term P 2 |φ| 2 does not alter quantitatively the branching ratios of the annihilation of φφ to SM particles relevant for the bounds derived from IceCube (see also footnote 8). Since we want to exploit the information collected by the experimental collaborations on the allowed partial decay widths of P , and we do not want to introduce too many free parameters in the benchmark choices we are going to illustrate, we decided to set c P P φφ = c HHφφ = 0 5 .
In order to explain the diphoton excess and ensure that it can be consistent with Run-1, we need at least c BB and one among c gg and y qP to be different from zero. In the benchmark scenarios that we illustrate below we always assume a coupling of P both to gluons and to light quarks. These allow respectively to improve the compatibility between Run-1 and Run-2 of LHC (given the stronger increase of gluon PDF with energy with respect to valence quarks), and to increase the elastic cross section between proton and DM, in order to achieve more easily the equilibrium of DM number density in the Sun.
The coupling of P to B B induces partial decay widths of P to the final states γγ, Zγ and ZZ, proportional respectively to cos 4 θ W , 4 cos 2 θ W sin 2 θ W , sin 4 θ W , where θ W is the Weinberg angle. We stress again that the compulsory coupling of P to γγ automatically implies, when requiring SU (2) L invariance, a coupling of P to ZZ and Zγ that can give a relevant signal for IC, given the typical hard spectrum of ν from Z decay.
We identify then three benchmark scenarios, where we always assume a coupling of P to DM, gluons, photons (and hence ZZ and Zγ), a light quark (for simplicity we assume that only y uP = 0, but the result would not differ much if also y dP and y sP were present), and possibly a heavy quark (top or bottom).
We choose the couplings c BB , c gg , y uP , y bP , y tP , and the coupling to DM in such a way to satisfy three requirements: 1) the corresponding partial decay widths of P satisfy the present collider bounds from Run-1 and explain the Run-2 diphoton excess [5], 2) the capture rate of DM in the Sun due to elastic scattering with the proton is high enough so that the DM number density in the Sun has reached equilibrium today, and 3) we obtain a reasonable compatibility between γγ searches in Run-1 and Run-2 6 .
The preliminary indication on the total decay width Γ P = 0.06 m P is included in most of the bounds shown in Ref. [5]. As we will see in section 3, the IC bounds depend only on the BR's of χχ (or φφ) into SM particles. With the use of Eqs. (A.13) and (A.14) (or (A.15) in the scalar DM case) one can rewrite the BR into, say, γγ, as follows, if P couples for example to BB, gg and uu: As stressed in appendix A, the dependence on Γ P and on the coupling of P to DM disappears rendering the IC flux predictions independent from both. Even more importantly, the ratios of partial widths appearing in Eq. (2.4) are the quantity directly constrained by LHC observations. Furthermore, the allowed ranges for these ratios do not alter significantly if one drops out the assumption Γ P = 0.06 m P (see for example Figs. 1 and 2 of Ref. [5]).
We now illustrate the three benchmark scenarios we have chosen 7 . 6 We quantify this latter requirement by means of the gain factors defined as the ratio r℘ between the cross section for the process pp → P → γγ at 13 and 8 TeV, computed assuming that the resonance is produced from a single parton ℘. Compatibility between Run-1 and Run-2 is achieved for r℘ ∼ 5, which favours ℘ = g, b, s, c and disfavours ℘ = u, d, γ. 7 For the rest of this section, for simplicity of notation we denote the DM particle by χ, but the discussion is identical for the case of scalar DM.     Table 3. Partial decay widths of P into SM channels for scenario C. Scenario A: P couples to B, g, u, χ. We choose the coefficients of the model leading to the partial decay widths of P into SM particles listed in table 1. The production of P at LHC at 13 TeV occurs mainly from gluons, while at low energies the elastic scattering with protons is mediated mainly by u. The resulting BR's for the annihilation of χχ into SM channels as a function of m χ are shown in Fig. 1. 8 The relevant channels for the production of hard neutrinos are ZZ and Zγ.
Scenario B: P couples to B, g, u, χ and b. The chosen partial decay widths of P are listed in table 2. The coupling of P to b quarks improves the compatibility between Run-1 and Run-2, although the coupling to photons is quite larger than the one to gluons. As far as LHC is concerned, the production of P occurs mainly from b partons. The resulting BR's for the annihilation of χχ into SM channels are shown in Fig. 2.
The bb channel yields soft neutrinos for IceCube, for which the main channels remain ZZ and Zγ.
Scenario C: P couples to B, g, u, χ and t. The chosen partial decay widths of P are listed in table 3. We choose the coupling to tt to be the main one. This channel does not contribute to the production at LHC, and amplifies the signal for IC. The production at LHC occurs via gluon fusion. The corresponding BR for χχ into SM channels are plotted in Fig. 3. Given its large BR, tt is the only important channel for IceCube.
We illustrate in the next section the main features of IceCube and its relevance in DM searches. Section 4 shows the resulting bounds from IC, DD experiments and γ-rays observations for the three scenarios just described.

Constraints on DM annihilations from the IceCube experiment
The IceCube experiment, located at the South Pole, is a neutrino telescope observing high energy neutrinos by detecting Cherenkov photons radiated by charged particles produced in their interactions [37]. Muons from ν µ and ν µ charged current interactions leave long visible tracks within the detector, which can be easily reconstructed to estimate the direction of the incoming neutrino. IceCube has an angular resolution of a few degrees for ∼ 100 GeV ν µ (and < 2 • for ∼ for a 1 TeV ν µ ), allowing it to search for an excess of GeV-TeV neutrinos from the direction of the Sun [38]. DM accumulates in celestial bodies if it loses some of its kinetic energy via elastic scattering and remains gravitationally bound within them. As the DM number density n DM increases, pairs of DM particles annihilate into SM ones, and some are lost due to evaporation (for further details, see Ref. [39] and references therein). The system eventually reaches an equilibrium, and n DM freezes. Neglecting the evaporation rate, which is a plausible assumption for m χ 10 GeV, at equilibrium the capture rate Γ cap is equal to 2Γ ann , where Γ ann is the annihilation rate of DM particles into SM particles and is proportional to n 2 DM . The capture rate Γ cap is proportional to the elastic scattering cross section σ pχ between a proton and the DM particle χ. The higher Γ cap , the faster the equilibrium is reached. When equilibrium is achieved, σ pχ can be directly related to the annihilation rate.
The key point to note is that once equilibrium is attained, in order to predict the neutrino fluxes searched for by IC only the ratios of the annihilation cross sections (or branching ratios) matter, and these ratios can be directly related to LHC measured quantities. Indeed, once we know the flux of neutrinos on Earth per DM annihilation, we can infer an upper bound on the rate of DM annihilation from the non-observation of an excess over the expected background. In other words, the neutrino fluxes depend only on branching ratios, not on the absolute value of annihilation rate.
The computation of the energy spectra of the neutrinos originated from the annihilation of DM particles is performed including the electroweak corrections, which can significantly alter the neutrino spectra when the mass of the DM particles is larger than the electroweak scale [40]. This is because soft electroweak gauge bosons are copiously radiated when the mass of the DM is larger than the gauge boson mass, and this opens new channels in the final states, including neutrinos, which otherwise would be forbidden if such corrections are neglected. In order to implement such electroweak corrections, we make use of the PPPC4DM ID code [41].
The current most stringent bounds on spin dependent DM-nucleon scattering cross section from IC comes from a search employing 3 years of IC data [35]. The IC collaboration presents constraints in terms of three benchmark cases: DM annihilation into bb, W + W − and τ + τ − . The least stringent constraints are for DM annihilating 100% into bb, a situation which produces few neutrinos and at energies much below the mass of the DM. DM particles annihilating 100% to W + W − or τ + τ − produce a significantly larger number of neutrinos at energies close to the DM mass and consequently stronger bounds on the scattering cross section. Neutrino flux predictions for these cases can be obtained from WimpSim [42]. The event rate expected from a differential (anti)-muon neutrino flux F(E) in an IC sample of effective area A eff is given by These quantities can be calculated for neutrino flux predictions calculated as described before, and also for the benchmark channel predictions, for each of the three IC samples described in Ref. [34]. Subsequently, the limit on the annihilation rate Γ ann for theoretical neutrino flux F(E) can be obtained by rescaling the benchmark limits according to the expression: 3) Figure 4. Upper limits on the DM annihilation rate in the Sun for the IceCube benchmark channels, as well as the three scenarios discussed above. The benchmark limits shown and used for rescaling conservatively correspond to the upper edge of the systematics band in [35] and are derived using the online flux conversion tool provided by the authors of WimpSim [43].
where Ψ(E) is the angular resolution of the IceCube sample at at neutrino energy E. This scaling is possible because the search is performed using the unbinned maximum likelihood ratio method [33,34], for which the sensitivity scales with Signal/ √ Background and the background level varies as Ψ 2 . In practice a scaling factor averaged over the three samples weighted by their exposure at the median energy is used. For a given F theory , σ theory can be calculated with respect to any of the three benchmark IC channels. The different calculations are consistent to within ∼ 30% and are thus averaged.
The limits on Γ ann for the IC benchmark channels can be obtained from the limits on σ using tools provided by WimpSim and DarkSuSy [43]. Figure 4 illustrates the limits on Γ ann for the scenarios discussed above, as well as for the IC benchmark cases. We notice that the exclusion bounds reported in Fig. 4 are not affected by the thresholds for the opening of new annihilation channels, such as m χ ∼ m top in scenario C for example. The reason is that the assumption of equilibrium in the DM number density in the Sun implies that IC constrains the DM capture rate in the Sun, rather than its annihilation rate. This is the reason why also the corresponding bounds in Fig. 6 and 7 do not display bumps in correspondence of kinematic thresholds.

Results and discussions
Having set our benchmark models describing the interactions of the pseudoscalar P with the SM particles and having addressed the way we deal with the IC physics, we now proceed to present our results in the cases in which the DM is a fermion and a scalar.

The case of fermionic DM
Using the solar capture rates evaluated in Ref. [44], the bounds on Γ ann shown in Fig. 4 can be interpreted as a bound on σ pχ for the operator O NR 6 = ( s χ · q)( s N · q). This operator originates in the non-relativistic limit in the case in which the DM is a Dirac fermion χ.
In comparison to direct detection, bounds from solar DM searches are generally more stringent in the scenario where the DM-nucleon scattering depends on the spin of the nucleus, since the Sun consists of mostly protons, in contrast to the target nuclei used in direct detection experiments which usually have no spin. Bounds on σ pχ for the operator ( s χ · q)( s N · q) are significantly weaker than for the common spin dependent operator due to the double velocity suppression.
The bounds from DD are obtained with the use of the software made available by the authors of Ref. [45], which allows to derive easily the bounds on a combination of non-relativistic operators for the interaction between DM and nucleons.
We show in Fig. 5 the upper bounds on the scattering cross section between proton and DM for the six experiments included in the software of [45], for the case of fermion and scalar DM. In the final comparison plots (Figs. 6, 7) we show only the convolution of these exclusion limits, corresponding to a combination of the bounds from LUX and XENON-100. Together with the current IC bounds obtained using the procedure explained above, we show a forecast for the bound that could be obtained by a similar neutrino telescope with 300 times the exposure. The sensitivity will scale with the square root of exposure [33]. While the exposure scales linearly with time and an improvement of ∼ 300 is unlikely for IceCube, future proposed neutrino telescopes such as KM3Net/ARCA [46] may achieve a similar improvement in sensitivity faster, using larger volumes, better angular resolutions or improved analysis techniques. However, it is not clear if future neutrino detectors will target the sub TeV energy range in primary neutrinos that is crucial for Solar DM searches.
We now consider the bounds from the observations of γ-rays in the sky. We recast the constraints from the searches for spectral lines in the spectrum of γ-rays from the centre of the Milky Way from Fermi-LAT [47] and HESS [48]. As benchmark choices for the DM density profile we select the Einasto and Burkert profiles [41], which present a high and a low density in the centre of the Galaxy, respectively. The two experiments are sensitive to different energy ranges, thus constraining complementary intervals in m χ . We show them with the same colour code in Figs. 6 and 7. We also examine the exclusion bounds from the the γ-ray continuum searches from a set of 15 Dwarf Spheroidal Galaxies (DSG) performed by Fermi-LAT [49]. A detailed description of our recasting procedure can be found in Appendix A.4.
We also show, as a tentative reference point, the expected σ pχ for a value of y χP = 1, from Eq. (A.3) and the lines corresponding to the points of the parameter space where the relic energy density of the DM candidate through the freeze-out mechanism turns out to be Ω DM h 2 = 0.1196 ± 0.0031 [50].
Our results are summarized in Fig. 6. Upper bounds on σ pχ from Direct Detection, IceCube and Indirect Detection experiments, for fermionic DM. The convolution of DD constraints is shown with a golden solid line. All the other constraints depend on the chosen scenario, which is identified by a specific dashing style. IC constraints are plotted with thick blue lines, and forecasts for a neutrino telescope with 300 times the exposure of IceCube are shown with a light blue band. With green and orange lines we show respectively the upper bounds from γ-ray lines and γ-ray continuum observations. To improve the readability, we use the same colour code for Fermi-LAT and HESS constraints on γ-ray lines which apply to different ranges of m χ , and we do not show the constraint for scenario C, being similar to the other two. The red lines show the prediction obtained by imposing that the relic density of χ and χ equals the observed one. We show also the expected signal for y χP = 1 (thin grey lines).
The conclusion we draw from Fig. 6 is that the double velocity suppression arising in the case of fermion DM and pseudoscalar mediator dramatically reduces the experimental exclusion reach of DD and IC. Stronger constraints from γ-ray observations turn out to fall in the region favoured by the calculation of DM relic density. This outcome further motivates us to investigate the case of scalar DM.

The case of scalar DM
In this case the non-relativistic operator for the interaction between the DM particle φ and nucleon is O NR 10 = i( s N · q), which is SD and suppressed by one power of DM velocity, instead of two as in the fermion DM case. This behaviour alters significantly the experimental bounds, as apparent from Fig. 7. As a reference point, we show the predicted cross section for the value A P = 1 TeV of the dimensionful coupling A P introduced in Eq. (2.3). The cross section σ pφ , reported in Eq. (A.5), is proportional to A 2 P . We conclude that, in the case in which the DM is a scalar particle, IC experimental bounds are important for two reasons: first, for masses of the DM larger than O(300) GeV, IC constraints are stronger than the bounds coming from DD experiments and, secondly, they bound A P to be not larger than a few TeV. Forecasts for a possible future neutrino telescope lower this bound to around 1 TeV in the region (m φ ∼ 0.1 ÷ 1) TeV, improving significantly the bounds from DD experiments in the mass region close to m P . The reach of the exclusion limits from the observation of γ-ray lines is affected by astrophysical uncertainties, mainly the DM density profile. In the region m φ 500 GeV, HESS yields are comparable or weaker to IC for DM profiles such as Burkert or Isothermal, and stronger for steeper profiles as Einasto. Upper limits from the observations of DSG, on the other hand, fall above IC ones.

Conclusions
Should future data collected by the LHC collaborations ATLAS and CMS confirm the existence of a new resonant state with a mass of around 750 GeV, the era of physics beyond the SM would start. Following this wishful route, one can imagine that the resonance acts a portal into the DM sector. By identifying three benchmark models, in this paper we have investigated the possibility that the resonance is caused by a pseudoscalar particle which also couples to either fermionic or scalar DM. Motivated by the fact that the bounds from solar DM searches are generally more stringent than DD experiments when the DM-nucleon scattering depends on the spin of the nucleus, we have analyzed the bounds coming from the search for neutrinos originated from DM annihilations in the Sun performed by the IC collaboration, and we have compared them with Direct Detection experiments and γ-ray observations. Our findings indicate that the IC data provide constraints stronger than the DD experiments for DM masses larger than a few hundred GeV. Furthermore, if the DM is a scalar particle, the IC data limit the cross section between the DM and the proton to values close to the predicted ones for natural values of the parameters.

A Relevant formulae for scattering with protons, annihilation of DM and computation of relic density
We report in this section the relevant formulae for the computation of cross sections for DD, for DM annihilations (which are relevant for IC), and of the relic density of DM.

A.1 Elastic scattering between proton and DM
Fermionic DM The effective operator for the interaction between χ and a nucleon N is (throughout this section we follow the notation of Ref. [45]) where the coefficient y χP is introduced in Eq. (2.2), and we define c N as are the nucleon spin form factors for the quark q (see Ref. [45] for a compilation of their numerical values according to various references).
The cross section for the elastic scattering χN → χN in the low velocity limit turns out to be (the same result holds for the process χN → χN ) where v DM is the DM particle velocity in the centre-of-mass frame (which we assume to be 220 km/s), and µ N = mχm N mχ+m N is the DM-nucleon reduced mass. We notice that this cross section is suppressed by the fourth power of the velocity of DM, as a result of the presence of the two pseudoscalar Lorentz bilinears in Eq. (A.1) which in the non-relativistic limit reduce each to the scalar product of the spin of the correspondent particle and of q.
Scalar DM The effective operator for the interaction between the DM particle φ and a nucleon is where A P is the parity-odd dimensionful coefficient introduced in Eq. (2.3), and c N is defined in Eq. (A.2). The corresponding cross section in the low velocity limit is We point out that in this case σ N φ is suppressed just by the second power of v DM , accordingly to the fact that O φN 2 contains only one pseudoscalar Lorentz bilinear.

A.2 Annihilation of DM pairs into SM channels
The squared matrix elements, summed over the polarisations of the final states, for the two-body decays of P are the following [13] (we denote by θ W the Weinberg angle): M P →qq 2 = 6 y 2 qP s , (A.10) M P →χχ 2 = 2 y 2 qP s , (A.11) The partial width of P into a final state composed of two particles ij is then given by where s ij is a symmetry factor equal to 1/2 if i and j are identical particles, or 1 otherwise. We can then write the cross section for the annihilation of χχ into a final state ij at a centre-of-mass energy s as (A.14) We highlight with a shade of grey the terms that are common to all final states, and simplify when computing the branching ratios for the annihilation of DM pairs into SM channels. These are indeed the relevant quantities for computing the neutrino fluxes for IC. This interesting simplification makes IC bound independent of the coupling of P to DM, of m P and of Γ P . The analogue formula to Eq. (A.14) for the case of scalar DM is: The centre-of-mass energy s, when applying Eq.s (A.3) and (A.5) to DM annihilations inside the Sun, has to be evaluated with the typical kinetic energy of DM particles in the Sun. We notice that all the squared matrix elements in Eqs. (A.6)-(A.12) in the low velocity expansion have a non vanishing constant term. We assume as a reference kinetic energy the thermal one inside the core of the Sun, around ∼ 1 keV (corresponding to a velocity 10 −4 for m χ ∼ 100 GeV).

A.3 Relic density of DM via freeze-out
We collect the main formulae needed to compute the actual relic abundance of DM. For a thorough discussion, we refer to Refs. [51,52].
By defining x = m DM /T , the expression for the thermally averaged cross section σv at a temperature T reads [51] σv = x 8 m 5 where K i is the modified Bessel function of order i. The relic abundance is then obtained from [52,53] where T f is the freeze-out temperature, T 0 is the present temperature, g is the degrees of freedom parameter as a function of the temperature and the factor of 2 accounts for the fact that the total DM density is the sum of the density of DM particles and antiparticles. Note that here we approximate g = g eff = h eff , where g eff and h eff are the number of degrees of freedom that enters the definition of the energy density and of the entropy density respectively [51]. The freeze-out temperature T f (or, equivalently, x f ) is defined to be the temperature at which the quantity Y = n/s differs from its equilibrium value by Y − Y eq = c Y eq . With a standard notation we indicate by n and s the number density of DM particles and the entropy density of the Universe. The value of x f is obtained by where g = 1 (scalar), 2 (spinor) accounts for the number of spin states of the DM particle. Following Ref. [51] we take c = 1.5.
The thermally averaged cross section appearing in the denominator of Eq. (A.17) can be expanded as a function of x. This approximation can be slightly inaccurate when a resonance is excited [54], because for m DM just below the resonance also higher orders of the expansion of σ ann (v) at low velocity can matter. For this reason, we perform our computation without expanding Eq. (A.16) in powers of v DM .
We report the results for the expansion of σ(χχ → ij) · v Møl in series of v Møl , to show that the term of order zero in the velocity is always non vanishing. Therefore, the final value for the cross section does not depend strongly on the numerical value assumed for DM velocity.
For the case of fermionic DM, the expanded cross sections read 23) In the case of scalar DM, the results are

A.4 Constraints from observations of γ-rays
Search for γ-ray spectral lines from the centre of the Milky Way In this section we provide further details about the procedure we adopted in order to recast the exclusion bounds on DM annihilations coming from the search for γ-ray spectral lines from the centre of the Milky Way. The annihilation of two DM particles 9 χ and χ into γγ, in the limit v DM → 0, provides two photons of energy E γ = m χ . The process χχ → Zγ gives only one photon, with The constraints provided by Fermi and HESS on σv however are derived from searches for monochromatic spectral lines. Nevertheless for m χ larger than 100 GeV, the relative separation in the energies of γs from the two processes is smaller than the energy resolution of the instrument and consequently, constraints can be derived using the following method. For m χ smaller than 100 GeV, the second term in the LHS of Eq. (A.29) is ignored to obtain a conservative bound.
The fluxes of photons obtained in our model for a DM mass m χ , to be compared with the exclusion limit σv limit γγ provided by Fermi-LAT or HESS which assume that only χχ → γγ occurs, is given by where the arguments of σv inside brackets refer to the DM mass to be inserted into Eqs. (A.20), (A.21). The Fermi-LAT analysis, performed with the data collected over 6 years by the Large Area Telescope hosted by the satellite Fermi [47], identifies different signal regions depending on the DM density profile under consideration. They select the region R16 (a cone with an opening angle θ = 16 o around the centre of the Milky Way) for the profile Einasto, and the region R90 (corresponding to θ = 90 o , i. e. half of the sky around the centre) for the profile Isothermal. We want to recast the limits assuming the DM profile Burkert, characterised by a lower DM density in the centre of the galaxy. Thus we compute the integral of the J-factors in the region θ ≤ 90 o with the tables provided by [55] (to which we refer for the definitions of the DM density profiles and of the J-factors), and we use it to rescale the bound computed by Fermi-LAT assuming the Isothermal profile.
The sensitivity of LAT to γ-rays of energies up to ∼ 300 GeV is complemented by the sensitivity of the telescope HESS, located in Namibia, which provides important bounds on the annihilation of DM into γ rays for m χ 500 GeV [48]. HESS observes a cone of 1 o around the centre of the Milky Way, and recasts the exclusion limits assuming the Einasto profile. Analogously to what we did in the previous case, in order to get the corresponding limits with a Burkert profile we rescale them by computing the ratio of the integrated J-factors.
Observation of the γ-ray continuum from Dwarf Spheroidal Galaxies Fermi-LAT performed [49] an analysis of the γ-ray spectrum coming from 15 DSGs. These particular galaxies are characterised by an higher density of DM particles than ordinary galaxies, thus they represent an ideal target to look for secondary γ-rays resulting from the primary products of DM annihilations.
A careful recast of the analysis of [49] would require the computation of the spectrum of γ-rays from DM annihilations (which can be done through the tables provided by [55]), and the knowledge of the exclusion limits on the fluxes of γ-rays. This information is not available, but we can make a good approximation in our scenarios by identifying, in suitable intervals of m χ , the leading annihilation channel providing secondary photons (for example, in scenario B we consider only the annihilation channels χχ → bb for m χ < 1.5 TeV and χχ → Zγ for m χ > 1.5 TeV, see Fig. 2). The exclusion limits are recast by [49] as limits on the thermally averaged annihilation cross section into some specific channels, assuming in every case that they are the only annihilation channels of DM. These include uu, which provides basically the same spectrum as gg, and W + W − , which yields analogous fluxes as ZZ, tt, bb, and the double of the flux of Zγ [55]. These considerations allow us to perform the recast of Fermi-LAT observations by imposing that our main annihilation channel (for a given scenario and range of m χ ) equates the exclusion limit of [49].