Complementarity of DM Searches in a Consistent Simplified Model: the Case of Z'

We analyze the constraints from direct and indirect detection on fermionic Majorana Dark Matter (DM). Because the interaction with the Standard Model (SM) particles is spin-dependent, a priori the constraints that one gets from neutrino telescopes, the LHC, direct and indirect detection experiments are comparable. We study the complementarity of these searches in a particular example, in which a heavy $Z'$ mediates the interactions between the SM and the DM. We find that for heavy dark matter indirect detection provides the strongest bounds on this scenario, while IceCube bounds are typically stronger than those from direct detection. The LHC constraints are dominant for smaller dark matter masses. These light masses are less motivated by thermal relic abundance considerations. We show that the dominant annihilation channels of the light DM in the Sun and the Galactic Center are either $b\bar b$ or $t\bar t$, while the heavy DM annihilation is completely dominated by $Zh$ channel. The latter produces a hard neutrino spectrum which has not been previously analyzed. We study the neutrino spectrum yielded by DM and recast IceCube constraints to allow proper comparison with constraints from direct and indirect detection experiments and LHC exclusions. Note that the original version of the paper contains an important error: the contribution of the $Z$-mediated processes was overlooked. This changes some of the results of the paper. For the details and the correct results please see the erratum attached to this file as an appendix.


Introduction
Although there is overwhelming evidence that Dark Matter (DM) comprises nearly a quarter of the energy budget of the Universe [1], the precise nature of the DM still remains unknown. Until now all the evidence for DM has come from gravitational interactions, while searches for other DM interactions have yielded negative results. However, nowadays there is a well-developed program of searches for DM via its presumed non-gravitational interactions. This program includes direct detection experiments, a plethora of various LHC searches, indirect searches for the annihilation of DM captured by the Earth and the Sun, as well as searches for DM annihilations in the Milky Way Galactic Center and in nearby dwarf Galaxies.
One of the main motivations for thinking about DM interactions beyond the gravitational is the idea of thermal freezeout. This idea finds its natural home in the so-called WIMP-miracle scenario: a particle χ with a mass of order ∼TeV, charged under a weak force, yields the measured relic abundance. Direct detection experiments have mostly been aiming to discover this particular DM candidate, which also looked very attractive as part of a solution to the SM naturalness problem, in particular via supersymmetry and extra-dimensional scenarios.
The strongest bounds on the spin-independent DM-proton scattering cross-section in most cases are coming from LUX [2] that constrains this scattering at the level of σ χp ∼ 10 −44 − 10 −45 cm 2 . These cross sections are much smaller than the natural cross sections one would expect from the WIMP miracle, and, in fact, they are more characteristic of "higgs-portal" DM [3,4]. These direct detection results are very suggestive and put the entire idea of WIMPs under pressure.
However, there are of course several noticeable caveats in this logic, with the most important being the possibility of spin-dependent (SD) DM. If the nucleon-DM scattering rate depends on the spin of the nucleus, the total cross section in direct detection experiments is reduced by orders of magnitude. The strongest direct detection bounds on SD DM, coming from PICO [5], constrain σ SD χp 10 −38 − 10 −39 cm 2 , perfectly within the naive ballpark of the WIMP DM. LUX constraints on the SD DM are comparable to those of PICO [6].
Interestingly, at least naively, the direct detection bounds are not the strongest bounds one can put on the SD nucleon-DM scattering. Even stronger exclusions have been claimed both by Fermi-LAT, IceCube and by collider searches at the LHC. However, these constraints are more model dependent and the comparison of these constraints to the direct detection operators is not straightforward.
The problem with comparing these results is twofold. Let us first consider the signal in neutrino telescopes. Assuming that the DM capture and annihilation rates are already in equilibrium in the Sun's core, the total neutrino flux is only proportional to the capture rate, and independent of the DM annihilation rate [7]. However, in order to know the predicted neutrino flux, one should also know the annihilation channels. As we will discuss later, these are highly model dependent. The IceCube collaboration reports its results [8,9] considering annihilations into W W , τ τ (optimistic scenarios with hard neutrinos) and bb (pessimistic scenario with soft neutrinos). However, there is no guarantee that in a full UV complete picture any of these channels is dominant. In fact, if the DM is a Dirac fermion, it will most likely annihilate into a pair of SM fermions, and assuming flavor universality, these will most likely be light flavors, leaving no distinctive signature at IceCube. On the other hand, if the DM is Majorana and the fermion current preserves chirality, the light fermion channels are velocity-suppressed, giving way to other channels: W W , ZZ, hh, Zh and tt. A priori, there is no way to know which of these channels dominates the annihilation process and sets the neutrino spectrum. The same problem holds for the bounds obtained by indirect detection of the diffuse γ-rays from the Galactic Center and the dwarf satellite galaxies. These bounds also depend on the choice of the primary DM annihilation channels that give rise to the final spectrum of diffuse γ-rays.
The second problem has to do with the comparison of the direct detection exclusions to the LHC searches that nominally claim the strongest bounds on SD DM interactions with nucleons. However, this claim is strongly model-dependent. The LHC searches are extremely sensitive to the ratio between the mediator couplings to the DM particles and to the SM particles. Moreover, the collider exclusions quickly weaken in the high DM mass region.
In this paper we analyze the complementarity of all these searches in one particular, but highly motivated scenario. Of course the question of complementarity is meaningless if done purely within the effective field theory (EFT) of DM. Firstly, the EFT may not include some relevant annihilation channels of DM into weak gauge bosons and the Higgs, if the particles integrated out directly interact with them. Secondly, the EFT of DM is inadequate for analyzing the results of LHC searches unless the mediator is very heavy, otherwise the typical momentum running in the diagram is at least comparable to the mass of the mediator [10][11][12][13][14][15]. This leads to the inevitable conclusion that the question of complementarity can only be addressed in the context of full models, even if these models are simplified.
In this particular paper we will address the question of the complementarity between searches in the case that the interaction between the DM and the baryonic matter is mediated by a heavy vector boson (Z ). This is probably the easiest scenario one can consider beyond the proper WIMP (namely, the DM charged under the SM SU (2) × U (1)). It can also also be part of more complicated scenarios, e.g. non-minimal SUSY or strongly coupled physics at the TeV scale. We explicitly calculate the annihilation rate into the electroweak (EW) bosons and Zh and show that, if the annihilations into SM fermions are helicity suppressed, the dominant channels (depending on the DM mass) are bb, tt and Zh. The latter have been neglected in previous studies, however we show that for the heavy DM it is a dominant source of hard neutrinos, which can be detected by IceCube. Because even relatively small BRs into the gauge bosons that can potentially be important sources of the energetic neutrinos that can be detected by IceCube, we calculate the rates into the EW gauge bosons at the NLO level. We further comment on this issue in Sec. 3, while the details of the calculations are relegated to Appendix A. We show that the one-loop corrections to the W W annihilation channel are indeed important and noticeably harden the neutrino spectrum that one would naively estimate at the tree level. We also reanalyze the collider constraints, both on the direct production of Z and on the production of DM associated with a jet and review the situation with the indirect detection from dwarf satelite galaxies.
Our paper is structured as follows. In Sec. 2 we review the basic features of SD DM, we describe in detail our setup and we review the constraints from direct searches for a Z at colliders and from the calculation of the relic density of DM. In Sec. 3 we illustrate the constraints on our model from the LHC, direct and indirect DM searches and from the IceCube experiment, and we discuss the branching ratios for DM pair annihilation into the SM channels (which are required in order to derive bounds from IceCube and Fermi). In Sec. 4 we illustrate our results, and in Sec. 5 we summarize our conclusions. Appendix A contains further details about the formulae for the annihilation of DM.

Z -Mediated Spin-Dependent Interactions
Spin dependent DM-baryon scattering is a fairly generic phenomenon. The non-relativistic (NR) operators that describe this scattering [16][17][18] have explicit dependence on the nucleus spin. These interactions can arise from different types of UV theories and studying all possible SD interactions case by case would a tedious and unilluminating exercise. However, an S χ · S N interaction, where S χ , S N are the spin of the DM particle and the nucleon respectively, is in many senses unique in the long list of NR SD operators. It is the only operator arising from an interaction mediated by a heavy particle, for which the NR limit of the DM-nucleon interaction is not also suppressed by halo velocity; that is, either powers of q or v in the language of the NR EFT. These extra suppressions often render the direct detection and solar capture rates too small for most practical purposes.
In the high-energy EFT this operator descends from the axial vector -axial vector interaction [19], namely L = χγ µ γ 5 χ N γ µ γ 5 N . (2.1) As we have described in the introduction, this relativistic EFT description is not very useful, if we would like to compare direct detection searches for DM with the constraints from the LHC and IceCube. The UV completion for this operator is needed. There are several options for how this interaction can be completed at the renormalizable level. The simplest and most economical possibility, from the point of view of model building, is to mediate the interaction via a massive neutral vector boson, which couples to the SM axial current and to the DM axial current. This might either be a Z-boson of the SM (corresponding to a standard WIMP scenario) or a new Z boson. The simplest realization of this scenario would be a "classical" WIMP, namely a Majorana fermion that is charged under the SM SU (2) × U (1) (but not under the electromagnetism) and couples to the heavy weak gauge bosons at the tree level. One can view a pure supersymmetric electroweakino as a prototype of this kind of DM -either a doublet or a triplet of SU (2) L , possibly mixed with a singlet (bino). 1 SUSY dark matter also satisfies the criteria that we have listed above: the direct detection cross-section is spin-dependent 2 and therefore the bounds from direct detection experiments are relatively weak. Because the generic SUSY DM candidate is a Majorana fermion (rather than Dirac), annihilations into the light quarks are helicity suppressed. In the neutralino case, the annihilation into the EW gauge bosons proceeds at the leading order and in this case the W W , ZZ and tt channels dominate the annihilation branching ratios. The only relevant NLO annihilation channels for the neutralino spin-dependent DM are γγ and γZ, which have been fully analyzed in Refs. [26,27].
Another straightforward yet less explored UV-completion is mediation of interaction (2.1) by a heavy Z , remnant of a new gauge symmetry at the TeV scale. In particular, we assume the following renormalizable interaction between the DM, Z mediator and the SM fermions:  [28]). The scenario of DM which couples to the SM via a TeV-scale Z has been addressed in a number of references e.g. , however it is still unclear how various direct and indirect searches compare in this case. The most important caveat in this comparison has to do with the signal from neutrino telescopes, because it is very sensitive to the annihilation branching ratios of the DM into the EW bosons. The annihilation cross sections of DM into the SM fermions are helicity suppressed (except for m χ ∼ m f , with f a SM fermion, where there is no suppression into f f ), and therefore bosonic channels are expected to be the dominant source of neutrinos that can be detected by IceCube from DM annihilation in the Sun.
Parenthetically we note that another option for mediating the axial vector is via diagrams with new scalars in the t-and u-channels [51] (see also [52]). In this work we will put these examples aside, because of the generic problem of flavor changing neutral currents that these scalars can mediate. However it would also be interesting to see what the dominant annihilation channels are in this case. 3 In this study we exclusively concentrate on the Z mediator in the s-channel, as it appears in Eq. (2.2). Of course not every Z with arbitrary charges would be a consistent UV-completion of the effective contact term (2.1). If the Z is a gauge boson of a broken symmetry at multi-TeV scale, we should make sure that we are analyzing an anomaly free theory. This is a crucial demand, because without specifying the full matter content of the anomaly free theory, the annihilation rates of DM into the gauge bosons remain, strictly speaking, incalculable. Indeed, some of them could occur via loop diagrams, whose value is only guaranteed to be finite if the theory is renormalizable.
Based on the demands of renormalizability and anomaly cancellation, one can relatively easily parametrize all possible flavor-blind U (1) models. Any such model will be a linear combination of the SM hypercharge and U (1) B−L (see Sec. 22.4 of [53]). In extreme cases we either get a Y -sequential theory or a pure U (1) B−L . As we will shortly see, the latter is quite a boring case for DM direct and indirect searches because it is halo-momentum suppressed. 4 Because the charges of the SM fermions under the new U (1) are a two-parameter 3 An obvious way to avoid dangerous FCNC would be to impose some flavor symmetries, similar to those one usually assumes in flavor-blind SUSY models. 4 So-called "anomalous" U (1) has also been discussed in literature, e.g. U (1) ψ [54]. Note that these models demand the introduction of spectator fermions to cancel the anomalies. One should necessarily take into account the contribution of spectators when calculating the DM annihilation branching ratios to avoid non-physical results. family, we will conveniently parametrize the generator of the new symmetry as where t Y and t B−L stand for the generators of the hypercharge and B − L symmetries respectively. To ensure that the fermion masses are gauge invariant, the SM fermion charges under the U (1) unambiguously determine the SM Higgs charge under the U (1) . For completeness we list all the charges under the gauge symmetries, including the U (1) , in Table 1. These charges have a strong impact on the DM phenomenology in this scenario. Because the SM Higgs couples to the Z , tree level couplings between the Z and the EW gauge bosons are induced after EW symmetry breaking. In this case Z mixes with the Z. This allows annihilations of the DM to EW gauge bosons at the tree level.

Direct constraints on Z from LHC searches
Here we review direct constraints on this Z from the LHC. In addition, we will consider monojet constraints on DM production in Sec. 3.2. The easiest way to spot a Z at a collider is via an analysis of the leptonic modes, unless they are highly suppressed. For these purposes we recast a CMS search for a narrow Z in the leptonic channel [55], which conveniently phrases the constraints in terms of and for the reference point we take σ(pp → Z) × BR(Z → l + l − ) = 1.15 nb at √ s = 8 TeV [56].
We show the results of our recast in Fig. 1 as a function of the mass m Z of the Z and the angle θ as defined in Eq. (2.3). Note that, for a given g Z and θ, all couplings to the SM are fixed, therefore both production cross sections and BRs are unambiguously determined by these values. In Fig. 1 we show the exclusions both for nominal cross sections, as we get from MadGraph5 [57], as well as for cross sections one gets by applying the flat kfactor k = 1.3 (similar to the suggestion in [58]). It is apparent from these lines that, in order to have a Z with O(1) couplings, its mass must be m Z 3 TeV. Note that the constraints from LEP, coming from the mixing between the Z and Z , which further affect the T -parameter, are much weaker than the direct LHC constraints.

DM annihilation to SM particles
The annihilation proceeds via Z into the SM fermions and the EW gauge bosons. The couplings of the Z to the EW gauge bosons, in particular to Zh and W + W − , arise at the tree level after EW symmetry breaking, when the Higgs Goldstone modes are "eaten" by the massive gauge bosons. Before EWSB one can write down the couplings of the Z to the SM as follows: with D µ = D SM µ − i 2 cos θ g Z Z µ and the vector/axial vector couplings of the Z to the SM fermions given in Table 2. The Lagrangian describing DM is where g χ is the coupling of χ to the Z . In this paper we do not consider kinetic mixing, and we neglect the effects of renormalization group equations that could mix Z and Z via loop effects. We also do not take into account the running of the operator coefficients due to the RG flow because the quantitative effect is expected to be mild (see recent Ref. [59] for the details). The latter is indeed a minor effect, if the model is such that the mixing is zero at a scale close to the electroweak one. We do not specify the dynamics of the spontaneous symmetry breaking sector of U (1) , and for our purposes we just assume that it provides a mass term 1 2 m 2 Z Z 2 . As we expand Eq. (2.5) we find that the Z mixes with Z, with a mixing angle ψ fully determined by the mass m Z of the physical mass eigenstate, the value of the coupling g Z and the angle θ. If we denote the mass of the lighter mass eigenstate by m Z = 91.2 GeV, then ψ turns out to be of order ψ ∼ g Z cos θ m Z m Z 2 1 in the regime m Z m Z that we consider in this work. For this reason, in the remainder of the paper we ignore the mixing for simplicity of notation, and denote both the interaction or mass eigenstate equivalently by Z , and the lighter mass eigenstate identifiable with the SM vector boson by Z . Table 2. Coefficients of the vector and axial vector bilinear currents for the SM fermions (cf Eq. (2.5)). With obvious meaning of the notation, the coefficients g V and g A are obtained from Due to the Z-Z mixing, the heavy Z couples at tree level to Zh with a vertex where g Z = g 2 + g 2 = 2m Z /v with v = 246 GeV, and with W + W − with a vertex equal to sin ψ times the SM vertex for ZW + W − , Notice that both (2.7) and (2.8) are proportional to cos θ, thus they vanish in the pure U (1) B−L limit. Interaction (2.8) turns out to be velocity suppressed, thus it gives a small cross section in the low DM velocity regime of DD and IC, and both (2.7) and (2.8) are proportional to cos θ, thus vanish in the pure U (1) B−L limit. Because of these suppressions, loop channels are also relevant at low kinetic energy, as discussed in Sec. 3.5. In particular, the annihilation of χχ to Zh occurs at tree level, except for the case in which the U (1) extension is a pure U (1) B−L gauge symmetry. At low velocity, annihilation into W + W − is instead dominantly driven by diagrams with a fermionic loop (see Appendix A).

Calculation of DM relic density
Now that we have all the tree level couplings of the Z , we are ready to calculate the thermal relic abundance of the Majorana DM. We calculate the annihilation rates and perform the thermal average using the procedure of [60]. The result for the thermally averaged self-annihilation cross section as a function of temperature T is where x = m χ /T , and K i is the modified Bessel function of order i. We do not approximate the thermally averaged cross section with a low DM velocity expansion, since close to the resonance m χ m Z /2, terms of higher order can also yield important contributions to the relic density [61].
Once we fix the values of the angle θ and m Z , we are left with g Z as the only free parameter. In Fig. 2 we show the value of g Z that yields the correct relic density as measured by [1]. The areas above (below) the lines correspond to points of the parameter space where the DM is under-(over-) abundant in the thermal scenario.
For the calculation of the relic density we relied only on tree level cross sections, which we confirmed are dominant at typical freeze-out energies. Varying the value of g χ (which for this computation was set to 1) while keeping the product g Z √ g χ fixed only very mildly affects the decay width Γ Z , and would have practically no effect on Fig. 2.
The lines in Fig. 2 stop at m χ = m Z : above that threshold, annihilation into Z Z opens up, and in principle one would need to specify the details of the spontaneous symmetry breaking sector of U (1) in order to compute the relic abundance precisely. This is not required for our purpose of understanding the plausible range of values for g Z that yield a relic abundance close to the observed one.

Overview of direct and indirect bounds
The constraints on the scenario that we have described above come from three different primary sources: direct detection, neutrino telescopes and the LHC. Since we assume a Majorana DM particle, all interactions that we get between the nuclei and the DM are either spin dependent or halo velocity suppressed, or both. We will comment on these interactions in detail in Subsec. 3.1. This particular feature renders the direct detection results much less efficient than in the case of spin-independent interactions, while other experiments, notably neutrino telescopes and the LHC, become competitive. In the following section we will carefully go through all of these different experiments and discuss the bounds that they produce.

Direct Detection Experiments
In this type of experiment, in a model with a Z mediator the effective DM theory is always valid, because the transferred momentum never exceeds hundreds of MeV, well below the mediator scale. The effective contact terms that one gets between the DM and SM quarks are with coefficients g given in Table 2, and 1/Λ 2 = g 2 Z g χ /m 2 Z . It is straightforward to translate these interactions to the more intuitive language of the NR effective theory. Using the dictionary of Ref. [19] and considering a nucleus N instead of the partons q we get For a generic Z (arbitrary θ angle in (2.3)) we get both axial current -axial current (AA) and axial current -vector current (AV) interactions with the coefficients being roughly of the same order of magnitude. However, the latter induces interactions that are halo velocity suppressed (O 8 ) and both nuclear spin and halo velocity suppressed (O 9 ). This usually renders the AV interactions smaller than the AA interactions for direct detection and solar capture, although we do include AV interactions when we derive our bounds. In fact, the halo-velocity suppression in O 8 is sometimes comparable to the suppression due to the spin-dependence in O 4 . The operator O 9 , which is both spin-dependent and halo velocity dependent via the exchange momentum q is completely negligible. The one-loop contributions may induce a spin-independent scattering cross section of the DM with the nucleons, but their quantitative impact is negligible in our model [62]. For direct detection, there is a special point in parameter space where the usual spindependent operator completely shuts down. This happens when our new gauge symmetry is exactly U (1) B−L , or, in our language, θ = π 2 . In this case the Z couples only to the vector current of the SM, and therefore the NR interaction operators proceed only via O 8 and subdominantly via O 9 .
We derive the exclusion bounds obtained by the experiments LUX, XENON100, CDMS-Ge, COUPP, PICASSO, SuperCDMS with the help of the tables made available by [63], and the bounds by PICO from [5]. The strongest constraints among them are shown in Figs. 7 and 8, together with the constraints from IceCube and monojet searches.

LHC monojet constraints
In this section, we explore the bounds on our U (1) model coming from LHC searches for DM, analyzing events with one hard jet plus missing transverse energy ( / E T ). Currently, the strongest exclusion limits are from the CMS analysis [64], which analyzes 19.7 fb −1 at a collision energy of 8 TeV.
Despite the fact that the use of EFT to investigate dark matter signatures through missing energy advocated in Refs. [65][66][67] is severely limited at the LHC energies [10,11,13,14], we can nevertheless use the EFT interpretation of the exclusion bounds, as it is consistent for m Z 2 TeV [10,11,13,14]. The effective operators describing the interaction between DM and quarks, in the high m Z limit, are where Λ ≡ m Z /(g Z √ g χ ) and the coefficients g V q , g A q are given in Table 2. At LHC energy scales, the occurrence of a vector or an axial vector current in the fermion bilinears in (3.4) and (3.5) does not affect the cross section for the production of DM. This is also apparent from the experimental exclusion limits reported in [64].
The CMS analysis recasts the exclusion bound as a function of the coefficient Λ of Eq. (3.5), assuming that 1) all the g A q coefficients are equal to 1, and 2) that χ is a Dirac fermion (with canonically normalized kinetic term). These two assumptions are not true for our analysis. The second assumption gives an overall factor of 2 in the cross section for the Majorana case relative to the Dirac case. We take into account both the first assumption and the convolution with the parton distribution functions (PDFs) of quarks, by means of a parton level simulation performed with MadGraph 5 [57]. We simulate the signal both with the EFT defined by CMS and with our model, and we compute for each value of θ (which determines g V q , g A q ) and m χ a rescaling factor that we use to rescale the nominal limit reported by the CMS analysis.
The final result 5 for the bounds on Λ from monojet searches are reported in Fig. 3.

Constraints from observations of γ-ray spectrum
We now examine the exclusion bounds that can be obtained from the analysis of the γ-ray continuum spectrum. Limits coming from γ-ray lines are irrelevant for our model because the γγ, Zγ and hγ channels are strongly suppressed. The most stringent and robust bounds on the γ-ray continuum spectrum come from the observation of a set of 15 Dwarf Spheroidal Galaxies (dSph) performed by Fermi-LAT [68,69]. The robustness of these bounds against astrophysical uncertainties comes mostly from the fact that the photon flux is integrated over the whole volume of the dwarf galaxy, and in this way the inherent uncertainty due to the choice of the DM profile is largely diluted (as a reference, results are presented for the Navarro-Frenk-White profile [70]). However, the bounds are practically independent on the profile choice and the variation of the bounds due to J-factor uncertainties typically does not exceed 30% [69]. 5 We remark the following. If the Z mass m Z is larger than a few TeV, the bounds shown in Fig. 3 and the following fall in a region in which the product g Z g 1/2 χ is necessarily 1. This is in contrast with the fact that our Z model has a rather large mediator width, and it must be g Z g Here we should also briefly comment on HESS searches for DM using diffuse γ-rays from the Galactic Center. These bounds, claimed by the HESS collaboration [71], are nominally much stronger than Fermi dSph bounds for heavy DM, m χ 1 TeV. However one should also consider the uncertainties on these bounds. Unlike Fermi-LAT searches for emission from the Galactic Center, which mask a large region around the Galactic Center, 6 HESS merely masks a tiny region of 0.3 • around the Galactic Center, mainly to avoid the cosmic ray photon background, which is of course not present for the spacebased Fermi-LAT. This makes the search much more vulnerable both to the astrophysical uncertainties and to the choice of the DM profile. HESS assumes cuspy DM profiles in its search (NFW and Einasto), and if the profile is cored, the bounds can be attenuated by orders of magnitude. This point was nicely illustrated in the context of a different (wino) DM candidate in Refs. [74,75]. Therefore we decided not to show HESS' bounds. In order to properly recast the results of [69] for our model two ingredients are necessary. The first one is a knowledge of the spectrum of γ-rays from DM annihilations, which can be computed using the tables provided in Ref. [76]. Results of this calculation are shown in Fig. 4, for three reference values of the DM mass. The second ingredient would be the exclusion limits on the flux of γ rays, information which is not provided by the Fermi-LAT collaboration. For this reason, we adopted a simplified recasting procedure. Firstly, we identified in each interval in m χ the leading annihilation channel providing secondary photons, and approximated the total annihilation cross section with the one into that particular channel (or, in the case of multiple relevant primary channels with a similar γ-ray spectrum, we considered their sum). Secondly, we used the results of [76] to compare the photon flux from our dominant primary channel to the benchmark fluxes, namely e + e − , µ + µ − , τ + τ − , uū, bb and W + W − , which are shown in the left panel of Fig. 5. We used the limit on σv from the channel with the most similar photon flux as the limit on our channel. Finally, the limit on σv is converted into a limit on Λ. Though rough, we expect our procedure to provide bounds with at least an order of magnitude accuracy on σ SD χp (which translates into a factor of 2 on Λ). As we will discuss in Sec. 6, there are two regions of interest, the first for m χ < m W and the second for m χ > m W . Leaving aside for the moment the peculiar case θ = π/2, for m χ < m W the dominant channels are bb and cc, which give a similar γ spectrum, so the Fermi-LAT limit on bb can be assumed. On the other hand, for m χ > m W the dominant annihilation channel is Zh, complemented by W + W − , ZZ and tt, all of which give a similar photon flux. Since the flux of photons in the Zh channel is similar (up to a factor 2) to that in the bb channel, we again picked the Fermi-LAT limit on the bb channel. In the peculiar case θ = π/2, the dominant channels are leptonic for m χ < m W and W + W − , ZZ for m χ > m W . Therefore, for m χ < m W we picked the τ + τ − channel (which, among leptons, gives the strongest bounds), while for m χ > m W we summed the W + W − and ZZ contributions and compared with the limits on the W + W − channel. Fig. 5 shows the result of this recast in the right panel.

Neutrino Telescopes -IceCube
A stringent constraint on the spin dependent WIMP-Nucleon cross section comes from IceCube, a Cherenkov neutrino detector in the deep glacial ice at the South Pole [77], through a search for neutrinos from WIMP annihilations in the Sun [78]. Let us briefly review the main points of the physics related to DM annihilations in the Sun. We refer to Ref. [79] for a more detailed discussion.
DM particles in the Galactic halo can scatter with atomic nuclei inside the Sun and lose some of their energy. If the final velocity of the DM particle is low enough, it can remain gravitationally bound in the Sun. This accumulation process is counterbalanced by the evaporation process, i.e. the escape of DM particles from the Sun, and by the annihilation of DM pairs 7 into SM particles. The evaporation process can be safely neglected, if the DM particle mass is above ∼ 10 GeV. The interplay between DM capture and annihilation drives the DM number density n DM towards an equilibrium. The final value of n DM , and the equilibrium time by which this value is attained, depend on the capture rate Γ cap and on the annihilation rate Γ ann , the latter being proportional to n 2 DM . When equilibrium is attained, then Γ cap = 2Γ ann , and the annihilation cross section can be directly related to the cross section of elastic scattering between DM and proton.
Neutrinos are among the final annihilation products of DM, even more so when including the electroweak (EW) corrections [80]. If the energy of the final states is above the weak scale, then electroweak interactions imply a considerable emission of EW bosons in the final state, further amplifying the neutrino flux.
We include the effects of EW corrections and of the propagation of neutrinos by means of the PPPC 4 DM ID code [76,81]. This code provides a set of interpolation functions for the neutrino flux at Earth due to annihilation in the Sun, for a range of DM annihilation channels, and includes the effects of electroweak corrections. The code takes into account the cascade of the primary annihilation products into neutrinos within the Sun, as well as the subsequent neutrino oscillation in the Sun, in vacuum, and within the Earth.
What matters, in particular, are the branching ratios (BR) of DM pairs into SM final states: once equilibrium is achieved, the neutrino flux at Earth depends only on the relative annihilation cross sections to SM particles. This property is particularly interesting from the phenomenological point of view, because in the computation of the branching ratios, interesting simplifications occur and even the excitement of a resonance may have only a modest impact on the neutrino flux (see Sections 3.5, 4 and Appendix A).
Since the flux and spectrum of neutrinos from WIMP annihilations depend upon the preferred annihilation channels of the WIMPs, these constraints have traditionally been provided for extreme scenarios of WIMPs annihilating 100% into W + W − or τ + τ − , 'hard' channels corresponding to many high energy neutrinos and consequently a very stringent bound on the σ SD χp , or bb, a 'soft' channel producing a very weak bound on σ SD χp [8]. These results can be recast for the scenario of known annihilation channels by the following method. The search utilizes the Unbinned Maximum Likelihood ratio method [82,83], for which the sensitivity improves as Signal/ √ Background. For an unbinned maximum likelihood search of variable resolution, the background level varies as Ψ 2 where Ψ is the median angular resolution [84]. For a given differential (anti)-muon neutrino flux F(E), the total number of signal events expected within a sample can be calculated as where A eff is the effective area from Fig. 3 of Ref. [83]. Since the fluxes and effective areas of ν µ andν µ are different, Eq. 3.6 has to be evaluated separately for ν µ andν µ . The median energy E med (F) is then defined through If the capture and annihilation processes in the Sun are in equilibrium, the neutrino flux, the capture/annihilation rate, as well as WIMP-nucleon scattering cross section all scale linearly with respect to each other. Thus the theoretical limit on the WIMP-nucleon scattering cross-section for a given flux prediction F theory can be derived as where the first term in the RHS accounts for the variation in the level of signal events while the second term accounts for the variation in background due to the shift in median angular resolution. An analogous scaling relation can also be used to obtain theoretical limits on the annihilation rate Γ ann . The bounds on Γ ann for the IceCube benchmark channels can be derived from the limits on σ by mean of the tools provided by WimpSim and DarkSUSY.
The IceCube limit on the neutrino flux F limit requires knowledge of the neutrino spectrum per annihilation as it would be observed at Earth. The first step is to calculate the branching ratio to all relevant final states. Results are discussed in Sec. 3.5.
In order to convert these branching ratios into the required neutrino spectrum, we use the PPPC 4 DM ID code. This is combined with the results for the branching ratios to determine the final spectrum of muon neutrinos and antineutrinos per DM annihilation event. The Zh, γh and γZ final states are not available in the code, and so we use the average of the two pair-production spectra for each of these final states. We assume that the differences in the kinematical distributions, due to the different masses of Z, h and γ, have a minor impact on the shape of the final neutrino flux.
For the theoretical flux prediction thus obtained, the number of expected signal events as well as the median energy can be obtained from the expressions in Eqs. (3.6) and (3.7) for each of the three IceCube samples described in Ref. [83]. The integrals are evaluated separately for ν µ andν µ . These quantities can also be evaluated for the IceCube benchmark channel flux predictions (F benchmark ) obtained from WimpSim 3.03 and nusigma 1.17. Subsequently, theoretical limits on σ and Γ ann can be obtained using the scaling relation (3.8) and the analogous one for Γ ann .
For a given F theory , σ theory can be calculated w.r.t any of the three benchmark IceCube channels. The different calculations are consistent to within ∼ 30% and are thus averaged.

Results for the Branching Ratios
As discussed in the previous section, in order to extract bounds from IceCube observations, the branching ratios for the annihilations of DM pairs into SM final states must be computed. In this section we present the final results for the BR's, based on the cross sections that are computed in detail in Appendix A. Fig. 6 shows the BR's into SM two-body final states as a function of m χ , for m Z = 10 TeV and for six different values of θ (defined in Eq. (2.3)). The BR's for the final states shown on the plots do not depend at all on g Z and g χ , because they cancel in the ratios of cross sections, as can be seen from the formulae in Appendix A. Leaving aside for a moment the pure U (1) B−L , the main annihilation channels are the heavy fermion pairs (tt, bb and τ + τ − ) and Zh. These are indeed the only tree level channels for which a ≈ 0 in the low velocity expansion σv DM ∼ a + b v 2 DM . The main annihilation channel below the kinematic threshold m χ = 108 GeV for Zh production is bb, while the BR's into cc and τ + τ − are less than 10% each. Even if its branching ratio is not the dominant one, the τ + τ − channel dominates the IC bound below the Zh threshold because it yields more energetic neutrinos than the bb one.
In the region m χ ∼ 80 GeV − 108 GeV annihilation into W + W − may overcome the one into bb, depending on the value of θ. Notice that, as will be explained at the end of this section, the one loop contribution to the W + W − cross section dominates over the tree level one, which is suppressed by the small Z − Z mixing angle and by the fact that, in the low velocity expansion, it has a = 0.
When the Zh channel opens up it overcomes all others and remains the only relevant channel unless m χ m top , where the cross section into tt is comparable to the one to Zh. At higher DM masses, the cross section to Zh in the low v DM limit is proportional to m 2 χ , while σ(χχ → tt) is basically constant in m χ . The former proportionality comes from the final state with a longitudinally polarised Z boson and a Higgs boson, and is ultimately due to the derivative coupling of would-be Goldstone bosons. This explains why Zh is the only relevant channel at large m χ .
Around the resonance σ(χχ → Zh) goes to 0 because the coefficient a in the low velocity expansion σv DM ∼ a + b v 2 DM vanishes. The reason is explained in Appendix A. Therefore, in a small window around the resonance, other channels dominate. The position of this window is basically the only way in which the branching ratios depend on m Z , as can be seen from Figs. 7 and 8.
The previous picture applies for all values of θ, except θ = π/2 which corresponds to the pure U (1) B−L case. In that case, there is no mixing between Z and Z , and the channel Zh disappears at tree level. Below the W + W − threshold m DM = m W , annihilation predominantly happens at tree level into fermionic channels. For m χ > m W the dominant channel is instead W + W − , with a O(10%) contribution from ZZ. Annihilation into these channels is due to a diagram with a triangular fermionic loop, as discussed in Appendix A. The fermion channels, in the zero velocity limit, have a cross section proportional to m 2 f c A 2 f , where m f is the fermion mass and c A f is the coupling of Z to the axial vector fermion bilinear. When U (1) is reduced to U (1) B−L the Z couples to the vector current only. Thus in this limit the σ(χχ → f f ) has a = 0. The coefficient b is not proportional to m 2 f (as it is a because of the helicity suppression), thus in the fourth plot of Fig. 6 the fermions contribute equally to the annihilation cross section (apart from a factor B 2 × (# colors) = 1/3 which penalises quarks with respect to leptons), unlike what happens for θ = π/2.
Let us conclude this section by explaining why the tree level contribution for χχ → W + W − has a = 0, which, together with the additional suppression by sin ψ, selects Zh as the main channel at low velocities. The initial state χχ, in the limit v DM → 0, has total angular momentum J = 0 and CP eigenvalue −1. The final state Zh is not a CP eigenstate, unlike W + W − . Now, a pair of vector bosons can have a CP eigenvalue −1 and a total angular momentum J = 0 only if they are both transversally polarized [85]. In this case, the tree level interaction Z W W (2.8) turns out to give a vanishing cross section. We notice that this argument does not apply to the W + W − and ZZ amplitudes when the triangular fermion loop is included (see Fig. 9). In those cases, the effective Z W + W − , Z ZZ vertices contain the terms where k 1 , k 2 are the four-momenta of the outgoing bosons (see also [86,87]). These terms lead to a = 0 in the cross section, making the W + W − and ZZ channels relevant despite of the loop suppression.

Summary of results
We show in Fig. 7 bounds on σ SD χp for the case θ = 0, comparing direct detection results with those coming from IceCube, LHC's mono-jet searches and γ-ray searches. Analogous results for different values of θ are presented in Fig. 8, in the plane Λ vs. m χ . As representative values of θ we choose θ = 0, π/4, π/2 and 3π/4. We did not consider θ = π because it is exactly equivalent to θ = 0. Since θ = π/2 is quite a peculiar point, we added two values of θ in its vicinity, namely 0.45π and 0.55π.
As is clear from Eq. 3.3 and Table 2, for a generic angle θ the DM-nucleon scattering is mediated by a linear combination of O 4 , O 8 and O 9 , with coefficients similar in magnitude. The contributions from O 8 and O 9 can be safely ignored for IceCube: given the composition of the solar environment, their nuclear form factors are between 100 and 1000 times smaller than the one for O 4 [88]. This is not the case for DD experiments, where the three operators give a similar contribution, and the scattering cross section is not exactly σ SD χp , which is defined as being given by the operator O 4 alone.
To obtain a bound in the usual σ SD χp vs. m χ plane, one needs to first obtain a bound on the total scattering cross section, which is easily done using the results of [63]. Then, this bound is translated into a bound on Λ = m Z /(g Z √ g χ ), using the expression for the Figure 6. Branching ratios for the annihilation of χχ into pairs of SM particles, for m Z = 10 TeV and a kinetic energy of the DM particles equal to the thermal one in the Sun core. There is no dependence on g Z and g χ . The six plots, from left to right and top to bottom, correspond to θ = 0, π/4, 0.45π, π/2, 0.55π, 3π/4. Only the channels with a BR greater than 10 −3 are shown.
cross section obtained from Eq. 3.1. Finally, the bound on Λ is translated into a bound on σ SD χp using the expression which gives the scattering cross section when only the O 4 operator is involved, where µ p is the reduced mass of the proton-DM system, and the dimensionless coefficient is given by cos θ · 1.35, with ∆ q parametrizing the quark spin content of each nucleon, and is assumed to be equal for protons and neutrons [63].
For θ = 0, we do not show bounds on σ SD χp but only on Λ. The reason is that, since σ SD χp Figure 7. Bound on σ SD χp from direct detection, LHC's monojet analysis, IceCube and Fermi-LAT, for θ = 0.
is defined as the contribution to the scattering cross section due to the operator O 4 only, given the limit on Λ we have σ SD χp ∝ cos 2 θ/Λ 4 lim . When θ gets close to π/2, the computed value of σ SD χp goes to 0 independently of Λ lim , resulting in a spuriously strong bound. The situation is slightly different when θ = π/2. In this case, the coupling of the Z to the vectorial current of the quark fields is identically 0, and therefore the coefficient of the O 4 operator in the NR expansion vanishes. While for IceCube the contribution of O 9 is subdominant with respect to that of O 8 and can be ignored, both of them have to be taken into account to obtain DD bounds. A recast of the bound on Λ in terms of the usual σ SD χp would make no sense in this scenario. The PICO experiment gives interesting direct detection bounds. For θ = 0 (i.e. for the usual spin-dependent operator O 4 ), bounds on σ SD χp can be read directly from Ref. [5], and translated into a bound on Λ. Bounds on Λ for other values of θ can not be obtained following the procedure we adopted for the other experiments, because PICO is not yet included between the Test Statistic functions given in [63]. Therefore, we obtained a conservative bound on Λ by rescaling the PICO limit on σ SD χp as cos 2 θ (since this enters the coefficient of the O 4 operator) and then applying Eq. (4.1).
Let us now comment briefly on the results shown in Figs. 7 and 8. Bounds coming from LHC searches are typically the strongest ones in the low mass region, up to m χ ∼ 400 − 700 GeV for large θ. IceCube searches have their maximal sensitivity in the region between a few hundred GeV and a few TeV. When θ is small, in this mass region they give a constraint on Λ which is stronger than the direct detection one. In particular, for θ = 0 the constraint from IceCube is the dominant one from m χ ∼ 100 GeV up to 10 TeV and beyond. The Fermi-LAT bound from dwarf spheroidal galaxies appears to be the dominant one for m χ mass above 200 − 300 GeV. This is mainly due to the fact that, in our model, the dominant annihilation channel is Zh, which produces a large photon flux thanks to EW corrections. We expect that, in a leptophilic model in which DM annihilates copiously into τ and µ pairs, energetic neutrinos produced in their decay would make IceCube bounds the dominant ones.
We notice that the bounds shown in Fig. 8 fall in the region where the DM is underabundantly produced via the freeze-out mechanism (compare with Fig. 2). The difference between the two values of Λ goes up to one order of magnitude for large m χ .
An important remark about IceCube results is that they are almost independent of m Z and Γ Z (together with g Z , g χ , as stressed in Sec. 3.5). There are two reasons for this: first, as explained previously, when the equilibrium between annihilation and capture in the Sun is reached, neutrino fluxes at Earth only depend on branching ratios, which are not modified dramatically by the resonance except for a very narrow region around it. Electroweak corrections further dilute the difference, and as a result IceCube bounds for different values of m Z are almost superimposed (except for a small bump around the resonant point, whose width is ∼ Γ Z ).
Previous studies of the constraints coming from the IceCube experiment were done in Refs. [89,90]. In these works, the authors examined simplified models without considering the mixing of Z and Z imposed by the Higgs boson charge under U (1) , and therefore they do not include annihilation channels which turn out to be the dominant ones. Moreover, they weight the nominal IceCube benchmarks (which assume 100% annihilation into one channel) by the annihilation cross sections into different channels computed in their model, and they do not take into account EW corrections. In our work, we compute the BR's into various SM channels in a complete and consistent model, and we compute the neutrino fluxes including the EW corrections with the help of PPPC4DM ID. We infer the exclusion bounds with a recast of IceCube limits, as explained in Sec. 3.4, using the full shape of the neutrino fluxes to obtain the new bound.

Conclusions
While the spin-independent WIMP scenario has been probed experimentally to very high precision and direct detection experiments basically disfavor this possibility, the bounds on SD DM are much milder. WIMP-strength interactions between spin-dependent DM and baryonic matter are still perfectly allowed. Moreover, often the strongest bounds on SD DM do not come from direct detection experiments, but rather from LHC searches (direct or indirect) and Fermi-LAT. The IceCube experiment also produces interesting bounds, which are typically stronger than the direct detection ones above a few hundred GeV.
In this work we analyzed and compared these constraints, coming from different experiments. In order to be concrete, we concentrated on a particular set of spin-dependent DM models, in which the DM-baryon interaction is mediated by a heavy gauge Z . If we restrict ourself to the models without spectators at the EW scale, the parameter space of gauge Z models can be conveniently parametrized by three quantities: the mass of the heavy Z , the gauge coupling and the mixing angle θ between the hypercharge and U (1) B−L generators.
Although this study is not completely generic, as it does not cover all possible consistent models of SD DM, there are good reasons to believe that these models capture important phenomenological features that are generic to WIMP-like SD DM. We identified the region of parameter space favored by the observed thermal relic abundance and we have reanalyzed the existing LHC constraints, both direct (from the monojet searches) and indirect, on the Z mass and coupling.
More importantly, we fully analyzed the low-temperature annihilation branching ratios of the Z -mediated DM, which is crucial to derive IceCube (and Fermi-LAT) constraints. In order to properly understand the expected neutrino fluxes from DM annihilation in the Sun's core, one has to know the annihilation channels of the DM in the Sun. Simply assuming that the dominant expected source of neutrinos is the W W or bb channel (as is done in the IceCube papers) is clearly insufficient. We found that, depending on the DM mass, one can divide the parameter space into three different regions. Below the Zh mass threshold (very light DM) the annihilations are indeed dominated by the bb channel, yielding very soft neutrino fluxes, although even in this case the IceCube constraints are dominated by small branching ratio to τ + τ − . Above this threshold Zh is a dominant annihilation channel and the dominant source of neutrinos almost in the entire parameter space. The third region is just above the tt mass threshold. If the DM mass sits in this "island", one usually gets comparable annihilations into Zh and tt, and both should be taken into account for the neutrino flux calculations. The W + W − channel can also become important, or even the leading channel, if the U (1) extension is very close to being U (1) B−L . In this work we have properly recast the existing IceCube bounds including electroweak corrections, in order to derive reliable exclusion bounds on the secondary neutrinos coming from all annihilation channels.
We find that currently the strongest bounds on the SD Z -mediated DM are imposed by LHC searches (for m χ 400 GeV) and by Fermi-LAT for heavier DM candidates, which are favored by thermal relic considerations. We also find that the best direct detection bounds come mostly from LUX (PICO becomes dominant only for very light ∼ 20 GeV DM particles). These bounds are subdominant with respect to LHC ones in the case of light DM, but can be comparable to IC bounds for heavy DM if U (1) is close to being a U (1) B−L extension. We have also computed the values that yield the observed DM abundance through the freeze out mechanism, and we found that experimental exclusion limits fall in the slightly underabundant region.
Finally we notice that it would be interesting to see similar works for other SD DM candidates. It would be also useful to have bounds on the DM, annihilating into Zh reported directly by the IC and Fermi collaborations. We stress that this channel immediately arises when considering a consistent anomaly-free model for a U (1) extension which is not a pure U (1) B−L . This important feature is not captured by the use of so-called simplified models if these issues are not considered.

Acknowledgments
We thank Maxim Perelstein, Stefan Pokorski, Benjamin Safdi and Sofia Vallecorsa for discussions. We also would like to thank the anonymous referee for his/her suggestion of considering the γ-rays constraints on our model, without which our work would have been less complete. D.R. and A.R. are supported by the Swiss National Science Foundation (SNSF), project "Investigating the Nature of Dark Matter" (project number: 200020 159223). A.K. is grateful to the organizers of KITP program "Experimental Challenges for the LHC run II" where this paper has been partially completed (partially supported by National Science Foundation under Grant No. NSF PHY11-25915).

A Details of the annihilation rate calculation
In this appendix we present the calculation of the annihilation cross sections of the DM into the SM in detail. We also go in detail over the one-loop order annihilation into the W W . The results of this calculations are summarized on Fig. 6. The Feynman diagrams for the possible annihilation channels of χχ into the SM particles are shown on Fig. 9. For each channel, we consider the leading order (tree level or one-loop), moreover we always restrict the calculation to the leading order in the mixing angle between Z and Z , ψ. At the tree level the DM annihilates into ff , and if θ = π/2 also to W + W − and Zh. Eqs. (A.1) to (A.6) summarize the annihilation cross sections into all these channels at the tree level, distinguishing between the polarization of the vector bosons in the final states through a superscript (T ) or (L) for transverse or longitudinal polarization, respectively.
A few clarifications are in order about the annihilation cross sections of the DM at tree level.
f f : We denote the number of colors of the fermion f by N f c , and its vector and axial vector couplings to Z by g V f , g A f respectively. The values of g V f , g A f are given in Tab. 2 In the zero velocity limit, corresponding to s = 4m 2 χ , the cross section is proportional to m 2 f , because of the helicity suppression for pairs of annihilating fermions. In that limit, σ ∝ (g A f ) 2 ∝ cos 2 θ, i.e. the a coefficient in the low velocity expansion σv a + bv 2 comes from the U (1) Y component of the U (1) extension.
Zh: The diagram on Fig. 9 contains the tree level vertex Z Zh of Eq. (2.7). In the zero velocity limit, the only contribution comes from the production of a longitudinally polarized Z, since in Eq. (A.3) the factor (s − 4m 2 χ ) 1/2 vanishes.
W W : We denote α W = g 2 W /(4π), where g W is the weak coupling constant. The amplitude is the sum of the two diagrams on Fig. 9: the annihilation occurs via the mixing of Z and Z and the SM trilinear gauge vertex ZW W , see Eq. (2.8). For each of the final polarization states, the cross section is proportional to (s − 4m 2 χ ) 1/2 .
Given that σ(χχ → W W )v DM (where v DM is the DM velocity) is suppressed by sin 2 ψ and by v 2 DM at tree level it is worth checking whether contributions arising at one loop can become important. The contribution to the amplitude of diagrams 5 and 6 on Fig. 9 is velocity suppressed because of the same argument reported at the end of Sec. 3.5. Therefore we only considered the contribution to the cross section coming from the sum of diagrams 3 and 4, also ignoring the interference terms with the other diagrams.
We also computed the cross sections for the annihilations into ZZ, γγ, γh, γZ and hh at one loop. It is worth computing these corrections because of the velocity suppression of some tree level channels, and because in the pure U (1) B−L case the tree level annihilations into W W and Zh disappear. As for the Zh channel, we computed the one loop cross section only in the pure U (1) B−L case (θ = π/2), in which the tree level amplitude vanishes, and the only remaining contribution comes from diagrams 2 and 3 in Fig. 9. The results of the loop calculation are: We do not report the formula for σ(χχ → W ± (T ) W ∓ (L) ) because this channel turns out to have a = 0 (therefore it is irrelevant for the annihilation process in the Sun) and the corresponding formula is too cumbersome to be reported here. (A. 18) In the cross sections involving a fermionic loop, we denoted by f a sum over all the SM fermion species. In the W W cross section, instead, we denoted by i a sum over the six fermion families (three families of quarks and three of leptons), with m u i , m d i the upper and lower component of the doublet, respectively, and with g L The functions A 0 , B 0 and C 0 are the standard Passarino-Veltman one loop one-, two-and three-points scalar integrals [91]: , (A. 21) sections receive a velocity suppression (σv bv 2 ) and precisely vanish in the zero velocity limit. Those are W (T ) W (T ) , Z (T ) Z (L) , γZ (L) , γh, hh, Z (T ) h and, for θ = π/2, Z (L) h. We do not explicitly show the analytical expansion around v = 0 because the velocity appears as an argument of the Passarino-Veltman functions. The only process, which acquires a relevant contribution at the one-loop level is W (L) W (L) . All the cross sections we computed, except for hh that has no s-channel exchange of a Z boson, vanish around m DM = m Z /2 in the v DM → 0 limit. Again, this is a cross-check of the correctness of our calculations. Indeed, close to the resonance the cross section for χχ → XX is proportional to the product Γ(Z → χχ) · Γ(Z → XX), but Γ(Z → χχ) vanishes if m DM = m Z /2, which is implied by a resonant production of Z with v DM = 0.