Global fits of the dark matter-nucleon effective interactions

The effective theory of isoscalar dark matter-nucleon interactions mediated by heavy spin-one or spin-zero particles depends on 10 coupling constants besides the dark matter particle mass. Here we compare this 11-dimensional effective theory to current observations in a comprehensive statistical analysis of several direct detection experiments, including the recent LUX, SuperCDMS and CDMSlite results. From a multidimensional scan with about 3 million likelihood evaluations, we extract the marginalized posterior probability density functions (a Bayesian approach) and the profile likelihoods (a frequentist approach), as well as the associated credible regions and confidence levels, for each coupling constant vs dark matter mass and for each pair of coupling constants. We compare the Bayesian and frequentist approach in the light of the currently limited amount of data. We find that current direct detection data contain sufficient information to simultaneously constrain not only the familiar spin-independent and spin-dependent interactions, but also the remaining velocity and momentum dependent couplings predicted by the dark matter-nucleon effective theory. For current experiments associated with a null result, we find strong correlations between some pairs of coupling constants. For experiments that claim a signal (i.e., CoGeNT and DAMA), we find that pairs of coupling constants produce degenerate results.


Introduction
Only about one sixth of the total matter in the observable Universe is made of known particles [1]. The remaining part is an invisible and unidentified cosmic component called the dark matter [2][3][4][5]. The particles forming the dark matter component have up to now escaped detection. Astronomical observations and numerical simulations [6] show that dark matter clusters in large astrophysical structures (dark matter halos). The clustering of dark matter particles has inspired a number of complementary methods to detect them (see, e.g., [7] for a recent review). In particular, dark matter particles in our own Milky Way galaxy can be searched for using the direct detection technique [8], which has recently played an important role in this context. Several direct detection experiments have indeed reached the sensitivity to probe the dark matter paradigm (see, e.g., [9] for a review). The goal of the direct detection technique is to measure the energy deposited in an underground detector by Milky Way dark matter particles scattering on a target material [10]. This detection strategy is therefore an ideal tool to probe the foundations of the dark matter-nucleon interaction.
There is an extensive literature devoted to the study of dark matter scattering on nuclei in an underground detector (for an overview of this subject, see for instance Refs. [11][12][13][14][15][16] and references therein). The vast majority of these analyses rely on the assumption that dark matter interacts with the detector nuclei either through the nuclear charge density operator or through the nuclear spin-current density operator. The former is commonly called "spin-independent interaction," the latter "spin-dependent interaction," although strictly speaking any other operator is either spin-dependent or spin-independent. This approach to the dark matter direct detection is motivated by its simplicity and by the fact that these two interaction operators are naturally generated in the most popular dark matter models. Important examples are those based on the Minimal Supersymmetric Standard Model or on many of its extensions [17][18][19]. On the other end, there is no empirical evidence supporting this assumption, and Nature might actually be more complex, allowing a broader spectrum of possible dark matter-nucleon interactions.
In the past few years alternative types of dark matter-nucleon interactions have been proposed and their exploration is now undergoing a very active and productive phase. Phenomenologically attractive extensions of the standard paradigm involve velocity and momentum dependent interactions [20], isospin violating couplings [21], and new long-range interactions [22]. The study of these theoretical frameworks is still in progress. In this context, interesting results have for instance been found in studying anapole and magnetic dipole dark matter [23], light dark matter candidates [24], dark matter capture by the Sun [25], and benchmark models designed for dark matter searches at the LHC [26]. Momentum dependent interactions have also been explored in the context of extracting the phase-space distribution of dark matter particles with direct detection experiments [27]. Halo-independent analysis based on velocity and momentum dependent operators can be found in Ref. [28].
Recently, Refs. [29][30][31] proposed the idea of studying the dark matter-nucleon interaction with a non-relativistic effective theory approach similar to the one used in the 60's for exploring weak-interactions. Ref. [29] extends the work of Ref. [32] to a systematic and complete classification of dark matter-nucleon non-relativistic interactions under Galilean transformations and conservation of energy and momentum. A general method to translate experimental limits into constraints on the dark matter-nucleon couplings in the non-relativistic effective theory has been devised [33]. Publicly available Mathematica packages to perform these calculations are also available [31,33]. Exclusion limits for velocity and momentum dependent interaction operators have been obtained from single experiments in [30,34,35].
So far, the analysis of the non-relativistic effective theory has (1) considered the different interaction operators separately and (2) analyzed distinct direct detection experiments independently. A global analysis of the full multidimensional parameter space defining the effective theory of the dark matter-nucleon interaction is still missing. To tackle this challenge is the main aim of this work. Here we present the first comprehensive analysis of the multidimensional parameter space of the dark matter-nucleon effective theory. In this study all the couplings and the dark matter mass are simultaneously considered as free parameters. In addition, we have combined in a single analysis many different direct detection data, including the recent LUX, SuperCDMS and CDMSlite results. To achieve these goals we have exploited state-of-the-art Bayesian/frequentist numerical tools to sample the posterior probability density function and the profile likelihood of the model parameters. Importantly, we find that present direct detection data contain sufficient information to simultaneously con-strain all the interaction operators present in the effective theory of the dark matter-nucleon interaction.
The paper is organized as follows. In Sec. 2 we briefly review the non-relativistic effective theory of the dark matter direct detection proposed in Ref. [29]. Sec. 3 describes the statistical methods used in our analysis, whereas the datasets to which they are applied are introduced in Sec. 4. Secs. 5 and 6 are devoted to the presentation of the results, and Sec. 7 contains our conclusions. Appendix A describes the dependence of our results on the astrophysical assumptions, whereas Appendix B contains a list of the dark matter response functions relevant for this paper.

Effective theory of the dark matter-nucleon interaction
In this section we review the basic concepts and equations defining the effective theory of the dark matter-nucleon interaction. For a more detailed introduction to this subject we refer the reader to the original literature [29][30][31][32].
From the point of view of relativistic quantum field theory, effective dark matter-nucleon interactions can be constructed from Lorentz-invariant combinations of dark matter and nucleonic bilinear operators. In the dark matter-nucleon non-relativistic effective theory the interactions are restricted by Galilean invariance, energy and momentum conservation, and hermiticity [29]. These requirements allow to construct a generating set of five non-relativistic operators for the algebra of χ-nucleon effective interaction operators (here χ denotes the dark matter particle): the identity 1 χ 1 N , the momentum transfer 1 q, the χ-nucleon transverse relative velocity operator v ⊥ χN (with matrix element equal to v χN − q/2µ N , where v χN is the initial χ-nucleon relative velocity and µ N is the χ-nucleon reduced mass), and the dark matter and nucleon spin operators S χ 1 N and 1 χ S N , respectively. The most general effective theory at the dark matter-nucleon level involves products of the five generating operators. In this paper, we restrict ourselves to the exchange of a heavy spin-0 or spin-1 particle, and following Ref. [29], we limit ourselves to the 10 operators listed in Tab. 1 are difficult to generate in explicit particle models. For spin-0 dark matter particles, the spin operator S χ is identically zero. For spin-1/2 particles it is equal to σ/2, where σ i , i = 1, 2, 3, are the Pauli sigma matrices acting on the χ-spinor. For spin-1 dark matter particles the components of S χ are spin-1 representations of the angular momentum generators acting on the χ-vector.
The most general effective Hamiltonian describing the dark matter interaction with a point-like nucleon is then given by the following linear combination of operators Here τ 3 is the third isospin Pauli matrix. The coupling constants c τ i (τ = 0, 1) have dimension (mass) −2 , and are analogous to the Fermi constant G F . 2 The constants c 0 i correspond to isoscalar dark matter-nucleon interactions, whereas the constants c 1 i describe the isovector 1 Our definition of momentum transfer q is the common one in dark matter direct detection studies, namely q = pχ − p χ , where pχ and p χ are the initial and final dark matter momenta. Ref. [29] defines q with opposite sign. This explains the minus signs in the q-dependent operators in Tab. 1 and in the expression for v ⊥ χN . 2 We define the c τ i constants following Ref. [31]. Other definitions exist in the literature. For example, Ref. [33] has c τ 1 = 4mχmN c τ 1 . Table 1. List of the 10 non-relativistic operators defining the effective theory of the dark matternucleon interaction studied in this paper. The operators O i are the same as in Ref. [31].
interactions. Equivalently, c p i = (c 0 i + c 1 i )/2 and c n i = (c 0 i − c 1 i )/2 are the coupling constants for protons and neutrons, respectively. In this paper we restrict our analysis to isoscalar interactions (often but improperly called "isospin-conserving" interactions), i.e., we set c 1 i = 0. The interaction Hamiltonian used to calculate the cross section for dark matter scattering on nucleons bound in a detector nucleus is obtained from Eq. (2.1) by replacing the point-like charge and spin operators with the corresponding extended nuclear charge and spin-current densities, as for instance in Eq. 27 of Ref. [31]. In this case the relative χ-nucleon transverse velocity operator v ⊥ χN is conveniently rewritten as [29], where the first term v ⊥ χT is the χ-nucleus transverse velocity operator (with matrix element equal to v χT − q/2µ T , where v χT is the initial χ-nucleus relative velocity and µ T is the χ-nucleus reduced mass), and the second term v ⊥ N T is the transverse relative velocity of the nucleon N with respect to the nucleus center of mass [29]. To simplify the notation and connect it to the usual notation in analyses of dark matter experiments, we write v without index for the relative χ-nucleus velocity v χT .
The differential cross section for dark matter scattering on a target nucleus of mass m T is given by where |M N R | 2 denotes the square modulus of the non-relativistic scattering amplitude M N R (related to the usual invariant amplitude M by M = 4m 2 T M N R ), and j χ and j N are the dark matter and nucleus spins, respectively. When averaged over initial spins and summed over final spins, |M N R | 2 gives a quantity P tot proportional to the total transition probability, which can be expressed as a combination of nuclear and dark matter response functions. In the most general case it takes the following form For completeness, we list the dark matter response functions and W τ τ ∆Σ can be evaluated for different target materials and isotopes using the Mathematica package of Ref. [31], or the approximate expressions provided in the appendix of Ref. [29], where the functions F τ τ IJ = 4πW τ τ IJ /(2j N + 1). The definition of these nuclear response functions is given in Eq. (41) of Ref. [31]. In our calculations we have rewritten the Mathematica package of Ref. [31] in FORTRAN, and used our own routines to calculate differential cross sections and scattering rates. In Eq. (2.3), y = (qb/2) 2 , where b is the oscillator parameter in the independent-particle harmonic oscillator model [31].
The differential rate of scattering events per unit time and per unit detector mass is obtained as where ξ T is the mass fraction of the nucleus T in the target material, ρ χ is the local dark matter density, and m χ is the dark matter mass. The angle brackets in Eq. (2.5) denote an average over the local dark matter velocity distribution, f , in the galactic rest frame boosted to the detector frame, namely where v e (t) is the time-dependent Earth velocity in the galactic rest frame, and v min (q) = q/2µ T is the minimum velocity required for a dark matter particle to transfer a momentum q to the target nucleus. In our calculations we consider two choices of f : a Maxwell-Boltzmann ) truncated at the local escape velocity v esc , and the anisotropic velocity distribution proposed in Ref. [36].

Statistical framework
In this section we introduce the statistical methods used to extract limits on the strength of the dark matter-nucleon effective interactions from present dark matter direct detection data. We present both a Bayesian approach [27,[37][38][39][40][41][42][43][44][45][46][47][48] and a frequentist approach. This allows a comprehensive analysis of the multidimensional parameter space studied in this paper.
In the Bayesian analysis, our efforts are concentrated on reconstructing the posterior probability density function (PDF) of the model parameters, P(Θ|d). The posterior PDF depends on the n datasets d = (d 1 , . . . , d n ) considered in the analysis and on an array of m free parameters, denoted by Θ = (θ 1 , · · · , θ m ). The posterior PDF describes our degree of belief in a certain hypothesis after having considered the data (e.g., the validity of a specific configuration in parameter space). It is related to the likelihood function L(d|Θ) by Bayes' theorem, In this expression π(Θ) is the prior PDF, which describes our degree of belief in a certain hypothesis before having seen the available data. The Bayesian evidence E(d) is an important concept in the model comparison. Being independent of the model parameters, however, it simply plays the role of a normalization constant when performing parameter inference, as in the present analysis. The parameter space explored in our investigations is spanned by the coupling constants c τ i , with τ = 0, 1 and i = 1, 3, . . . , 11, the dark matter mass, m χ , and a set of "nuisance" parameters η introduced to model different sources of uncertainties affecting the interpretation of the data. There are two types of nuisance parameters relevant for the present analysis. A first type concerns the local dark matter space and velocity distribution. In our investigations -which focus on the form of the dark matter-nucleon interaction -we set the astrophysical nuisance parameters in a configuration known as the "standard dark matter halo" [49]. This is characterized by a truncated Maxwell-Boltzmann dark matter velocity distribution (in the galactic rest frame), with v 0 = 220 km s −1 , escape velocity v esc = 544 km s −1 and a local dark matter density ρ χ = 0.3 GeV cm −3 .
In two examples, presented in appendix A, we relax this assumption, considering a more general astrophysical setup characterized by 8 parameters describing a galactic dark matter component with an anisotropic velocity distribution [36]. The second type of nuisance parameters introduced in our statistical analysis is instead related to the presence of poorly known experimental quantities affecting the calculation of the expected dark matter direct detection signals (e.g., quenching factors, threshold effects, etc.). These parameters are introduced in the next section, describing the datasets included in this study. They are always treated as free parameters. Tab. 2 summarizes the free parameters considered in the following analysis. Following [31], we have introduced the mass scale m v = 246.2 GeV and the dimensionless quantities c τ i m 2 v . Among the many interesting pieces of information that can be obtained from the knowledge of the posterior PDF, we concentrate on 1D and 2D marginal posterior PDFs. These are calculated integrating the posterior PDF over the other model parameters. For instance, the 2D marginal posterior PDF of the parameters θ 1 and θ 2 (e.g., c 0 1 and m χ ) can be obtained integrating, i.e., marginalizing, over the remaining parameters as follows Limits on the coupling constants c τ i are then expressed in terms of x% credible regions, defined as the portions of the parameter space containing x% of the total posterior probability and such that P marg at any point inside the region is larger than at any point outside the region.
When the likelihood function is well approximated by a multivariate Gaussian and contains more information than the prior PDF, the associated credible regions tend to favor the portion of parameter space where the likelihood function is near its maximum L max . However, this is not true in general. For datasets containing an insufficient amount of information like those studied here, the integral in Eq. (3.2) Table 2. List of model parameters and nuisance parameters. Together with the type of prior, we also report the prior range and the reference from which this range has been taken. Gaussian prior PDFs are characterized by a mean lying at the center of the prior range and a standard deviation given by σ. The nuisance parameters q Na , ξ Xe , a COUPP , a PICASSO and a SIMPLE are introduced in Sec. 4 to model various types of detector uncertainties, whereas t max is the peak date used to describe the CoGeNT modulation signal [54]. Regarding the lower bound of the dark matter mass prior range we have considered 0.1 for PICASSO, SuperCDMS and CDMSlite, and 0.5 for the other experiments. Following [31], we have expressed the coupling constants in units of statistical indicator that is insensitive to these "volume effects" is the D-dimensional profile likelihood, which in the 2D case is defined as follows [50] L prof (d|θ 1 , θ 2 ) ∝ max While the profile likelihood does not admit a formal interpretation in terms of a probability density function, it can conventionally be used to construct approximate frequentist confidence intervals from an effective chi-square defined as ∆χ 2 eff ≡ −2 ln L prof /L max . Wilks' theorem guarantees that under certain regularity conditions the distribution of ∆χ 2 eff converges to a chi-square distribution with m degrees of freedom [50].
Within this approach to data analysis, all the experimental information is encoded in the likelihood function, and, to some extent, in the choice of the prior PDF (when calculating the posterior PDF), if specific assumptions are made in order to give more weight to certain portions of the parameter space. If not otherwise specified, in the analysis we use a Poisson likelihood to model the distribution of the observed data. This is an appropriate choice when the datasets consist of a small sample of k events, as for the recoil events detected (or searched for) by current dark matter direct detection experiments. Therefore, neglecting an irrelevant (for the parameter inference) constant term, our "default" choice for the likelihood function is [27,[38][39][40][41][42][43][44][45][46][47][48] − ln L(d|m χ , c, η, where µ S (m χ , c, η) represents the expected number of scattering events. For every experiment, we calculate µ S (m χ , c, η) within the effective theory of the dark matter-nucleon interaction. µ S (m χ , c, η) depends on the dark matter mass m χ , the coupling constants c = (c 0 1 , c 1 1 , . . . , c 0 11 , c 1 11 ) and a set of nuisance parameters η characteristic of the experiment under analysis. The likelihood in Eq. (3.4) also depends on the expected (or measured) number of background events µ B . For some of the experiments considered here, this background is a stochastic variable with a Gaussian distribution of variance σ 2 B and averageμ B . In this case, one can marginalize over the experimental background analytically, obtaining in the limit σ B μ B the form of the effective likelihood actually implemented as default choice in our analysis, namely In the second line we have expanded the Poisson factor around µ B =μ B and kept only the leading terms in this expansion. This procedure is justified by the fact that in the limit σ B μ B , the Gaussian factor in the integrand of Eq. (3.5) tends to the Dirac delta δ(µ B −μ B ).
Within our investigations, we employ log-priors both for the dark matter mass and for the coupling constants c, namely where Θ H is the Heaviside theta-function and θ min i and θ max i are the extrema of the prior ranges shown in Tab. 2. This assumption allows to sample the posterior PDF varying the model parameters within prior ranges spanning several orders of magnitude.
To sample the multidimensional likelihood surface and therefore reconstruct the posterior PDF and the profile likelihood of the model parameters, we use the Multinest program [55][56][57]. We use our own routines to calculate the scattering rates predicted by the dark matter-nucleon effective theory and to evaluate the likelihood function. Figures have been produced using the programs GetDist [58], Getplots [59] and Matlab. When calculating the profile likelihood we set the Multinest parameters to n live = 20000 and tol = 10 −4 , producing approximately 3 × 10 6 likelihood evaluations.
notice that in a real experiment the theoretical rate in Eq. (2.5) is not the quantity directly observed. In many cases the measurable energy E O is only a fraction of the true nuclear recoil energy E R deposited by a dark matter particle in the detector. Scintillators are an important example of detectors with this property. Moreover, the finite energy resolution and the limited efficiency E of a real detector can affect the observed direct detection rates. For a Gaussian energy resolution function (see Eqs. (4.6)-(4.7) for a non-Gaussian example), we write the observable differential rate of scattering events per unit time and per unit detector mass as follows In this expressionÊ O is the actually observed energy, whereas E O is the energy potentially measurable. The latter coincides with the former only in the limit of infinite experimental resolution. The energy dispersion σ is in general an energy dependent quantity. When not otherwise specified, we consider the time average of Eq. (4.1). The total number of events µ S (m χ , c, η) in the signal region [Ê 1 ,Ê 2 ] is then calculated integrating Eq. (4.1) over this energy range and multiplying the result by the experimental exposure M T (in kg-days, e.g.),

CDMS-Ge
The Cryogenic Dark Matter Search (CDMS II) experiment employs 19 germanium and 11 silicon detectors operating at cryogenic temperatures to search for dark matter through the observation of phonons and ionization. The final exposure of the experiment, in an analysis focused on a subset of the operating germanium detectors, is of 612 kg-days, corresponding to four periods of stable data taking between July 2007 and September 2008 [60]. In this period, the CDMS collaboration has observed two events in the acceptance region 10 -100 keV at recoil energies 12.3 keV and 15.5 keV, with an expected number of background surface events equal to 0.9 ± 0.3. We include this information in our analysis using the default likelihood (3.5), with k = 2,μ B = 0.9 and σ B = 0.3. To evaluate this expression we assume a maximum experimental efficiency of 32% at 20 keV, linearly decreasing towards lower and higher energies, reaching the value of 20% at 10 keV and at 100 keV. The energy resolution adopted in the calculations features a dispersion σ = 0.2. For the CDMS experimental apparatus

CDMS Low Threshold
The CDMS collaboration has also performed a low-threshold analysis of the data collected during six runs between October 2006 and September 2008. In their analysis the recoil energy threshold was lowered to 2 keV, while keeping the upper bound of the signal region at 100 keV [61]. Below 10 keV the discrimination between nuclear and electron recoils degrades and leads to a higher expected number of background events. Within this analysis, one of the germanium detectors, namely the fifth detector in the first tower of the CDMS detector array (T1Z5), has observed 38 candidate events within the signal region (36 of which between 2 keV and 20 keV; see Fig. 2 in Ref. [61]). The CDMS collaboration has identified three possible sources of background contamination, namely zero-charge events, surface events, and bulk events, which can explain 75% of the observed candidate events [61]. We include the results of this analysis in our investigations adding to the total likelihood a term of type (3.5) with k = 36 (considering the interval 2 -20 keV as signal region, as in Ref. [14])μ B = 36 × 0.75 and σ B = 0. Within our analysis we assume the efficiency shown in the inset of Fig. 1 in Ref. [61], an energy dependent energy resolution featuring σ = 0.293 2 + (0.056E O ) 2 , an exposure for T1Z5 of 241/8 kg-days, and E O = E R .

SuperCDMS
The SuperCDMS experiment is an upgrade of CDMS II which features new hardware devices interfaced with fifteen 0.6-kg cylindrical germanium crystals forming five towers containing three crystals each. The SuperCDMS collaboration has recently presented data recorded between October 2012 and June 2013 by a subset of 7 germanium detectors, corresponding to a total exposure of 577 kg-days [62]. This analysis has identified 11 dark matter candidate events passing the three levels of data-selection criteria introduced by the collaboration to discriminate candidate signal events from background events within the predefined signal region 1.6 -10 keV. Tab. 1 of Ref. [62] provides details regarding the number of events recorded by the 7 germanium detectors, and the number of background events expected for each detector separately. We include these data in our analysis adding a contribution of the form (3.5) to the total likelihood for each SuperCDMS detector, except for the detectors T5Z2 and T5Z3 for which the estimated number of background events seems to be significantly smaller than the number of actually observed candidate events (contrary to the other detectors).
To evaluate these contributions to the likelihood function we set µ B and σ B as in Tab. 1 of Ref. [62] and calculate µ S (m χ , c, η) for each detector using Eq. (4.1) with σ = 0.3 and the detector efficiency shown in Fig. 1 of Ref. [62], assuming an average exposure per detector equal to 577/7 kg-days.

CDMSlite
In a previous analysis the SuperCDMS experiment has collected data during 10 live days of dark matter search, using a single iZIP detector operating in a different mode (compared to previous studies) that yielded significantly better sensitivity to dark matter candidates of mass less than 10 GeV. This new operating mode is called CDMS Low Ionization Threshold Experiment, or simply CDMSlite [63]. The nuclear recoil energy threshold associated with this experimental configuration is of 170 eV ee , corresponding to 841 eV nr as one can see solving the non-linear equation [63] which relates the true nuclear recoil energy E R (measured in keV nr ) to the observable energy E O (measured in keV ee ). In this expression eV b = 69 eV and ε γ = 3 eV, whereas Y (E R ) is the ionization yield. The Lindhard model predicts for the latter [63] Y where g(ε) = 3ε 0.15 + 0.7ε 0.6 + ε, ε = 11.5E R Z −7/3 andk = 0.157 for a germanium target.
No background subtraction was applied to the collected data. This analysis has found that the average rate of nuclear recoils in the CDMSlite detector is 5.2 ± 1 counts/keV ee /kg-day between 0.2 and 1 keV ee , and 2.9 ± 0.3 counts/keV ee /kg-day between 2 and 7 keV ee . To include this information in our analysis, we have first calculated the expected average count rates in the CDMSlite detector using Eq. (4.1), with the efficiency reported in the inset of Fig. 1 in Ref. [63] and σ → 0. Then, we have added a term to the total likelihood given by where R [0.2,1] and R [2,7] represent the average rates between 0.2 and 1 keV ee , and between 2 and 7 keV ee , respectively.

XENON100
The XENON100 experiment uses liquid xenon to search for dark matter through the detection of ionization and scintillation signals produced by dark matter interactions in the active volume of the detector. in Ref. [64] the XENON100 collaboration has presented data collected in 13 months during 2011 and 2012, with an effective exposure of 34×224.6 kg-days. The ionization signal (S2) and the direct scintillation signal (S1) are both detected by arrays of photomultipliers (PMTs), and measured in numbers of photoelectrons (PE). The expected number of S1 photoelectrons ν(E R ) produced by a nuclear recoil of energy E R is given by 28±0.04 PE/keV ee is the light yield, S ee = 0.58 and S nr = 0.95 are the electric field scintillation quenching factors for electron and nuclear recoils, and finally, L eff (E R ) is the energy dependent scintillation efficiency. The actually produced number n of photoelectrons for a given energy E R is subject to Poisson fluctuations around ν(E R ), and to uncertainties related to the scintillation efficiency, which has not been measured at energies below 3 keV nr . We model these uncertainties introducing a nuisance parameter ξ Xe , first proposed in Ref. [41] to logarithmically extrapolate the scintillation efficiency toward low energies. This gives is the best-fit scintillation efficiency reported in Fig. 1 of Ref. [65]. Importantly, the observed number of photoelectrons S1 does not coincide with n when the finite resolution of the detector photomultipliers is taken into account. The analysis of Ref. [64] uses S1 to reconstruct the recoil energy of the candidate dark matter signal events and the ratio S2/S1 to discriminate signal events from background events. In that analysis, XENON100 observed 2 candidate signal events in the pre-defined nuclear recoil energy range 6.6 -30.5 keV nr (corresponding to S1 in the range 3 -30 PE). This observation is consistent with the background expectation of 1.0 ± 0.2 events. To include XENON100 in our analysis, we calculate the differential spectrum of the variable S1, first converting the recoil energy spectrum (2.5) into the spectrum of the expected number of photoelectrons n, namely and then convolving the resulting expression with a Gaussian filter, to model the resolution of the detector photomultipliers: Gauss(S1|n, Following Ref. [66], in the first step we have used a Poisson PDF of mean ν(E R ) to sample the number of actually produced photoelectrons n associated with a given recoil energy E R . The Gaussian filter employed in the second step has mean n and variance nσ 2 PMT , with σ PMT = 0.5. The detector efficiency E(S 1 ) relevant for this analysis is shown in Fig. 2 of Ref. [65]. Integrating Eq. (4.7) between 3 and 30 photoelectrons and multiplying the result by the experimental exposure, we finally obtain the expected number of signal events, µ S (m χ , c, η). This prediction is then used to evaluate the likelihood (3.5) with k = 2,μ B = 1 and σ B = 0.2.

XENON10
In a second study the XENON collaboration reanalyzed the data from a 12.5 live day dark matter search, collected between August 23 and September 14 2006, using the S2 signal only to measure the nuclear recoil energy of the detected events [67]. The relation between the S2 variable (measured in PE) and the observed nuclear recoil energy is Fig. 1 of Ref. [67]. This type of analysis allows a very low recoil energy threshold (about 1.4 keV nr ), increasing thus the detector sensitivity to low mass dark matter candidates, even if the detector ability in discriminating and rejecting electromagnetic background events is reduced within this setup. The experimental exposure corresponding to these data is 12.5×1.2 kg-days. Within this analysis XENON10 observed 23 candidate events in the signal region 1.4 -10 keV nr . XENON10 has also observed several dozens of single S2 electron events at lower energies, the origin of which is not clear yet. In our analysis we treat the 23 events in the signal region as possible dark matter candidates and estimate the expected nuclear recoil energy spectrum in XENON10 as (see Eq. 5.5 in Ref. [68]) dR assuming a constant efficiency E(Ê O ) = 0.94. In addition, we also include the possibility that a source of background events of unknown origin and characterized by a large error contributes to the observed events. The XENON10 contribution to the total likelihood is then estimated using Eq.

LUX
The Large Underground Xenon (LUX) experiment is a dual-phase (liquid and gas) timeprojection chamber with 250 kg of active volume. Similarly to XENON100, LUX searches for dark matter through the observation of prompt scintillation (S1) and ionization electrons, extracted into the gas portion of the detector, where they produce electroluminescence (S2) [69]. In the case of LUX, the conversion between nuclear recoil energy (in keV nr ) and number of photoelectrons can be extracted from the panel (b) of Fig. 3 in Ref. [69]. The recent first data release consists of 85.3 live days of dark matter search data, collected between April 21 and August 8 2013. In this period, LUX observed 160 events with S1 between 2 and 30 PE, only one of which is (slightly) below the mean of the Gaussian fit to the nuclear recoil calibration events reported in Fig. 4 of Ref. [69], with an expected number of background events in that region of 0.64 ± 0.16 events. We include this information in our analysis adding a term of the form (3.5) to the total likelihood, with k = 1,μ B = 0.64 and σ B = 0.16. We calculate µ S (m χ , c, η) integrating Eq. (4.7) between 2 and 30 PE, and assuming σ PMT = 0.37, an exposure of 250 × 85.3 kg-days, and the experimental efficiency reported in Fig. 9 of Ref. [69], multiplied by an additional factor 1/2, corresponding to the 50% nuclear recoil acceptance quoted by the LUX collaboration.

COUPP
The Chicagoland Observatory for Underground Particle Physics (COUPP) experiment seeks to observe bubble nucleations arising from dark matter scattering in a superheated liquid. The latest results from a 4.0 kg CF 3 I bubble chamber operating from September 2010 to August 2011 at the SNOLAB deep underground laboratory have been reported in Ref. [51]. For this experimental apparatus, the probability that an energy E R nucleates a bubble above a threshold energy E th is given by [51] This probability depends on the target nucleus. The constant α T , with T =C,F,I, is determined by fitting the above expression to the observed rates of single, double, triple and quadruple bubble events in test runs performed using neutron sources. We calculate the expected number of dark matter scattering events with an energy larger than E th in the COUPP detector as follows [51] µ S (m χ , c, η) = (E th ) where (E th ) is the threshold dependent experimental exposure multiplied by the bubble detection efficiency. The values relevant for the present analysis are (7.8 keV) = 55.8 kgdays, (11 keV) = 70 kg-days and (15.5 keV) = 311.7 kg-days [51]. COUPP operated in three different experimental configurations, corresponding to bubble nucleation threshold energies of 7.8 keV nr , 11 keV nr and 15.5 keV nr respectively. With these threshold energies, COUPP observed k = 2, k = 3 and k = 8 events respectively for each experimental configuration. The corresponding estimated number of background events associated with α-decays in the materials surrounding the CF 3 I volume is µ B = 0.8, µ B = 0.7 and µ B = 3 events, respectively. For each threshold energy, we add a term of type (3.5) to the total likelihood if µ S (m χ , c, η)+ µ B > k, and a large negative constant otherwise, which corresponds to considering the value of k as an upper bound only. In all cases we set σ B = 0. Regarding the parameter α i , following [33], we assume α I → +∞ (i.e., perfect efficiency for bubble nucleation), α C = 0 (i.e., no bubble nucleation) and finally, α F ≡ a COUPP , treating the latter as a nuisance parameter with a Gaussian prior, as shown in Tab. 2. We do not consider dark matter scattering on Carbon, since for this nucleus nuclear form factors are not available for all the nuclear responses emerging in the effective theory of Ref. [31].

PICASSO
The PICASSO experiment searches for dark matter using superheated liquid droplets made of C 4 F 10 [52]. In the last run, PICASSO operated in eight different experimental configurations, corresponding to the following bubble nucleation threshold energies (in keV nr ): 1.7, 2.9, 4.1, 5.8, 6.9, 16.3, 38.8 and 54.8. For each energy, the collaboration reported the observed rateR i (for i = 1, . . . , 8) of dark matter scattering events above threshold including the associated experimental errors σ i (see Fig. 5 of Ref. [52]). We calculate the expected scattering rate R i in the energy range [E th , +∞] using Eq. (4.10) but setting = 1 and assuming α F = a PICASSO , where a PICASSO is the nuisance parameter described in Tab. 2. Similarly to the COUPP experiment, we focus on dark matter scattering off fluorine only, since carbon form factors are not available for all the nuclear responses proposed in Ref. [31]. The PICASSO contribution to the total likelihood is then obtained assuming a Gaussian likelihood function for the eight PICASSO data points, namely (4.11)

SIMPLE
Similarly to COUPP and PICASSO, the Superheated Instrument for Massive ParticLe Experiments (SIMPLE) uses superheated liquid detectors made of C 2 ClF 5 to search for bubble nucleations produced by dark matter scattering in the detector volume [53]. We focus here on the Stage 2 data, collected with an effective experimental exposure of 6.71 kg-days, after cutting data obtained at pressures greater than 2.20 bar. In this run the experiment operated with a bubble nucleation threshold energy of 8 keV nr . Analyzing these data, the SIMPLE collaboration observed 1 candidate dark matter event. This result is consistent with an expected number of background events equal to 2.2 ± 0.3. To estimate the expected number of dark matter scattering events in the SIMPLE detector, we use an equation analogous to Eq. (4.10), with the parameter α F ≡ a SIMPLE treated as a nuisance parameter with a Gaussian prior, as shown in Tab. 2. To use the density matrix approach of Ref. [31] in the calculation of the nuclear form factors, we restrict our analysis to dark matter scattering off fluorine (i.e., we set α Cl = α C = 0). The SIMPLE contribution to the total likelihood is then of type (3.5) with k = 1,μ B = 2.2 and σ B = 0.3.

DAMA
DAMA uses highly-radiopure thallium-doped NaI scintillators to detect dark matter through the observation of an annual modulation in the measured nuclear recoil energy spectrum of the target sodium and iodine nuclei [70]. The modulation signal is expected as a consequence of the Earth's motion through the Milky Way dark matter halo, which sinusoidally modulates the flux of dark matter particles impinging on the DAMA detector, with a period of one year [49]. In the conventional dark halo model, the modulation is expected to be at a maximum on June 2 and at a minimum on December 2. Combining the data collected over 7 annual cycles by the DAMA/NaI configuration of the experiment with the data of the 6 annual cycles recorded by its upgrade DAMA/LIBRA, the total exposure of the DAMA experiment reaches 1.17 ton×year. As for any direct detection experiment relying on scintillators, only a certain fraction of the true nuclear recoil energy deposited by dark matter particles in the DAMA NaI crystals is directly accessible to the photomultipliers installed in the experimental apparatus. The measurable scintillation energy E O (in keV electron equivalent units, i.e., keV ee ) is related to the true nuclear recoil energy E R (in keV nuclear recoil units, i.e., keV nr , or simply keV) by a quenching factor, denoted by q Na and q I for sodium and iodine nuclei, respectively. The value of the quenching factor is uncertain, and different choices have been explored in the literature. In our fits we set q I = 0.09, and treat q Na as a nuisance parameter varying with a Gaussian prior in the range shown in Tab. 2. In the latter case E O = q Na E R . To evaluate Eq. (4.1) we also need the DAMA energy resolution, which is In this paper we use the annual modulation data reported in Fig. 6 of Ref. [70], including the first 12 energy bins only, since for energies larger than 8 keV ee DAMA has not observed any statistically significant modulation effect. For these data we assume the likelihood function where the expected annual modulation amplitude as a function of the energy bin lower bound E i O is calculated as follows [49] S (4.14) Here σ i are the errors associated with the 12 datapointsŜ m (Ê i O ) and ∆Ê O = 0.5 keV ee is the width of the energy bins.

CoGeNT
The CoGeNT experiment searches for a dark matter signal employing p-type point contact germanium detectors. After 1136 live days of data taking, the collaboration has recently published a new measurement of the observed nuclear recoil energy spectrum and of its time dependence [54]. The data show an exponential-like irreducible background of events that cannot be associated with cosmogenically activated nuclei decaying via L-shell or K-shell electron capture. In addition, the data collected also provide moderate evidence in favor of a modulation signal, with a phase subject to large uncertainties implying a peak date of t max = 106 ± 24 days. Using the same data, Ref. [71] has found that the ratio between the modulation amplitude and the total unmodulated signal (i.e., the fractional amplitude) is equal toÂ = (12.4 ± 5)%, with a dark matter signal estimated to be 35% of the total unmodulated signal. For the CoGeNT experimental apparatus, the relation between the observable energy and the true nuclear recoil energy is E O = 0.199 × E 1.12 R . In the fits we employ this relation and the recoil energy spectrum extracted from Fig. 10 in Ref. [54] as follows. We find the observed energy spectrum dR/dÊ O at the energyÊ i O by subtracting the best fit background estimate in Fig. 10 to the observed datapoints (labeled here by an index i). We denote by σ i the error associated to this spectrum. When fitting these data, we have set E = 1 and σ → 0 in Eq. (4.1). For this data sample we have assumed a multivariate Gaussian contribution to the total likelihood, namely Here σ A = 5%, and the second term in this expression takes into account the information on the CoGeNT annual modulation. When calculating the expected fractional amplitude A theory , we have treated the peak date t max as a nuisance parameter with a Gaussian prior as shown in Tab. 2.

Limits on the dark matter-nucleon interaction strength
We now compare the predictions of the dark matter theory introduced in Sec. 2 to the data of the previous section, using the statistical tools summarized in Sec. 3. We assume for definiteness that the dark matter particle has spin j χ = 1/2. In this section, we focus on experiments compatible with a null result. In the next section we study experiments with a dark matter signal, using the same theoretical and statistical frameworks.

Limits from single experiments
We start with a detailed analysis of the LUX data, to illustrate the many physical effects and computational challenges which can be encountered when studying a parameter space of large dimensionality, like in this work. The top-left panel of Fig. 1 shows the results of a fit where we have considered as free parameters c 0 1 and m χ only, setting to zero all the remaining couplings. This corresponds to the standard case in which the LUX data are interpreted in terms of spin-independent interactions, with one important difference, however, namely the fact that instead of presenting Confidence Levels (CL) in the m χ -σ SI p plane, where the latter is the χ-proton cross section, we present the 2D posterior PDF and its associated 99% Credible Region (CR) in the related m χ -c 0 1 plane. The connection is, from Eqs. (2.2)-(2.3) keeping the isoscalar part only, where µ N = m χ m N /(m χ + m N ) is the reduced χ-nucleon mass. For reference, the familiar spin-dependent cross section σ SD N is related to c 0 4 by In the top-left panel of Fig. 1 we observe a smooth 99% CR contour and a posterior PDF that grows below this contour to reach a plateau of approximately constant posterior probability. The calculation to produce this plot required O(10 4 ) likelihood-function evaluations, even when additional nuisance parameters are included in the fit to model various sources of uncertainty (as done for instance in the case of XENON100). In a second analysis, we fit the LUX data varying the ten couplings c 0 i , with i = 1, 3, . . . , 11, and setting to zero the analogous isovector couplings, assuming isospin-conserving interactions. The top-central panel of Fig. 1 shows the 2D marginal posterior PDF and the associated 99% CR in the m χ -c 0 1 plane. Important differences are observed between this PDF (marginalized over 9 parameters) and the previous one in the top-left panel (with the 9 parameters set to zero). The marginalized PDF is peaked at low masses, with a 99% CR contour consisting of two disconnected "islands," one, more pronounced, at small masses, and the other at large masses. Qualitatively, the appearance of the "islands" can be analytically understood as a volume effect emerging during the marginalization process (see discussion around Eq. (3.3)). In fact, if correlations between different c 0 i can be neglected, the 2D marginal posterior PDF P marg (c 0 i , m χ ) factorizes as follows In the first line,  Fig. 1, where the 2D marginal PDF P marg (c 0 j , m χ ) is highly peaked at low masses, has a minimum at m χ ∼50 GeV, and slowly increases at larger masses. The volume effects are however uncomfortable.
As already mentioned, a statistical indicator insensitive to the volume effects is the profile likelihood. The bottom-left panel of Fig. 1 shows the 2D profile likelihood, together with the associated 95% CL contour, obtained from a fit of the LUX data where we vary c 0 1 and m χ only. Below the 95% CL contour, the profile likelihood increases, it reaches a region of maxima corresponding to an expected count rate µ S (m χ , c, η) + µ B ∼ 1, and finally, it decreases towards a plateau associated with µ S (m χ , c, η) = 0. Importantly, when varying a single coupling, the 2D posterior PDF and the 2D profile likelihood in the m χ -c 0 1 plane have a very similar shape. The bottom-central panel shows the 2D profile likelihood extracted from a fit of the LUX data performed varying all the couplings and the dark matter mass simultaneously. Contrary to the case in which a single coupling is varied, the region of maxima extends everywhere below the 95% CL contour, except at low masses, where the expected dark matter signal is below the experimental threshold and µ S (m χ , c, η) = 0 independently of c. The presence of an infinite plateau of maxima can be explained as follows: even when c 0 1 is very small, the value of µ S (m χ , c, η), which would be tiny if all the other couplings were zero, can be sufficiently large to satisfy the maximum condition µ S (m χ , c, η) + µ B ∼ 1, because of contributions associated with the other couplings. Only at sufficiently low values of m χ , µ S (m χ , c, η) = 0 independently of c because of the already mentioned threshold effects.
The top-right panel of Fig. 1 illustrates the two 1D posterior PDFs obtained marginalizing over the coupling constants, and the two 1D profile likelihoods associated with the four CR/CL contours described in the previous paragraphs. This figure shows the flatness of the profile likelihoods (modulo threshold effects) and the location of the peaks of the 1D marginal posterior PDFs.
The conclusions, illustrated here in detail in the important case of the LUX experiment, also apply -of course with obvious quantitative differences -to the other experiments considered in this paper, and to coupling constants different from c CDMSlite, COUPP, SIMPLE and PICASSO data. 3 To extract these contours we analyze each dataset independently and sample the posterior PDF by varying m χ and one of c 0 i at a time, in addition to the relevant nuisance parameters. We have verified that the limits obtained in this way in the planes m χ vs c 0 1 and m χ vs c 0 4 , match very well standard results usually presented in the planes m χ vs σ SI n and m χ vs σ SD n . In addition, from Figs. 2 and 3 one  In Figs. 2 and 3 we superimpose scatter plots of σ SI N to the 99% CR regions, using the color code illustrated in the figure to distinguish different values of this cross section. In each panel we construct the corresponding scatter plot as follows: (1) we sample 1,000 points in parameter space from the posterior PDF resulting from a fit of the LUX data (in Fig. 2) or of the COUPP data (in Fig. 3), where the free parameters are m χ and all the 10 couplings c 0 i (volume effects explain why the sample points are concentrated at low masses); (2) we calculate σ SI N = µ 2 N |c 0 1 | 2 /(4π) for each point; (3) we plot the projection of each point in the 2-dimensional space where the CR contours are shown. Since σ SI N ∝ |c 0 1 | 2 , we find an obvious correlation between σ SI N and c 0 1 . In all the other panels, however, the existence of a correlation between c 0 1 , and therefore σ SI N , and the other coupling constants c 0 i (i = 1) is not as obvious. The important point of correlations between the c i 's is investigated in the next section.  we focus here on experiments which led to a null result. The DAMA and CoGeNT experiments, which claim a signal, are discussed in the next section. Before describing our numerical results, we present some semi-analytic considerations on the correlation between the coupling constant c 0 i . The expected number of dark matter events µ S (m χ , c, η) is a quadratic function of the constants c 0 i , as follows from Eqs. . If we consider the correlation for one of these pairs (c 0 i and c 0 j , say) setting the other coupling constants to zero, the contours of constant µ S (m χ , c, η) for a given experiment at a given m χ and η are ellipses in the c 0 i -c 0 j plane. These ellipses can be obtained without random sampling in parameter space by writing

Global limits
where µ Sconst is the desired value of µ S (e.g., its upper limit) and the coefficients a ii , a ij ,  We exploit the Multinest program to explore the multidimensional parameter space of the dark matter-nucleon effective theory by simultaneously varying the 11 model parameters and the 4 additional nuisance parameters listed in Tab. 2. Our analysis is based on about 3 million likelihood evaluations.  wit, the black line in Fig. 6 is the graph of an ellipse in a log-log plane, namely the LUX upper-limit ellipse in Fig. 4 (top-left). The boundary of the dark region follows the black line in Fig. 6, and is itself an ellipse in the c 0 1 -c 0 3 plane. The correlation between c 0 4 and c 0 6 is not visible in Fig. 5 because the Multinest analysis is restricted to positive values of the c 0 i . It is comforting that the semi-analytic considerations at the beginning of this section and the full Multinest analysis give the same correlation pattern for the c 0 i 's. An interesting and important result is summarized in Fig. 7, which shows the 2D marginal PDFs and the 2D profile likelihoods from our global analysis of the direct detection

Analyzing a signal: DAMA & CoGeNT
This last section is devoted to a Bayesian analysis of the DAMA and CoGeNT data. Contrary to the analysis illustrated in the previous section, we now concentrate on the interpretation of two candidate dark matter signals. We approach this problem within the theoretical and statistical frameworks introduced in Secs. 2 and 3, respectively, focusing on the operators O 8 and O 11 only, and leaving the general case for future work. We start with an analysis of possible degeneracies between different coupling constants.   obtained from the same analysis of the DAMA-Na and CoGeNT data described above, i.e., varying m χ , c 0 8 and c 0 11 simultaneously (dotted contours). The associated 68% and 90% CR contours are characterized by long tails extending toward the direction of zero coupling constants. These tails are related to the existence of the two limiting solutions to the problem of fitting the data described in the previous paragraph. In addition, Fig. 9 also shows the 2D marginal posterior PDFs resulting from an analysis of the DAMA-Na and CoGeNT data where we have separately considered as free parameters either c 0 8 only or c 0 11 only (solid contours). The solid contours are at the end of the dotted regions.
In this example, the regions favored by DAMA-Na and CoGeNT in the plane dark matter mass vs interaction strength are well separated, discouraging therefore a global fit of the two datasets. An analysis of the remaining couplings is left for future work.
We conclude this section investigating whether the annual modulation signal observed by DAMA can be due to dark matter scattering on iodine. Also in this case we focus on the operators O 8 and O 11 as an illustrative example.  Figure 10. Fits of the DAMA data performed varying a single coupling and m χ . We have assumed dark matter scattering on iodine only. We plot the resulting 2D marginal posterior PDF in the planes m χ vs interaction strength for two operators studied in this paper. In both cases the interpretation of the DAMA results in terms of dark matter scattering on iodine is ruled out.
annual modulation signal cannot be ascribed to dark matter interactions with iodine nuclei described by the operators O 8 and O 11 .

Conclusion
In this paper we have presented the first comprehensive analysis of the dark matter-nucleon effective interactions where the coupling constants, the dark matter mass, and additional nuisance parameters, are simultaneously considered as free parameters. To study experimental constraints on the multidimensional parameter space of coupling constants, we have implemented a Bayesian and a frequentist approach to extract credible and confidence regions from a varied sample of complementary direct detection data, including the recent LUX, CDMSlite and SuperCDMS results. We have extracted an upper bound on the 10 coupling constants characterizing the theory of heavy spin-0 and spin-1 mediators as a function of the dark matter mass, and in the limit of isospin-conserving interactions. We have calculated the posterior PDF and the profile likelihood of the model parameters, and shown credible regions and confidence levels in the planes dark matter mass vs interaction strength, marginalizing (in the Bayesian approach) and profiling out (in the frequentist approach) the uncertain or irrelevant model parameters. For the still limited experimental data currently available, the Bayesian and frequentist methods turn out to be complementary statistical indicators, in the sense that the Bayesian method is faster but subject to artificial volume effects that depend on the prior, while the frequentist method is free of volume effects but is computationally slower.
We find that present direct detection data contain sufficient information to simultaneously constrain not only the familiar velocity-and-momentum-independent interactions (i.e., the spin-independent operator O 1 = 1 χ 1 N and the spin-dependent operator O 4 = S χ · S N ), but also the remaining velocity and momentum dependent couplings predicted by the dark matter-nucleon effective theory. The interaction most severely constrained by the current data is O 1 , followed by the interaction O 11 and then O 3 , O 4 , and O 8 (see Fig. 7). Notice that the relatively strong constraints on O 11 have been observed in [32] and indirectly in [33] (through their relativistic operator O 2 ), but have not been considered in other studies [24].
In addition, we have found that strong correlations exist between c 0 1 and c 0 3 and between c 0 4 and c 0 6 , associated with the interaction operators O 1 = 1 χ 1 N , O 3 = −i S N · q × v ⊥ /m N and O 4 = S χ · S N , O 6 = S χ · q S N · q/m 2 N . Other correlations between the c 0 i 's are negligible either because there is no interference term between the corresponding operators or because they are suppressed by the smallness of the nuclear structure functions and/or momentum transfer.
Presenting our results we have also described the difficulties found when exploring the effective-theory parameter space of large dimensionality. For instance, we have found that the marginalization process introduces important volume effects, which significantly alter the shape of the resulting 2D marginal posterior PDFs. It is therefore important, in order to assess reliable upper limits on the strength of the dark matter interactions, to calculate the profile likelihood as well, though this calculation is in general computationally very demanding and in our case it has required about 3 million likelihood evaluations.
In a last part of the paper we have studied two candidate dark matter signals, namely those reported by the DAMA and CoGeNT collaborations. We have approached this problem within the same theoretical and statistical frameworks used in the study of the exclusion limits. In this analysis we have considered two interaction types only for DAMA and Co-GeNT, leaving the general case for future work. Analyzing a dark matter signal by varying two coupling constants and the dark matter mass, degeneracies between different couplings are apparent. They are associated with the existence of distinct solutions to the problem of fitting the data. For the two interactions types considered in this study, the regions favored by CoGeNT and DAMA in the plane dark matter mass vs interaction strength are well separated, discouraging therefore a global fit of the two datasets performed within this setup.
In summary, we have proposed a systematic approach to the study of dark matternucleon effective interactions. Our approach is based on the calculation of the posterior PDF and/or the profile likelihood in the full multi-dimensional parameter space characterizing the theoretical framework. This strategy allows the extraction of global limits on the model parameters accounting for a variety of theoretical and experimental uncertainties through the introduction of nuisance parameters, which are then marginalized or profiled out during the calculation. In addition, this approach allows an interpretation of the direct detection data which is not biased by having assumed a priori the form of the dark matter-nucleon interaction. We are confident that this general and flexible approach to the analysis of dark matter direct detection data will be particularly fruitful to exploit the results of the next generation of direct detection experiments. Standard Halo MCMC Figure 11. 99% credible regions in the planes m χ vs c 0 1 (left panel) and m χ vs c 0 5 (right panel) resulting from an analysis of the XENON100 data assuming the galactic model of Ref. [36], instead of the standard dark matter halo (in this analysis, the coupling constants not shown in the figures are set to zero, and the PDFs are marginalized over the 8 astrophysical model parameters). Changing the astrophysical assumptions modifies the 2D marginal posterior PDFs only moderately.

A Changing the astrophysical assumptions
In this paper we have assumed the standard dark matter halo to extract limits on the couplings c 0 i from present direct detection data. We check here with two examples to which extent our conclusions would have been affected by having assumed a different astrophysical configuration. In this appendix, we assume the galactic model studied in depth in Refs. [36,72,73]. It features 8 parameters describing a galactic bulge, a stellar disk and a spherical dark matter halo. This model [36] has been generalized to include an anisotropic velocity distribution for the Milky Way dark matter particles.
In this appendix we calculate the 2D marginal posterior PDFs in the planes m χ vs c 0 1 and m χ vs c 0 5 associated with the XENON100 data assuming the more complex astrophysical configuration in [36] and marginalizing over its 8 astrophysical parameters. For these parameters we assume the prior PDFs shown in Fig. 4 of Ref. [36]. Having computed the 2D marginal posterior PDFs, we then compare the corresponding 99% CR contours with those obtained assuming a standard dark matter halo. The results of this analysis are shown in Fig. 11. The left panel refers to the operator O 1 , as an example of velocity/momentum independent operator, whereas the right panel refers to the operator O 5 , which is instead a velocity/momentum dependent operator. As one can see from this figure, changing the astrophysical assumptions modifies the 2D marginal posterior PDFs only moderately. Our interpretation of present direct detection data is more sensitive to the assumptions made regarding the underlying dark matter-nucleon interaction. Notice however that the mean galactic model found in Ref. [36] and the standard dark matter halo have very similar local dark matter densities. A more drastic modification of the local density would have induced more significant changes in the 2D marginal posterior PDFs reported in Fig. 11.

B Dark matter response functions
In the following we list the dark matter response functions used in the calculations presented in this paper. These have been obtained from the ones derived in Ref. [31] setting to zero the couplings c τ 12 , . . . , c τ 15 , with τ = 0, 1: