ournal of C osmology and A stroparticle hysics J On the importance of direct detection combined limits for spin independent and spin dependent dark matter interactions

. In this work we show how the inclusion of dark matter (DM) direct detection upper bounds in a theoretically consistent manner can aﬀect the allowed parameter space of a DM model. Traditionally, the limits from DM direct detection experiments on the elastic scattering cross section of DM particles as a function of their mass are extracted under simpli-fying assumptions. Relaxing the assumptions related to the DM particle nature, such as the neutron to proton ratio of the interactions, or the possibility of having similar contributions from the spin independent (SI) and spin dependent (SD) interactions can vary signiﬁcantly the upper limits. Furthermore, it is known that astrophysical and nuclear uncertainties can also aﬀect the upper bounds. To exemplify the impact of properly including all these factors, we have analysed two well motivated and popular DM scenarios: neutralinos in the NMSSM and a Z ′ portal with Dirac DM. We have found that the allowed parameter space of these models is subject to important variations when one includes both the SI and SD interactions at the same time, realistic neutron to proton ratios, as well as using diﬀerent self-consistent speed distributions corresponding to popular DM halo density proﬁles, and distinct SD structure functions. Finally, we provide all the necessary information to include the upper bounds of SuperCDMS and LUX taking into account all these subtleties in the investigation of any particle physics model. The data for each experiment and example codes are available at this site http://goo.gl/1CDFYi, and their use is detailed in the appendices of this work.


Introduction
Nowadays, the dark matter (DM) nature is one of the greatest mysteries of modern physics. The possibility of non gravitational interactions of DM with ordinary matter is one of the key pieces to understand the DM problem. In this sense, the experimental community is putting a great deal of effort to provide such a valuable information.
Direct searches for DM aim at detecting the scattering of DM particles off nuclei inside a low-background target detector. Once a DM particle hits a nucleus, the recoiling nucleus will release a certain amount of energy that, above a given threshold, can be measured using different and sophisticated techniques. This field is currently undergoing an exciting period, with several experimental collaborations, using different targets and techniques, reporting potential signals of DM that are, nevertheless, in conflict with null results in other experiments. Most remarkably, the results obtained by the SuperCDMS [1] and LUX [2] Collaborations have placed the most stringent upper bounds on the elastic scattering spin independent (SI) cross section up to date, for DM masses above 4 GeV. Particle physics models for DM are then tested using these results, which in general, reduce considerably the allowed region of the parameter space of a given model.
We are living in the precision data era. There is a huge amount of experimental data from colliders, direct detection and indirect detection experiments, that must be taken into account to determine the viability of a model of physics beyond the Standard Model (BSM).

JCAP03(2016)019
In this regard, theoreticians are making a special effort to reduce the uncertainties from the theoretical side, for instance, in the calculation of the Higgs mass, which improves notably the interpretation of the experimental data in the context of different BSM scenarios.
Direct detection experimental collaborations adopt a set of assumptions about the DM nature and the DM halo which are useful in order to compare different results within a unified framework. More specifically, their results are traditionally presented considering only SI interactions 1 and the same interaction strength of DM with protons and neutrons, i.e. with a neutron to proton ratio f n /f p = 1. The speed distribution of DM particles is taken to follow a Maxwell-Boltzmann distribution, a consequence of assuming the halo to be an isothermal and isotropic sphere, the Standard Halo Model (SHM). Then, these results are used to constrain the parameter space of models, like supersymmetry (SUSY), in which theoretical calculations are performed with a high degree of accuracy. Furthermore, in many models the assumptions that experimental collaborations consider do not apply neither for the DM particle considered nor for the DM halo profile. 2 This is a consequence of the aforementioned over-simplifying assumptions, as realistic particle physics models of DM (e.g. SUSY) allow for different couplings to neutrons and protons and similar contributions from SI and SD interactions. Then, in principle, it is not fully consistent to employ the upper limits quoted by experimental collaborations to constrain general particle physics models of DM. However, to take into account all these subtleties is not an easy task, since it would require the simulation of a direct detection experiment, which makes very difficult their inclusion when scanning the parameter space of a model.
All in all, it is usual to compare the predictions of a given model with direct detection upper limits at face value, this is, considering the SI and SD interactions separately, 3 an equal contribution from protons and neutrons, i.e. f n /f p = 1, and assuming the SHM. Nonetheless, it has been pointed out that either the neutron to proton ratios [3][4][5][6][7][8] or the SD structure functions [9] and speed distributions [10][11][12][13][14][15][16] chosen can vary the interpretation of direct detection data substantially. In fact, to account for the latter effect, several halo independent methods have been proposed recently [17][18][19][20][21][22][23][24][25][26][27][28]. In this article, we present a set of tabulated data that contain the necessary information to calculate the expected number of DM induced events in the SuperCDMS and LUX experiments taking into account different possibilities of the mentioned ingredients. These data allow a fast check of the validity of a given model solution considering both, the SI and SD interactions, for arbitrary values of f n /f p and a n /a p , using realistic DM speed distributions and including the uncertainties related to the SD form factors. We provide the data and example codes at this site http://goo.gl/1CDFYi.
To this end, we have carried out a decomposition of the SuperCDMS and LUX detector signals in intensity (the part related to the cross sections) and shape (the part dependent on the energy) in such a way that the total number of predicted events in these experiments can be easily calculated even if the particle physics and/or DM profile assumptions made by collaborations do not allow a consistent and serious test. The tabulated data we have derived are very similar to the scaling functions introduced in ref. [29]. 4 Here, we extend this information to include more realistic descriptions of the DM halo profile, and different SD structure functions. These data can be used to evaluate at which confidence level a certain point of the parameter space of a model is excluded or allowed, without any computationally intensive calculation. Such an observable can then be translated into upper limits on the properties of the DM candidate. The aim behind this work is to provide the particle physics community with enough information to include direct detection bounds in terms of the DM expected number of events, in a complete and consistent manner.
To exemplify the impact that the combined limits (using the SI and SD interactions at the same time), together with distinct assumptions about the speed distribution and the SD structure functions, have on the allowed solutions in a model parameter space, we have chosen as case studies two well motivated and popular models: the Next-to-Minimal Supersymmetric Standard Model (NMSSM) and a Z ′ portal with Dirac DM. We show explicitly step by step how the previous considerations can affect the phenomenological viability of the solutions found. Moreover, the data provided in this work for different DM halo profiles, can be used to perform a consistent analysis of a model in light of DM direct and indirect detection experiments using the same DM halo density profile assumptions. This paper is organised as follows. In section 2, we present the basic expressions that allow to compute a DM signal in the LUX and SuperCDMS detectors. Then, in section 3 we derive the upper limits from these experiments under different assumptions in the velocity distribution and SD form factors, emphasizing their impact on the predicted number of events. In section 4 and 5, we show how the inclusion of the upper bounds in terms of the number of events, considering both the SI and SD interactions, realistic neutron to proton ratios and different speed distributions, has a notable impact on the allowed parameter space of the NMSSM and Z ′ portal models, respectively. Then, we repeat this exercise using different SD structure functions, in such a way that we can quantify the effect on the parameter space of these scenarios. Finally, the conclusions are left for section 6. In appendix A, we give detailed information about the extraction of the upper bounds and how to apply them to a generic model, and in appendix B we briefly describe the data files and example codes attached to this work.

Direct detection of dark matter
Let us start by summarising the basic expressions that describe the differential rate in a DM direct detection experiment [30] (for a recent review, see ref. [31]). The differential event rate for the elastic scattering of a WIMP with mass m χ off a nucleus with mass m N reads 5 where ρ 0 is the local DM density, f ( v) is the DM speed distribution normalised to unity and velocities are expressed in the detector reference frame. The integration over the DM velocity is computed from the minimum DM speed needed to induce a recoil of energy where µ N is the WIMP-nucleus reduced mass, to infinity. The DM-nucleus differential cross section dσ/dE R is computed from the Lagrangian that describes the interaction of a given WIMP with ordinary matter and encodes the microscopic interactions of DM particles and quarks inside nucleons. 5 Typically given in units of counts/day/kg/keV.

JCAP03(2016)019
As mentioned before, direct detection experiments usually assume the SHM, which describes the Galaxy as an isotropic isothermal sphere. This model provides a useful and common framework to compare different experiments, but definitely it is not a realistic description of our Galaxy. Due to this fact, we will use the mean speed distributions of ref. [38] derived separately for three possible choices of the DM density profile, namely a Navarro-Frenk-White (NFW) profile, an Einasto one and a Burkert profile. 6 These distributions are extracted from self-consistent DM phase space distributions that are in agreement with the astrophysical observables of the Galaxy and thereby, we expect them to give a more realistic description of the speed distribution of DM particles in the Solar neighbourhood. For completeness and to make it easy for the reader to compare with other results in the literature, including those provided by the experimental collaborations, we will use the SHM as well.
The DM-nucleus scattering cross section can be split into a SI and a SD contribution. The total cross section can be then calculated by adding these terms using nuclear wave functions. Thus, the differential cross section is given by where σ SI, N 0 and σ SD, N 0 are the SI and SD components of the so-called cross sections at zero momentum transfer, and the form factors F 2 SI, SD (E R ) account for the coherence loss which leads to a suppression in the event rate for heavy nuclei in the SI and SD contributions (see ref. [32] for a complete description of these prescriptions). Thanks to the Fermi's Golden Rule in the Born approximation, we can factorize out all the energy (momentum) dependence of the scattering of DM off a nucleus inside a form factor F 2 (q). This allows us to express σ SI, N 0 and σ SD, N 0 as [32,33] where µ p is the proton-WIMP reduced mass; S p and S n are the expectation values of the spin of protons and neutrons respectively, in their z-projection evaluated in the state with the maximal value of the angular momentum; f p , f n and a p , a n are the effective WIMP couplings to protons and neutrons in the SI and SD case, respectively. 7 The target material is defined by the following parameters, the atomic number Z, the mass number A, and the total nuclear spin J. In these expressions, we have defined the WIMP-proton cross sections (σ SI p ,σ SD p ) as, where G F is the Fermi coupling constant and the spin of the DM particles is considered equal to that of the proton. Note that eq. (2.3) reduces to these expressions when the quantities 6 The mean values of the escape velocity, Sun's velocity, and local DM density for each profile can be found in tables 2, 3 and 4 of ref. [38]. 7 Notice that for simplicity we have not included here a possible vector coupling (corresponding to non-Majorana DM particles), which would lead to an extra contribution in the expression for σ SI, N referred to the nucleus are substituted by those of the proton, i.e. using S p = 1/2, S n = 0, J = 1/2, A = Z = 1.

Differential rate in SuperCDMS
The actual energy observed by the SuperCDMS experiment is not the recoil energy presented in eq.(2.1), but the total phonon energy produced by a collision event. To take into account such a change, we follow the prescription given in ref. [34]. The phonon energy, E p , can be expressed as the sum of the recoil energy and the Neganov-Luke phonon energy, produced proportionally to the charge energy, where ∆V = 4V is the bias voltage applied to the detectors, e is the electron charge, and finally E Q is the mean charge for nuclear recoils. The latter is calibrated as a function of E p and then parametrised as follows, where the parameters α i depend on the corresponding detector and can be found in ref. [34]. Finally, to convert the differential rate terms of the nuclear recoil energy, given in eq.(2.1), to the differential rate in total phonon energy, we apply the following change of variables, (2.7) In figure 1, we show the differential rate in a Ge experiment as a function of the recoil energy (left panel) and the total phonon energy (right panel) for a SI cross section off nucleons σ SI p,n = 10 −6 pb, neglecting contributions from the SD interactions. In the right panel, we have used the charge model of the T2Z1 detector as it is given in ref. [34]. We have generated the differential rates for distinct DM masses and different speed distributions to show explicitly the impact that the velocity distribution choice has in the differential rate. Our SHM results (solid lines) are in very good agreement with that of the SuperCDMS Collaboration when comparing these results with figure 2 of ref. [34]. As it is shown, generally, the variation of the speed distribution tends to predict higher rates in the low mass region (all DM masses used in this figure would fall in such category) with respect to the SHM case. This would induce more stringent bounds as we will see later on.

Differential rate in LUX
We address now the LUX experiment. In general, the observed signal in liquid Xe-based detectors is the number of photoelectrons (PE) S1. To simulate it in LUX, we have followed the description given in ref. [35], where the number of expected PE, ν(E R ), reads   where Q γ = 0.14 is the photon detection efficiency for LUX [2]. L eff (E R ) is the absolute scintillation yield, which we have extracted from ref. [36]. The signal rate in number of photoelectrons n is then given by where Poiss (n|ν(E R )) is a Poisson distribution with expectation value ν(E R ), and the differential rate in E R is given by eq. (2.1). Taking into account the finite average singlephotoelectron resolution σ PMT = 0.37 PE of the photomultipliers [37], the resulting S1spectrum is given by Gauss(S1|n, where ζ(S1) is the acceptance of the applied cuts. Gauss(S1|n, √ nσ PMT ) denotes a normal distribution with mean n and standard deviation √ nσ PMT .
In figure 2, we show a set of differential rates for different DM masses and speed distributions, as a function of the recoil energy (left panel) and S1 (right panel) in the DM search energy window. The non-monotonic behaviour of the differential rates in the right panel is due to the efficiency of the LUX detector. As we can note, now the spectra are flatter than for Ge (see figure 1) as a consequence of the heavier Xe isotopes. The same tendency as before is observed, the predicted rates are higher for the self-consistent speed distributions than that predicted for the SHM. In this case, since the DM masses used here are heavier in some cases than that considered for SuperCDMS, we see that for a DM of 50 GeV the effect of the speed distribution starts to be negligible. This effect will also appear in the LUX upper limits.

Upper limits for SuperCDMS and LUX
To extract the upper limits, we use the Yellin's maximum gap method [39] at 90% C.L., which has been shown to produce good results for both LUX [40][41][42][43][44][45][46][47][48] and SuperCDMS [45][46][47][48]. This statistical tool allows to extract the upper bounds in presence of an unknown background in terms of the expected DM events in the signal region. To such end, we calculate the total number of expected DM events using the procedure outlined in appendix A for different speed distributions as well as SD structure functions. Notice that this method does not permit to combine the information from both experiments at the same time, however, there has been other proposed methods, like the ones presented in refs. [25,48], that allow such a combination. In the following, we provide the information about the SuperCDMS and LUX experimental setups used to calculate the tabulated data attached to this work, and therefore, the upper limits.
• SuperCDMS. The SuperCDMS Collaboration has reported a first search for DM using the latest background rejection capabilities for an exposure of 577 kg-days [1]. This analysis relies on a low-threshold performance of the detector, and hence, it was carried out only for DM masses below 30 GeV. Eleven events were observed after unblinding, allowing to set the most stringent upper bound on the SI DM-nucleon cross section for masses approximately below 6 GeV. For this analysis, only seven detectors (T1Z1, T2Z1, T2Z2, T4Z2, T4Z3, T5Z2 and T5Z3) had significant sensitivity to very low energy recoils, and then these were employed to set the upper limit. The remaining detectors were used to veto events that scatter in multiple detectors instead, since multiple scatter event topology has a negligible probability to occur for DM.
We calculate the expected signal for each detector separately using the efficiencies, exposures and charge models specified in ref.
[34] in terms of the total phonon energy given in eq. (2.7). We also include an estimation of the energy resolution, considering the 1.3 keV activation line at a total phonon energy of approximately 3 keV from figure 3 of ref. [1], to be energy independent and ∆E p = 0.25 keV. To find the upper limits using the maximum gap method we have taken into account the eleven observed events in the signal region E p =[2-12.5] keV.
• LUX. The LUX Collaboration reported null results in the search for DM particles for an exposure of 10065 kg days, which allowed them to place the most stringent upper limits on the SI interactions of DM off nucleons above 6 GeV approximately. The signal region was set to be 2-30 PE in S1. The acceptance of the detector is shown in the bottom of figure 1 of ref. [2] and we add an extra 1/2 factor to account for the 50% of nuclear recoil acceptance. We use the S1 single PE resolution σ P M T = 0.37 PE [37], a 14% of photon detection efficiency, and the absolute scintillation efficiency digitized from ref. [2].
In figure 3, we show the upper bounds on the SI proton cross section versus the DM mass for SuperCDMS and LUX experiments. In black (grey), we have plotted the official limit of the LUX [2] (SuperCDMS [1]) Collaboration, whereas in blue and green, respectively, we represent our findings of the upper limits at 90% C.L. using the Yellin's maximum gap method. For consistency, our results have been obtained using the same speed distributions as in the case of the collaborations. In the case of LUX, to extract the upper bounds we have considered zero candidate events observed, even though the results of the collaboration showed one candidate event marginally close to the background region in the log 10 (S2/S1) − JCAP03(2016)019 S1 plane. This fact is reflected as a slight disagreement between our results and those from the collaboration for masses below 30 GeV. On the contrary, for heavier masses, the agreement between our findings and the official limit is excellent. For SuperCDMS, the small discrepancies between our results and those of the collaboration are only due to the statistical method applied, since the latter used the Yellin's optimum interval method [39].

Limits on the SI cross section
Let us start with the results for the SI interactions of DM off nucleons. In figure 4, we show the upper limits for the LUX (blue) and SuperCDMS (green) experiments. All these curves have been derived under the usual assumptions, f n /f p = 1 and σ SD p,n = 0. Different type of curves correspond to different choices of the DM speed distribution. For the SHM (solid curve), we have used the same values of the speed distribution as in refs. [1,2] with a local DM density ρ 0 = 0.3 GeV/cm 3 . On the other hand, for the remaining profiles, our results have been obtained with the speed distributions from ref. [38] and are shown as dashed-dotted curves for the NFW profile, dotted for the Einasto profile, and dashed for the Burkert profile. In addition, we have rescaled all the upper bounds according to different values of the local DM density. We depict in darker colours the upper limits in which we have used the same value of the local density as in the SHM case, in order to show how the shape of the DM velocity distribution impacts on the upper bounds. In contrast, curves with lighter colours correspond to upper limits derived using the mean values that the authors of ref. [38] obtained for each of the above mentioned profiles, namely, ρ NFW 0 = 0.41, ρ Ein 0 = 0.42 and ρ Bur 0 = 0.41 GeV/cm 3 . As expected, the impact of the assumed DM density profile on the upper bounds is specially notable in the low mass region, where small changes in the tail of the speed distribution function can enhance or decrease prominently the number of expected DM events above the threshold [11][12][13][14][15][16]. Therefore, in the m χ < 30 GeV region, our upper limits using well-motivated DM profiles and a self-consistent extraction of the speed distribution tend to be more stringent than using the SHM, as expected from figures 1 and 2. These differences can be as high as one order of magnitude for masses below 10 GeV. Notice that the mean values of v esc we have used are in the three cases below 544 km/s, the value used for the SHM. This means that the difference in the upper limits cannot be ascribed to the different escape velocities assumed. In fact, this shows that the slope of the speed distribution at high velocities is smaller in the case of the NFW, Einasto and Burkert profiles. Needless to say that the values used for these distributions are subject to important uncertainties which actually in some cases would push these upper bounds closer to the SHM result.

Limits on the SD cross section
Now, we move to the case of SD interactions. Our results for the upper limits on the SD cross section of DM off protons (left panel) and neutrons (right panel) for the LUX and SuperCDMS experiments are displayed in figure 5. We show the SHM (solid) and NFW (dashed-dotted curves) cases, where we have included different SD structure functions. Namely, we have employed the results obtained with chiral effective field theory (EFT) currents [49] for Ge and Xe isotopes (thicker lines), we have also used the Bonn A potential [50] for Xe and the Dimitrov et al. calculation [51] for Ge (medium lines), and finally, the Nijmegen potential [50] for Xe and the Ressell et al. calculation [52] for Ge (thinner lines). This kind of uncertainties on direct detection experimental results has been shown to be important in ref. [9], which motivates us to include the upper bounds for the different calculations of the SD structure functions in the literature. We have checked that these results are in good concordance with those of ref. [53] which were performed using a profile likelihood analysis.  Figure 6. LUX combined SI and SD upper limits using a NFW profile for different masses, and a n /a p = 1 (left) and a n /a p = −1 (right). The colour convention is as in figure 2. Horizontal and vertical lines represent the limits for SI and SD interactions only, respectively.
We note from figure 5 that due to the nuclear spin of both Ge and Xe, the constraints from the neutron component of the cross section is substantially more stringent than that of protons. Nonetheless, the difference between distinct calculations of the SD structure functions is smaller in the neutron case. These two facts together imply that the exclusion of models from LUX and SuperCDMS results is not going to undergo strong variations from the choice of the SD structure functions. By contrast, from the left panel of figure 5 one can see that for both, LUX and SuperCDMS, the change in the σ SD χ−p owing to different SD structure functions can be of the order 10. Of course, this might be particularly important for models in which the DM particle couples to protons rather than to neutrons.

Combined limits for SI and SD interactions
In this section, we derive the LUX limits for a specific choice of the DM mass, and then we calculate the upper bounds on σ SI χp , σ SD χn and the combined limits. We have assumed f n /f p = 1 and a n /a p = ±1. We perform this exercise to evidence the difference of taking into account both cross sections at the same time instead of considering each one separately.
In figure 6, we show the LUX combined bounds for a NFW profile, using ρ 0 = 0.3 GeV/cm 3 , and for the aforementioned neutron to proton ratios for the SI and SD interactions. The structure functions used in the SD case are those from the chiral EFT [49]. Besides, we have fixed the DM mass to: 10 GeV (red), 50 GeV (green) and 100 GeV (magenta curves). Horizontal lines are the SI limits considering zero contributions from the SD interactions, while vertical lines represent the opposite, i.e. SD limits with no contribution from the SI interactions. The rectangles delimited by these two kind of lines represent the cases in which the limits are considered for each interaction separately. Conversely, when one takes into account both interactions at the same time, each rectangle bends and a new excluded-to-be region enclosed by the rectangles and the new limits, labeled by SI + SD, arises. Of course, a source of variation of these combined limits is given by the change of -11 -JCAP03(2016)019 the ratios f n /f p for SI interactions, and a n /a p for the SD ones. For the latter, in figure 6 we have plotted two different values, a n /a p = 1 in the left panel, and a n /a p = −1 in the right panel. The difference among these two results highlights the importance that these quantities have in the upper bounds. Regarding the SI interactions, the effect can be also remarkable, specially when f n /f p gets close to −0.7, i.e. the Xe-phobic region, but also for |f n /f p | ≫ 1.
Indeed, the aforesaid excluded regions are not taken into account if one applies the LUX results for the SI and SD interactions separately. Obviously, to extract information to rule out these regions it is necessary to compute the combined upper limits, but since the SI contribution is normally dominant, experimental collaborations only show results for one cross-section at a time. However, our results pose a question about the importance of the inclusion of such information to probe consistently DM models. As already stated, usually the SI contribution to direct detection experiments is dominant over that of the SD, because the former is summed coherently for all nucleons, hence scaling as A 2 (see eq. (2.3)). Nevertheless, there are scenarios in which the SD contribution can be of the same order than that of the SI, or even dominate over it [54]. In the following, we will show that proper combined limits are an important ingredient to take into account when analysing models. To such end, we will study two well motivated and popular DM models, the NMSSM and a Z ′ portal.

Next-to-Minimal Supersymmetric Standard Model
The NMSSM is one of the most popular supersymmetric models (see for instance ref. [55] for a review and ref. [56] for a DM related analysis). In this model, the DM candidate is the lightest neutralino (χ 0 1 ) which is, in general, a mixture of the superpartners of the gauge bosons and the Standard Model Higgs, plus a singlet scalar which is added in order to solve the µ-problem. Since neutralinos are Majorana particles they contribute to both, SI and SD interactions, and thereby are perfect candidates to search for the aforementioned cases.
To look for phenomenologically viable solutions of this model, we have carried out a series of scans over the parameter space with MultiNest 3.9 [71][72][73]. To this aim, we have defined a likelihood function whose inputs are the neutralino relic density, the SM Higgs mass, BR(B s → µ + µ − ), and BR(b → sγ), which are taken as gaussian probability distribution functions around the measured values with 2σ deviations. The gluino soft mass is fixed to M 3 = 1500 GeV, the trilinear parameters to A U = 3700 GeV, A D = 2000 GeV, and A E = −1000 GeV, and the soft scalar masses of sleptons and squarks to mL i = mẼ i = 300 GeV and mQ i = mŨ i = mD i = 1500 GeV, respectively, where the index i runs over the three families. This conservative choice of the squark masses is motivated by the LHC null searches. Also note that despite the high trilinear term A U , the instability against charge-and/or colorbreaking minima is avoided since the squark soft masses are at the TeV scale [74].
We have performed three scans over the parameter space where the different input parameters are varied according to table 1. The neutralino relic abundance has been calculated with micrOMEGAs 3.6.9 [75]. We have used NMSSMTools 4.1.2 [76][77][78] to compute the NMSSM mass spectrum, the masses of the Higgs bosons including full two-loop contributions, and the relevant low-energy phenomenology observables. LHC measurements  of the Higgs properties are included by constructing the ∆χ 2 distribution from ref. [79] and allowing 3σ deviations. Limits on the velocity averaged cross section of DM particles from the gamma-ray null observation of dwarf spheroidal galaxies (dSphs) are also taken into account, following the procedure sketched in ref. [80] for the Pass 8 data [81].

Combined limits for different speed distributions
In this section, we will discuss how the inclusion of direct detection upper limits taking into account the SI and SD interactions at the same time can have an impact on the viability of DM scenarios in the NMSSM that otherwise will be allowed. More precisely, we will use the results from previous sections to illustrate how different ways of including direct detection bounds can affect the parameter space of the NMSSM. We will also show the impact of JCAP03(2016)019  the astrophysical uncertainties from the speed distributions, leaving to the next section the analysis of the effect of these uncertainties on the SD structure functions.
To show schematically how different solutions are excluded depending on the assumptions considered when evaluating the limits, we will use the following procedure: (1) SI: we will apply first the SuperCDMS and LUX bounds only for the SI interactions and f / f p = 1 to the viable solutions found, using the SHM. This bound correspond to the usual bound published by the SuperCDMS and LUX Collaborations.
(2) SI/SD: then, we will add to the previous limits the bounds that only take into account the SD interactions, also for the SHM. This step is equivalent to apply the combined limits generating the rectangles shown in figure 6.
(3) SI + SD: next, instead of the previous considerations we will calculate the bounds taking into account all the contributions, SI and SD and the specific neutron to proton ratios. This is also done for the SHM.
Notice that with the purpose of isolating the effect that speed distributions have on the exclusion regions, in none of the previous steps we have varied the value of the local density, ρ 0 . In fact, if we change the local density value in agreement with those extracted in ref. [38], according to figure 4, the upper bounds become more stringent.
The set of plots shown in figure 7 depicts how the previous steps affect the exclusion regions of neutralinos in the NMSSM. All the plots show the points excluded by the Super-CDMS and LUX results in the plane ξσ SĨ In the first row left-hand side plot, we display the solutions excluded when the SI proton cross section weighted by ξ is above the nominal SuperCDMS and LUX limits. This corresponds to step (1), as labelled in the plot. In the same row, the right-hand side plot, includes also the upper bounds both for SD proton and neutron interactions. As it can be seen, when compared with respect to the left-hand side plot, a new population of excluded points has appeared around mχ0 1 ≈ 30 GeV. These points correspond to solutions in which the SD neutron cross section is above the LUX bound. 9 Next, we move to the middle row of figure 7. In the left panel, we note that many solutions previously allowed by LUX are now ruled out. These points fall precisely in the regions labelled by SI + SD in figure 6. All the points that have appeared at this step have a similar contribution from the SI and SD interactions to the expected number of events in LUX. Therefore, all these solutions require a careful calculation of the LUX upper bound to be correctly excluded. This nicely highlight the importance of using combined limits and showcases how the upper limits should be compared point per point in the parameter space. This task can be achieved very quickly with the help of the tables provided in this work.
The remaining plots correspond to step (4). As previously discussed, the impact of considering different velocity distributions is notable for neutralino masses below approximately 50 GeV, as expected from the results shown in figure 4. This light mass region is actually one of the favoured regions that explain the Fermi-GeV excess in terms of annihilating DM particles [91][92][93][94][95][96][97][98][99][100]. Moreover, this excess prefers a NFW-like profile; hence, the use of realistic speed distributions according to the Galactic DM profile is required to suitably constraint possible DM explanations of this excess with direct detection upper bounds [103]. 9 Remind that Xe-based experiments are much more sensitive to the neutron contribution of the SD interactions since the total spin of these nuclei are dominated by neutrons.

JCAP03(2016)019
Thus, our findings point out that in order to explore robustly the parameter space of the NMSSM, one has to include properly the SuperCDMS and LUX limits. It is worth noting that steps (1) and (2) assume a fixed neutron to proton ratio for both interactions, which does not coincide exactly with those of neutralinos, when we apply step (3) this subtlety is considered using the exact values that neutralinos have for these quantities. Therefore, these ratios can affect the exclusion of some of the previously allowed points.

Combined limits for different SD structure functions
As seen before, many solutions found in our scans predict a SD contribution similar to that of the SI. In section 3.2, we have shown that the use of distinct structure functions lead to differences in the LUX upper bounds that can be significant. Therefore, we must expect that the choice of the structure functions is going to impact on the excluded regions of neutralinos in the NMSSM.
In figure 8, we depict the excluded points assuming different SD structure functions (columns) and speed distributions (rows). For the sake of completeness, we have plotted in light blue those points in which the contribution to the DM expected events in LUX from the SD interactions is dominant. As we can observe, in figure 8 from left to right panels, light blue points are generally affected by the choice of the structure functions and, in agreement with figure 5, the chiral EFT structure function is able to rule out smaller regions of the parameter space. It is noteworthy that the BonnA potential model of the SD structure functions provides the most stringent limits of Xe-based experiments, highlighted as a higher population of light blue points in the middle column of figure 8 with respect to the other two (compare also steps (2) and (3) in figure 5). Notice that the SD dominant solutions we have found accumulate around neutralino masses of 30 GeV, while those with similar SI and SD contributions in LUX range from 20 to 100 GeV approximately (see also in figure 7 the change from step (2) to step (3)). This means that the choice of the SD structure function has a notable effect on the allowed (excluded) solutions for NMSSM neutralinos below 100 GeV, which is precisely the region where most of the explanations of the Fermi-GeV excess in this model lie [66,[68][69][70]104].

Z ′ portal
We now turn our attention to another popular and well-motivated extension of the SM, the Z ′ portal. In this scenario, a hidden sector communicates with the SM via a gauge boson, provided that the SM is enlarged with an extra abelian gauge group [105]. The phenomenology of these constructions is very rich, and ranges from colliders to direct and indirect searches for DM [106][107][108][109][110][111][112][113][114][115][116][117][118].
We will consider the same construction as in ref. [116], which is based on a certain type II string compactifications with intersecting branes. The key point of these scenarios is that gauge bosons acquire a mass through the Stückelberg mechanism [119,120], except of course the one corresponding to the hypercharge, which remains massless in the phase of unbroken electroweak symmetry.
We will assume throughout this section that the DM in this scenario is a Dirac fermion, ψ, which lives in the hidden sector. The phenomenology of this kind of DM will be driven by the interaction with the SM fermions through the exchange of a Z ′ boson. The couplings of the Z ′ to the SM fermions will be determined by specific combinations of four parameters, a, b, c and d [116]. We have performed a scan over the parameter space, with the parameter ranges given in table 2, where h is the coupling of DM to the Z ′ boson. To this end, we have -17 -

JCAP03(2016)019
Parameter Range implemented the model in micrOMEGAs 3.6.9 [75]. We have taken into account that there exist certain combinations of the these parameters that are not allowed by the orthogonality and normalization conditions imposed on the Z ′ mass eigenvectors [116]. In general, U(1) extensions of the SM can be constrained using resonances at colliders, since the coupling of the Z ′ boson to leptons and quarks contributes to the appearance of dimuon, dielectron and dijet resonances. Here, we have applied the ATLAS results for high mass resonances decaying into a µ + µ − or an e + e − pair at a center of mass energy √ s = 8 TeV and luminosities of 20.5 fb −1 and 20.3 fb −1 for dimuons and dielectrons resonances, respectively [121]. Furthermore, we have used the searches for dijet resonances and monojets plus missing energy that receive additional contributions from the presence of a Z ′ boson, both at the LHC and Tevatron colliders [122][123][124]. In order to implement these constraints, we have followed the procedure of ref. [116]. It is worth mentioning that we do not impose any bound on the relic abundance of ψ. This is a consequence of the potential complexity of the gauge structure and matter content of the hidden sector, which might induce other contributions that account for the correct abundance such annihilations of ψ into hidden sector matter. Of course, this will affect the predictions of the model for indirect dark matter searches. To deal with this, we have defined ξ as in the NMSSM, but in this case if Ω ψ h 2 > 0.13 then we set ξ = 1. 10 Finally, we have imposed the upper bounds on the annihilation cross section of ψ from dSph galaxies using the Pass 8 data of the Fermi-LAT satellite [81]. These limits are presented for specific SM final states, while in this model each point of the parameter space entails DM annihilation into a combination of different SM particles. This circumstance prevents us to use these Fermi-LAT results directly, but we have to weight each upper limit by the corresponding annihilation percentage instead.

Combined limits for different speed distributions
As in the NMSSM case, presented in section 4, in this section we will evaluate the impact of the speed distribution on the experimentally allowed solutions of the Z ′ portal with Dirac DM. We will follow the same steps described before for the NMSSM.
In figure 9, we show the results for the Z ′ portal model using the same steps as in figure 7. Notably, from step (2) to (3) many new excluded solutions appear below the LUX SHM limit (blue curve). Unlike in the NMSSM, all these ruled out points correspond to JCAP03(2016)019  solutions in which the ratio f n /f p is different from 1, and in many of these cases, the neutron contribution to the SI cross section is much higher than that of the proton. For this reason, many points that lie far below the LUX bound on the SI-proton cross section are excluded. In addition, some of the solutions that were excluded in step (2) are now allowed. These solutions correspond to those with f n /f p ≈ −0.7, the region in which the sensitivity of Xebased experiments substantially decreases, the so-called Xe-phobic region. Also, solutions with |f n /f p | < 1 are now allowed, since the neutron contribution does not reach the size of the proton one, and hence, the upper bound weakens.
Regarding the impact of the velocity distribution, one can see in figure 9 that from step (3) on, more points are excluded in the low mass region. Although in this case, the difference is less visible than in the NMSSM, because the density of points for masses below 30 GeV is not very high. Furthermore, the fact that in this model f n /f p can virtually take any value, makes more difficult to identify points where the effect of the speed distribution is remarkable.

Combined limits for different SD structure functions
In figure 10, we show the results for the Z ′ portal model using the same speed distributions and SD structure functions as in figure 8. Notice that the absence of light blue points (like those appearing in figure 8 for neutralino DM) means that in this case we have not found any point of the parameter space in which the SD contribution to the total number of expected events in LUX dominates. Consequently, we observe that the change in the SD structure functions employed to impose the limits does not affect the exclusion of the solutions found.

Conclusions
In this paper, we have performed a thorough study of the impact of the assumptions usually made to compute DM direct detection upper bounds. We have used the NMSSM and a Z ′ portal DM models to illustrate how the use of direct detection limits, including SI and SD interactions at the same time, different neutron to proton ratios, distinct DM density profiles as well as different SD structure functions, has important implications on the excluded parameter spaces.
We have extracted the LUX and SuperCDMS upper bounds under many different assumptions about the speed distribution and SD structure functions. For the former, we have used the velocity distributions derived in ref. [38] for three different choices of the DM density profile. In ref. [38], the authors computed the speed distribution in a self-consistent way from the gravitational potential of the Milky Way. We have employed these distributions, to compute upper limits that rely on a realistic description of the DM halo, and hence, produce results consistent with a wide range of astronomical observations. As expected, the differences in the upper bounds from the speed distribution focus mainly on the light mass range, m DM 50 GeV. Another source of uncertainty for direct detection experiments are the SD structure functions. By using the different results in the literature, we have explicitly shown the variation in the upper limits due to the choice of the SD structure functions.
In order to incorporate all these ingredients when analysing models, in each of the mentioned cases, we have carried out a decomposition of the expected DM signal in both experiments, LUX and SuperCDMS, in terms of several expressions that can be used to derive upper limits using statistical tools based on the total number of expected events. In this paper, we have used the maximum gap method, and obtained tabulated data, that is suitable for this method, for each of the different combinations of these ingredients. This  -21 -data set allows to constrain any generic model with both SI and SD interactions, for protons and neutrons. The data and python scripts that exemplify their use can be downloaded from this site http://goo.gl/1CDFYi.
As case studies to apply our findings, we have chosen the NMSSM and the Z ′ portal with Dirac DM. We have performed a series of scans over the parameter space of both models imposing the latest experimental limits. More specifically, we have applied to the resulting points the constraints from the LHC and Tevatron, Planck, low energy observables and Fermi-LAT bounds from dSph galaxies. Then, we have imposed how the inclusion of the upper limits from LUX and SuperCDMS taking into account the SI and SD interactions, and different assumptions about the halo profile and SD structure functions, to analyse the effect of each of these factors on the allowed/excluded regions of the parameter space. Our findings point out that a more careful implementation of direct detection limits is needed when analysing models in light of the current experimental constraints. More specifically, in the NMSSM, the nature of the solutions found requires a combined implementation of the SI and SD bounds at the same time, since many solutions have similar contributions from both types of interaction. For this reason, in this model, not only the halo profile choice, but also the SD structure functions used to derive the SD contribution play an important role. In the case of the Z ′ portal, the necessity of applying combined limits arises from the neutron and proton contributions to the SI interactions. This model, with a Dirac DM candidate, shows strong deviations from the f n = f p equality. This situation is translated into huge variations of direct detection upper limits that should be considered for a careful analysis of the parameter space.
The amount of experimental data nowadays is enormous. Huge experimental and theoretical efforts are being made to interpret these data in terms of many different models. With this purpose, theoreticians are performing very fine calculations of observables that can be constrained in light of the current data. In this sense, our results provide a missing piece to this effort, making easier the implementation of direct DM detection limits using more realistic assumptions; namely the SI and SD contributions at the same time, speed distributions in agreement with astrophysical observations, and different SD structure functions.
Note added. After this work was completed a new bound on the SI interactions from LUX collaboration was publicly released [125]. This new analysis uses the 2013 data set but incorporates several advances such as calibration and background models, which results in a sensible improvement of the upper limit specially at very low masses. Since this work shows the importance of using different speed distributions and SD structure functions to derive upper limits, and the impact of these on different models, the conclusions of this work does not depend on the upper limit used, and the new LUX bound does not affect to our analysis. However, we will update the tabulated data attached to this work with the new LUX upper limit as soon as possible.

JCAP03(2016)019
A Extraction of upper limits using the maximum gap method In this appendix, we give detailed information about the extraction of the direct detection bounds for SuperCDMS and LUX using the data provided in this work. Notice that this data set can be used with two different, nevertheless equivalent, purposes. Its first application is to calculate the upper limits under a number of assumptions, in the same way the collaborations present their results for the SI proton cross section as a function of the mass (see for instance figure 4). The second application is to test if a given point in the parameter space of a model is allowed or excluded, as we have done to obtain figures 7-10.
As stated in the Introduction, by means of these tabulated data one can determine the expected number of DM events in SuperCDMS and LUX in an easy and quick way. Then, with these number of events, we can estimate the upper bounds on DM interactions with protons and/or neutrons, using the Yellin's maximum gap method [39], which produces excellent results (see figure 3 for a comparison between our results and those of the experimental collaborations). With the attached data, one can calculate the expected number of events in specific energy windows and the total number expected for a given experiment. The former and the latter are necessary to apply this method to SuperCDMS (for LUX, we assume zero observed candidate events and then we will only need the total number of expected events).
Let us start considering a generic direct detection experiment, denoted by Z. It is straightforward to show that the number of expected DM events for the SI interactions can be written as, where the F p,n,pn Z functions enclose all the information about the DM halo, form factors, energy resolution, as well as the transformation from the actual measured energy to recoil energy (see section 2). These SI functions have been normalized to a reference cross section of 10 −8 pb. Note that these functions have been calculated for a specific energy window, and then, are not valid for arbitrary energy intervals. Furthermore, we have reabsorbed all the constants, different for each function, namely Z 2 , (A − Z) 2 and Z(A − Z), and the exposure of the detector in the definition of the functions.
For the SD interactions, one can also perform the same exercise yielding, N SD = σ SD p 1 + a n a p 2 + 2 a n a p F 00 Z (m χ ) + 1 + a n a p 2 − 2 a n a p F 11 Z (m χ ) where the superindices of the F 00,01,11 Z functions correspond to the isoscalar and isovector decomposition of the SD structure functions. These SD functions have been normalized to a reference cross section of 10 −4 pb. Therefore, the total number of events in a generic DM experiment can be calculated just by summing the SI and SD contributions. Up to some constants, these functions are the same as those presented in ref. [29] for O 1 and O 4 in the framework of the effective field theory for DM interactions.
The maximum gap method finds the most constraining gap between observed events for setting an upper limit. This is, the energy range between two observed events with the -23 -

JCAP03(2016)019
highest number of events predicted by a model. Then, to calculate the bounds we determine the probability of this gap being smaller than a particular number of events, x [39], where µ is the total number of events in the whole search energy window, this is N SI + N SD from eqs. (A.1) and (A.2), and m is the greatest integer ≤ µ/x. First, note that if there is any observed candidate event for the experiment considered in the search window then x < µ, and the functions introduced in eqs. (A.1) and (A.2) to calculate x and µ are not the same, because the energy ranges used in each case are different. Second, for zero observed events µ = x by definition of the maximum gap [39], and the previous equation reduces to, In the following, we will explain how to use the above mentioned expressions to find the upper limits and to test points in the parameter space of a given model, for the specific cases of SuperCDMS and LUX.
• SuperCDMS bounds. Since SuperCDMS observed eleven candidate events in the DM search region [1], we have to define the functions for different energy ranges. All the functions have been extracted accordingly to the experimental setup of ref. [1] (see section 3 for more details) and following the prescription given in section 2. To evaluate the limits up to DM masses of 30 GeV, two different mass ranges must be considered m χ < 5.4 GeV and m χ ≥ 5.4 GeV. For each range, there is a different maximum gap and consequently we have to compute a distinct x, denoted as x 1 and x 2 . Now, to calculate µ and x 1,2 (x 1 or x 2 , depending on the DM mass), one must use eqs. (A.1) and (A.2) with the proper tabulated functions given in the provided datafiles (see appendix B).
To test if a specific particle physics input characterized by σ SI p , f n /f p , σ SD p , a n /a p and m χ is allowed or excluded at e.g. 90% C.L. by SuperCDMS, we compute µ and x 1,2 and then evaluate C 0 from eq. (A.3). If C 0 ≥ 0.9, the input is allowed, otherwise is excluded. On the other side, if we want to determine the upper limit at 90% C.L. for the SI interactions and f n /f p = 1, one has to calculate µ and x 1,2 using only the SI functions and choosing f n /f p = 1. Afterwards, for a specific DM mass, one finds the σ SI p such that C 0 = 0.90.
• LUX bounds. The calculation of LUX bounds is easier because we have considered zero candidate events. In this case, the calculation of C 0 only requires the computation of the total number of expected events, µ, for the whole energy range and valid for all DM masses. Hence, the suitable expression for C 0 is eq. (A.4). To derive the necessary functions, we have used the experimental setup described in ref. [2] (see section 3 for more details) and followed the prescriptions given in section 2. To perform the same tasks as for SuperCDMS, we follow the same procedure, the only difference being the definition of C 0 .
For LUX, we also give results for DM masses above 200 GeV, for different SD structure functions, and using only the SHM, since as we have seen, different speed distributions do not impact on the results for masses above 50 GeV approximately. These files contain an extra _extended in the filename, but they must be used in the same way as the other files.
In the header of each file, the distribution of the columns is given. E.g. in the file SCDMS_SI_SHM.txt, we find the tabulated functions needed to calculate µ, x 1 and x 2 (denoted by the suffixes, _Total, _MG1 and _MG2, respectively) using eq. (A.1), taking into account only SI interactions and for the SHM. Finally, we also provide three example codes for LUX and SuperCDMS, namely LUXexample.py, SCDMSexample.py, and LUXexample2.py. These python scripts implement the expressions of appendix A and illustrate how to use them together with the attached datafiles. The first two codes evaluate if a point of the parameter space, whose specific characteristics are included in the scripts as an example, is allowed or excluded at 90% C.L. for each experiment separately. The third code reads 1000 inputs from the NMSSM parameter space (not all constraints have been imposed to this data), saved in testdata.dat, checks if they are excluded or not by LUX, and then plots the allowed and excluded points in the σ SI p−χ -m χ plane. This is useful when one wants to compare a scan over the parameter of a given model with the combined LUX bounds.