Marco Cirelli
email: marco.cirelli@cea.fr

Eugenio Del Nobile
email: delnobile@physics.ucla.edu

Eugenio Del Nobile

Paolo Panci
email: panci@cp3-origins.net

Tools for model-independent bounds in direct dark matter searches (updated including LUX 2013 data) (updated including SuperCDMS 2014 data)

published or not. The documents may come

Tools for model-independent bounds in direct dark

matter searches

Introduction

Since several decades, cosmological and astrophysical evidences accumulate on the existence of Dark Matter (DM), thanks to its gravitational influence at different scales: at the galactic scale, with the flattening of rotation curves, at the galaxy cluster scale, with weak lensing measurements, and the cosmological scale, with CMB (Cosmic Microwave Background) and LSS (Large Scale Structure) observations. See e.g. [1] for a review of all these aspects. On the other hand, an explicit (non gravitational) manifestation of DM is yet to be identified. As well known, there are three main approaches: direct detection searches aim at detecting a feeble kick to an atomic nucleus in extremely low background environments; indirect detection techniques aim at unveiling possible excess cosmic rays produced by the annihilations or decay of DM particles in the Galaxy or beyond; collider searches aim at identifying signatures of the production of the DM particles at high energy particle accelerators, most notably the Large Hadron Collider (LHC).

The observables which are relevant in each one of these approaches, of course, all pertain to the interactions of the unknown Dark Matter with the ordinary matter, which is what we measure 'in the laboratory'. In direct detection (DD), the important quantity is the cross section for scattering of DM particles off a nucleon, in a deeply non-relativistic regime (since the speed of the incoming DM particle is of the order of v ∼ 10 -3 c, with c the speed of light) characterized by a momentum transfer of a few (tens to hundreds of) MeV. As a consequence, as stressed recently by [2][3][4], the most useful framework for DD analyses is in terms of a basis of non-relativistic operators (given in Eq. ( 1) below), which parameterize most efficiently the actual interactions among DM and ordinary matter in terms of a very limited number of non-relativistic degrees of freedom. This approach turns out to be very generic and powerful, given that, in spite of the many possible DM candidates one can envision (complex or real scalar, Dirac or Majorana fermion, vector), and of the different interactions they can have with matter, the scattering process can always be expressed in terms of a very limited number of non-relativistic degrees of freedom. As a matter of fact, one can write a basis of non-relativistic (NR) operators; the matrix element of any process of elastic scattering between the DM and the nucleon, regardless of the high energy interaction responsible for the scattering, can then be expressed as a combination of these NR operators. Moreover, most of the field operators usually encountered in both renormalizable and effective theories reduce to only a subset of the NR operators, given in Eq. ( 1) below.

In order to make contact with the actual experimental observables (for instance: the number of events in a detector like Xenon100), however, one needs to perform two additional steps. Firstly, one has to compute the nuclear response caused by each one of the NR operators, or a combination of those. This can be done using the formalism clearly spelled out in Fitzpatrick et al. in [3]. Their work provides form factors F (N,N ) i,j which take into account the non-relativistic physics of the DM-nucleus interaction, and encode all the nuclear information as well as the dynamics of the DM-nucleus interaction, for a number of different nuclei of interest in DD experiments. We will heavily rely on their results, as we will discuss in the following. Secondly, one has to concretize the computed nuclear response to a specific experiment, considering detection efficiency, cuts, acceptances and the energy resolution of the detector. Moreover, one has to make astrophysical assumptions on quantities such as the density and velocity distribution of DM particles in the local halo, and fold them in the computation. This is where we come in: we will present a formalism Experiment Section Reference Updated as of Xenon100

4.1 [START_REF] Aprile | XENON100[END_REF][START_REF] Aprile | XENON100 Coll[END_REF] 03 in terms of integrated form factors

F(N,N ) i,j
which allow to compute the constraints from DD experiments taking into account the concrete experimental details and the astrophysical assumptions. This constitutes the main result of our work. Now, Dark Matter being a particle physics problem, any 'model' of the interactions of DM with ordinary matter will be expressed in terms of operators involving the degrees of freedom proper of the high-energy theory, i.e. quarks and gluons. If one works in the framework of a specific theory comprising a DM candidate (e.g. Supersymmetry or Extra Dimensions), one generally knowns the explicit form of such terms, however complicated and parameter-loaded they might be. Conversely, if one is ignorant or agnostic about the underlying theory, one can express the interaction 'Dark Matter -ordinary matter' in terms of effective operators, which are descriptions of particle interactions at energies lower than the masses of the interaction mediators. In both cases, anyway, the high-energy physics will necessarily reduce to NR operators. Indeed, in the second part of our paper, we make this reduction explicit within the relativistic effective operator approach, reminding how one passes from a set of general high-energy effective operators to the low-energy NR operators.

Summarizing, the goal of our paper is to describe a streamlined method and to provide the needed tools to compute the bounds from the main current DD experiments without having to perform explicitly the steps described above. The starting point can be a description in terms of relativistic effective operators (in which case we provide a self-contained guide to obtaining the bounds) or any high energy theory, in which case it is up to the model builder to reduce it to the relevant NR operators on which our tools can be employed. The concrete output consists of a set of Test Statistic (TS) functions, one per each featured experiment, and a set of scaling functions (based on the integrated form factors) allowing to translate a benchmark bound into a bound for an arbitrary combination of operators. These numerical products are provided on this website: www.marcocirelli.net/NRopsDD.html. It is worth stressing that these results apply mainly to the case of DM-nuclei elastic scattering. However, our method is quite general and in principle similar tools can also be developed for other classes of models (e.g. inelastic DM scattering). We encourage therefore independent extensions in this direction. Before moving on, three important remarks concerning the experiments are in order. (i) In this paper we focus exclusively on bounds from DD experiments which have reported null results. Of course there is another set of DD experiments whose results can be interpreted as claims for an allowed DM region. While a method similar to ours can be constructed for the analysis of these experiments, we leave this to future work. (ii) The TS functions that we provide for each experiment reflect the current status of results: they will change whenever new results are released and we will update them accordingly. On the other hand, the integrated form factor (and, a fortiori, the scaling functions) change only if major changes occur in the experimental set-ups, e.g. in efficiencies, cuts, and thresholds. We do not foresee frequent updates for those, but we will consider them when necessary. (iii) Among the many very well performing experiments, we choose to focus on six of them, since these currently provide the strongest bounds. Our method, however, obviously can be applied to any set-up and indeed we encourage independent extensions to other experiments. Table 1 (updated in the arXiv version of this paper) lists the experiments we consider with the corresponding most updated reference.

The rest of this paper is organized as follows. In Sec. 2 we lay down the formalism of the NR operators and illustrate how the experimental observables are expressed in terms of it. In Sec. 3 we discuss how one derives constraints based on the experimental results and how one can rescale a general benchmark bound into a bound on a given DM model. In Sec. 4 we describe in some detail the experiments that we consider, in order to present the specific form that the integrated form factors take for each one of them. In Sec. 5 we discuss the reduction from the set of high-energy effective operators to the low-energy NR operators. In Sec. 6 we illustrate all the formalism with some explicit examples. Finally in Sec. 7 we conclude.

2 Phenomenology: from the NR operators formalism to the experimental observables

In this Section we introduce the formalism of NR operators, following [2,3], and we describe how to compute the experimental observables (essentially the number of events in the energy bins of a certain experiment) in terms of it. More precisely, we will write the differential event rate as a linear function of a manipulation of the form factors

F (N,N ) i,j
(provided by [3]), which take into account the non-relativistic physics of the DM-nucleus interaction, and encode all the nuclear information as well as the dynamics of the DM-nucleus interaction.

In a NR description of the elastic scattering of a DM particle χ with a nucleon N , the relevant degrees of freedom are the DM-nucleon relative velocity v, the exchanged momentum q, the nucleon spin s N and the DM spin s χ (if different from zero). The scattering amplitude will then be a rotationally invariant function of these variables; invariance under Galilean boosts is ensured by the fact that these vectors are by themselves invariant under Galileo velocity transformations, and translational symmetry is also respected given the absence of a reference frame/point in space. In this regard, a basis of 16 rotationally invariant operators can be constructed with v, q, s N , and s χ [5], which include all possible spin configurations. The scattering amplitude can then be written as a linear combination of these operators, with coefficients that may depend on the momenta only through the q 2 or v 2 scalars ( q • v = q 2 /2µ N by energy conservation, with µ N the DM-nucleon reduced mass). Before introducing these NR operators, however, let us notice that, instead of v, the variable v ⊥ ≡ v -q/2µ N is somehow more suitable to write the amplitude. v ⊥ is Hermitian, in a sense explained e.g. in Ref. [3], while v is not, and moreover one has v ⊥ • q = 0. Following Ref. [3] we will therefore use, in the description of the NR operators, v ⊥ instead of v. The NR operators considered in this work are

O NR 1 = 1 , O NR 3 = i s N • ( q × v ⊥ ) , O NR 4 = s χ • s N , O NR 5 = i s χ • ( q × v ⊥ ) , O NR 6 = ( s χ • q)( s N • q) , O NR 7 = s N • v ⊥ , O NR 8 = s χ • v ⊥ , O NR 9 = i s χ • ( s N × q) , O NR 10 = i s N • q , O NR 11 = i s χ • q , O NR 12 = v ⊥ • ( s χ × s N ) , (1)
where we follow the numbering adopted in Ref. [3,4]. As in [3], we do not consider the full set of independent operators (for instance, as apparent, we do not consider the operator labeled O NR 2 in [3], nor those above the 12 th ); however, as we will see in Sec. 5, the operators listed above are enough to describe the NR limit of many of the relativistic operators often encountered in the literature. The form factor for the operator O NR 12 was obtained from the authors of [3,6].

Given a model for the interaction of DM with the fundamental particles of the SM, we can build the non-relativistic effective Lagrangian describing DM-nucleon interactions as follows. Starting from the fundamental Lagrangian, the matrix element for a scattering process at the nucleon level1 can be expressed as a linear combination of the operators (1):

M N = 12 i=1 c N i (λ, m χ ) O NR i . (2)
The coefficients c N i , where N = p, n can be proton or neutron, are real functions of the parameters of the model, such as couplings, mediator masses and mixing angles, (collectively denoted) λ, the DM mass m χ and the nucleon mass m N . For example, if the scattering between a fermonic DM χ and the nucleon N is described by the (high-energy) scalar operator g N /Λ2 χχ N N , the only non-relativistic operator involved is O NR 1 , and its coefficient is

c N 1 = 4 g N m χ m N /Λ 2 .
The general way to determine the coefficients entering the matrix element (2), starting from high-energy effective operators, will be described explicitly in Sec. 5.

As anticipated above, the c N i can in principle also depend on the exchanged momentum squared q 2 ; in this case we factorize the momentum dependence outside of the coefficients and redefine the c N i as independent from q. The most notable cases of q dependence is featured perhaps in long range interactions, where the exchange of a massless mediator is responsible for the interaction between the DM and the nucleons. The differential cross section displays in this case negative powers of q, thus enhancing the scattering rate at lower exchanged momenta. Assuming that the massless mediator responsible for the interaction is the Standard Model photon, 2 the most relevant cases are a DM with small but nonzero electric charge, electric dipole moment or magnetic dipole moment. As we shall see in more detail in Sec. 5, these interactions all feature a 1/q 2 -dependence [4]. In addition to those in Eq. ( 1), we will therefore consider also the following long range operators:

O lr 1 = 1 q 2 O NR 1 , O lr 5 = 1 q 2 O NR 5 , O lr 6 = 1 q 2 O NR 6 , O lr 11 = 1 q 2 O NR 11 . (3)
According to Eq. ( 55) of [3] we can then write the spin-averaged amplitude squared for scattering off a target nucleus T with mass m T as

|M T | 2 = m 2 T m 2 N 12 i,j=1 N,N =p,n c N i c N j F (N,N ) i,j . (4)
The

F (N,N ) i,j
(v, E R , T ) are the form factors provided in the appendices of [3], and depend critically on the type of scattering nucleus T ; they are also functions of m χ , v and the nuclear recoil energy E R = q 2 /2m T . We can then construct the differential scattering cross section, which reads, in the nonrelativistic case, dσ

T dE R (v, E R ) = 1 32π 1 m 2 χ m T 1 v 2 |M T | 2 . ( 5
)
To write the scattering rate we need to take into account the general case in which the detector is composed of different nuclides (these can be different isotopes of the same specie, as well as different kind of nuclei). We take the numeric abundances of different nuclides used in direct detection searches from Table II of [7], and convert them into mass fractions3 ξ T for each type of target nucleus T , with mass number A T , in the detector. The differential rate for DM scattering off a specific target, expressed in cpd (counts per day) per kilogram per keV, is then

dR T dE R = ξ T m T ρ m χ v min (E R ) d 3 v v f E ( v) dσ T dE R (v, E R ) , (6)
where ρ 0.3 GeV/cm 3 is the DM energy density at the Earth's location and f E ( v) is the DM velocity distribution in the Earth's frame. v min (E R ), the minimum velocity with which a DM particle can scatter off a nucleus with a given recoil energy E R , also depends on the target nucleus via the relation v min = m T E R /2µ 2 T (for elastic scattering), where µ T is the DM-nucleus reduced mass. Using Eq. ( 4) we get the following expression for the differential rate:

dR T dE R = X ξ T 12 i,j=1 N,N =p,n c N i (λ, m χ ) c N j (λ, m χ ) F (N,N ) i,j (E R , T ) , (7)
where we defined

X ≡ ρ m χ 1 32π 1 m 2 χ m 2 N . (8)
and

F (N,N ) i,j (E R , T ) ≡ v min (E R ) d 3 v 1 v f E ( v) F (N,N ) i,j (v, E R , T ) . (9)
In this work we use the customary Maxwell-Boltzmann distribution for the DM velocity. We fix the velocity dispersion to v 0 = 220 km/s and we truncate the distribution at v esc = 544 km/s [8]. We refer to Appendix A for definitions and details. The important thing to notice here is that, up to first order in the non-relativistic expansion of the cross section, one can encounter only two types of velocity dependence, namely dσ T /dE R ∝ v -2 and dσ T /dE R ∝ v 0 (the latter being actually a velocity in-dependence). For these two cases, we use the two integrals I 0 (v min ) and I 1 (v min ) defined in Eq. ( 82) and (83), respectively. Our formalism also allows to treat the case of a linear combination of these two velocity dependences, that one finds e.g. in scattering of DM candidates with magnetic dipole moment. The linearity of dR T /dE R in the form factors is a fundamental ingredient that allows to parametrize a scattering rate in terms of few functions

F (N,N ) i,j
. All the operations we will perform on this quantity will preserve such linearity and will enable us to provide a simple recipe to 'scale' a bound for a specific DM-nucleus interaction to the appropriate bound for another type of interaction.

To actually compare Eq. ( 7) with the experimental rate, we now have to take into account detection efficiency, cuts acceptance and energy resolution of the detector. The set of operations that allows to go from the differential rate (7) to the measured rate is very different from one experiment to the other, and will be discussed in detail later on for each of the experiments considered here (see Sec. 4). For the sake of introducing our results, however, let us sketch now an illustrative procedure. The differential rate must be convolved with the (target-dependent) probability (E ) K T (E R , E ) that a recoil energy E R is measured as E , taking also into account the quenching. (E ) is understood to be the detector's efficiency and acceptance, while K T (E R , E ) reproduces the effect of the energy resolution in every bin, and it is usually assumed to be a Gaussian distribution (with possibly energy-dependent width). In the end we have the differential rate as a function of the detected energy E :

dR dE = T (E ) ∞ 0 dE R K T (E R , E ) dR T dE R . ( 10
)
After convolving with all the experimental effects, the recoil rate of Eq. ( 10) must be averaged over the energy bin of the detector. For each energy bin k of width ∆E k , therefore, we define the number of events predicted by the theory to be

N th k = w k R k ≡ w k ∆E k dE dR dE , (11)
with w k the exposure (expressed in kg • days) and R k the expected rate in the k th energy bin. This quantity can then be directly compared with the detected number of events in the same bin, once the expected number of background events is taken into account. Collecting all the elements in the previous equations we expand Eq. ( 11) and write

N th k = X 12 i,j=1 N,N =p,n c N i (λ, m χ ) c N j (λ, m χ ) F(N,N ) i,j (m χ , k) . (12)
Here F(N,N ) i,j

(m χ , k) is a sort of integrated form factor that encapsulates all the information related to nuclear physics, astrophysics and the detector dependency of the rate. There is one of these factors for each energy bin k of each experiment under consideration, and for each choice of pair of operators i, j and a pair of nucleons N, N . In the approximative determination of the experimental rate sketched above, this would read explicitly

F(N,N ) i,j (m χ , k) = w k T ξ T ∆E k dE (E ) ∞ 0 dE R K T (E R , E ) F (N,N ) i,j (E R , T ) . (13)
While this definition has only illustrative purposes, our analysis will take into detailed account the process of deriving the F(N,N ) i,j

's for each experiment as shown in Sec. 4. At this point one can write the total number of events predicted by the theory as

N th = k N th k = X 12 i,j=1 N,N =p,n c N i (λ, m χ ) c N j (λ, m χ ) F(N,N ) i,j (m χ ) , (14)
where we defined

F(N,N ) i,j (m χ ) ≡ k F(N,N ) i,j
(m χ , k). Due to the linearity of the operations performed to get to Eq. ( 14), the quantity N th inherits from Eq. ( 7) the possibility of being expressed as a linear combination of few ingredients, here the 'integrated' form factors F(N,N ) i,j . The meaning of Eq. 14 is that, once one has computed or is provided with the integrated form factor functions F(N,N ) i,j

, one can straightforwardly derive the expected number of events from possibly any model of DM interactions, whose particle physics is entirely encapsulated in the c N i coefficients, and to immediately compare it with the experimental results. This is the subject of the next Section.

Deriving and rescaling constraints

In the previous Section we have seen how to write down the prediction for the experimental observables (the number of events per bin) in terms of NR operators and the corresponding form factors. In this Section we show how to esplicitly use that formalism in order to derive a bound, from experimental data, on the physics parameters λ entering in the formalism (see e.g. Eq. ( 14)) i.e. ultimately on the DM particle physics. First we derive a bound for a benchmark model. Then we discuss how to translate that result into a bound for an arbitrary choice of operators, i.e. an arbitrary model. Before all this, however, it is useful to review the basics of the statistical method employed to derive constraints, specialized to the case at hand. For experiments that do not see anomalies in their counting rate, the bounds on the free parameter(s) of a given DM model arise by comparing the theoretical number of events with the data, taking also into account the predicted background. As already said, in general direct detection experiments can employ more than one detector (module), or be sensitive to nuclear recoils in different energy bins. For each module or energy bin labelled by k, we can write the total number of expected events as

N k (λ, m χ ) = N th k (λ, m χ ) + N bkg k , (15)
where N bkg k is the expected number of background events. The predicted number of events from DM N th k depends both on the DM mass m χ and on the parameters λ, including e.g. coupling constants and mediator masses. 4To compare the theoretical model with the data, we use a standard Likelihood approach. Since for null result experiments the number of events is very low, data are distributed according to a Poisson distribution. Upon assuming a flat prior, the likelihood of obtaining the set of experimentally observed data N obs given a certain value of the unknown parameter(s) λ is

L( N obs | λ) = k N k (λ, m χ ) N obs k N obs k ! e -N k (λ,mχ) . ( 16
)
It is more convenient to use the log of the Likelihood function, that writes

-2 ln L( N obs | λ) = -2 k N obs k ln N k + 2 k ln N obs k ! + 2 k N k . ( 17
)
Supposing that λ = 0 corresponds to zero rate, i.e. no DM signal, we define the background likelihood L bkg ≡ L( N obs | 0). To compare the DM model with the background-only hypothesis, we adopt the Likelihood ratio Test Statistic

TS(λ, m χ ) = -2 ln L( N obs | λ)/L bkg , (18)
which according to the theory has an approximated χ 2 distribution with a number of degrees of freedom equal to the number of free parameters λ of the DM model. We can then extract, for any given m χ , the maximal value of the parameter(s) λ allowed by the dataset N obs . Explicitly one has

TS(λ, m χ ) = -2 k N obs k ln N th k (λ, m χ ) + N bkg k N bkg k + 2 k N th k (λ, m χ ) . (19)
Bounds on the parameter(s) λ at a given confidence level (CL) can be determined by solving TS(λ, m χ ) = χ 2 CL , where the right-hand side is the χ 2 value corresponding to the desired CL; for instance a 90% CL bound is obtained by choosing χ 2 90% = 2.71, for one single parameter λ. A contour plot in the (m χ , λ) parameter space is the standard way to numerically solve this equation, yielding at the same time a graphical presentation of the result.

Let us now fix a benchmark model of DM interaction with the nuclei. λ B will denote the free parameter of this model, while λ B (m χ ) will denote the bound on it to a given CL, for each value of the DM mass. The total number of expected events for the benchmark model will be denoted as N th,B . Any model can be used as benchmark, and we choose here the simplest possible one:

M p,B = λ B O NR 1 (benchmark), (20)
i.e. c p 1 = λ B while all the other c N i = 0. This results obviously in N th,B = X λ 2 B F(p,p) 1,1 (m χ ). From this it is possible to solve TS(λ B , m χ ) = χ 2 CL for a given experiment to a given CL, thus determining the bound λ B (m χ ). On the website we provide the user with the functions TS(λ B , m χ ) for the four experiments we consider as Mathematica interpolating functions; for a chosen CL, a contour plot in the (m χ , λ B ) plane draws the bound λ B (m χ ) as a function of m χ , as illustrated in Fig. 1. 20) as a function of the DM mass from Xenon100 (solid black), Cdms-Ge (dashed blue), Coupp (dashed-dotted orange), Picasso (dotted red), Lux (long-dashed green) and SuperCdms(short-dashed dark cyan).

We can now scale this benchmark bound to a bound on the parameter λ of another DMnuclei interaction, at the same CL. First one can notice that, once the bound λ B (m χ ) for the benchmark model is given, a bound on λ (as a function of m χ ) for another interaction is trivially obtained by equating

TS(λ, m χ ) = TS(λ B (m χ ), m χ ) , (21)
where the right-hand side is computed within the benchmark model for any value of the DM mass. For an experiment counting the total number of recoil events, the general and unique solution to this equation is N th = N th,B . We can immediately see that this apparently trivial result is actually very powerful: in fact, this translates into

12 i,j=1 N,N =p,n c N i (λ, m χ ) c N j (λ, m χ ) Y (N,N ) i,j (m χ ) = λ B (m χ ) 2 , ( 22
)
where the functions

Y (N,N ) i,j (m χ ) ≡ F(N,N ) i,j (m χ ) F(p,p) 1,1 (m χ ) (23)
are provided on the website as Mathematica interpolating functions. With also the TS(λ B , m χ ) interpolating functions in hand, a bound on λ is promptly obtained by drawing the contour plot of TS(λ B (m χ ), m χ ) at the desired CL, using the expression for λ B (m χ ) in Eq. (22). The Y

(N,N ) i,j
functions are plotted for illustration in figures 2 to 6 below. Notice that, due to the interplay of form factors with different signs for different nuclides, the integrated form factors (and thus the Y's) can be positive or negative depending on the DM mass. In order to accommodate the Y's in logarithmic plots, we therefore show only their absolute value in the figures.

For experiments measuring a spectrum in the recoil energies, instead of the mere total event number, rigorously Eq. ( 21) translates into an equality between the total number of events N th = N th,B (summed over all bins) only if the predicted spectrum is flat, i.e. the expected number of signal events is the same in each bin. However, Eq. ( 22) still provides an excellent approximation to the true bound computed with Eq. ( 19) for the experimental data sets employed here. We verified this explicitly for the operators presented in Sec. 2 for which the maximum difference in number of events in different bins is expected.

Summarizing, on the website we provide:

• The function TS(λ B , m χ ) (as a function of the DM mass m χ over the range 1 GeV → 10 TeV and of λ B over the range 10 -15 → 10 5 ) that allows the user to compute the bound λ CL B (m χ ) at the desired confidence level.

• The functions Y

(N,N ) i,j
(m χ ) for each experiment that we consider (Xenon100, Cdms-Ge, Coupp, Picasso), for each pair of indices (i, j) and each pair (N, N ), as functions of the DM mass m χ over the range 1 GeV → 10 TeV.

• A sample Mathematica notebook which illustrates the usage of the above files.

With these ingredients, Eq. ( 22) allows to set a bound on a parameter λ for any possible interaction type (meaning any possible choice of the c N i ). If λ consists of a set of several parameters, one can analogously draw constraints on one parameter with the others fixed, on combination of parameters (if factorizable) or even obtain multidimensional bounds.

A frequent case is when λ is an overall multiplicative parameter that sets the scale of the cross section (and therefore of the rate). This is, for instance, often the case when λ is a coupling constant or a product of coupling constants, or when λ = Λ -n with Λ a scale of new physics or the mass of the interaction mediator (if one process only dominates). This means that we can set c N i (λ, m χ ) = λ c N i (m χ ) in Eq. ( 2), so that Eq. ( 22) reads

λ B (m χ ) 2 = λ(m χ ) 2 12 i,j=1 N,N =p,n c N i (m χ ) c N j (m χ ) Y (N,N ) i,j (m χ ) . (24)
If the analytic form of the benchmark bound λ B (m χ ) is known, this equation allows to obtain the analytic bound on λ(m χ ). However, the procedure sketched above for obtaining a bound on λ by evaluating the TS(λ B (m χ ), m χ ) with Eq. ( 22) is fully general, and provides fast numerical results if one is able to draw the contour plot of the TS at the desired CL.

In Sec. 6 we will illustrate how to apply Eq. ( 22) in the context of different models and different choices of the parameter λ.

Description of the experiments

In this Section we describe in detail the experiments that we consider, in order to present the specific form that the integrated form factors F take for each one of them. In practice, the following subsections articulate for each experiment the qualitative discussion presented in the second part of Sec. 2. The Lux and SuperCdms experiments are discussed in the (m χ ) for contact interaction and for i = j. The ones shown here are those relative to NR operators for which the effective interaction is mostly spin-independent.

Addenda at page 43 and 45, respectively. The uninterested reader can skip this Section and just consider the F functions (and consequently the Y functions) as black boxes to be plugged into Eq. ( 22).

Xenon100

The Xenon100 detector, located in Gran Sasso National Laboratory in Italy, is a two-phase time projection chamber enclosing 62 kg of target mass. The ratio between the scintillation signal (S 1 ) due to a particle interaction in the liquid xenon, and the subsequent ionization signal (S 2 ), allows for an excellent discrimination of electromagnetic background events. Moreover, the ability to reconstruct the three-dimensional coordinates of each event further enables background reduction by exploiting the self-shielding of liquid xenon in the fiducial volume. Due to their large mass number, Xe nuclei are an excellent target for WIMPs (Weakly Interacting Massive Particles) with spin-independent interactions. However, the unpaired neutron of the 129 Xe and 131 Xe isotopes also makes the experiment sensitive to n-WIMP spin-dependent interactions.

In [START_REF] Aprile | XENON100[END_REF], the collaboration reported the results of the last run, a blind analysis with an exposure w = 34 × 224.6 kg • days which yielded no evidence for DM interactions. The two candidate DM events in the pre-defined nuclear recoil energy range of 6.6 -43.3 keV nr are (m χ ) for contact interaction and for i = j. The ones shown here are those relative to NR operators for which the effective interaction is mostly spin-dependent.

in fact consistent with the background expectation of N bkg = 1.0 ± 0.2 events.

The experiment detects primary (S 1 ) and secondary (S 2 ) scintillation light, that is converted into photoelectrons (PE) by the photomultiplier tubes (PMT). So what is actually measured is a number of electrons instead of the energy of the event. The steps to convert the theoretical event rate into a quantity that is closer to what the experiment actually measures are illustrated in [START_REF] Aprile | XENON100 Coll[END_REF]. The event rate in number of photoelectrons n (summing over the Xe isotopes indexed by T ) is given by

dR dn = ∞ 0 dE R Poiss(n|ν(E R )) T dR T dE R , (25)
where (m χ ) for contact interaction, in the case of interference among NR operators with i = j. Notice that the lines for (i, j) = (4, 5) are not visible as they are superimposed with some of the other ones.

ν(E R ) = E R L eff (E R ) L y S nr /S
single-photoelectron resolution σ PMT = 0.5 PE of the photomultipliers, and the acceptance of the applied cuts (S 1 ) from [START_REF] Aprile | XENON100[END_REF], the resulting S 1 spectrum is

dR dS 1 = (S 1 ) ∞ n=1 Gauss(S 1 |n, √ n σ PMT ) dR dn . (26)
Finally, we compute the total rate as

R = S max 1 S min 1 dS 1 dR dS 1 = S max 1 S min 1 dS 1 (S 1 ) × ∞ n=1 Gauss(S 1 |n, √ n σ PMT ) ∞ 0 dE R Poiss(n|ν(E R )) T dR T dE R , (27)
where S min 1 = 3 PE and S max 1 = 30 PE corresponding to the total energy range 6.6 -43.3 keV nr used in the analysis.

Substituting the expression for dR T /dE R , Eq. ( 7), within Eq. ( 27), and matching with (m χ ) for long range interaction and for i = j.

Eq. ( 14) one finally gets the integrated form factors for the Xenon100 experiment:

F(N,N ) i,j (m χ ) Xenon100 = w T ξ T S max 1 S min 1 dS 1 (S 1 ) × ∞ n=1 Gauss(S 1 |n, √ n σ PMT ) ∞ 0 dE R Poiss(n|ν(E R )) F (N,N ) i,j (E R , T ) . (28)

Cdms-Ge

The Cryogenic Dark Matter Search (Cdms II) experiment is located at the Soudan Underground Laboratory in Minnesota. 19 Germanium and 11 Silicon detectors measure the energy deposited by incident particles in the form of phonons and ionization. The ratio between these two signals provides event-by-event rejection of electron recoils produced by incident electrons and photons. Due to the reduced ionization collection in the external part of the detectors, electron recoils occurring near the surface ('surface events') are more likely to be misidentified as nuclear recoils; however, these events can be discriminated and rejected using phonon pulse timing, further lowering the background from electromagnetic events.

In the final data run of the experiment [13], a blind analysis with an exposure w = 612 kg • days has yielded no significant evidence for DM interactions. Only Ge detectors were used;6 due to its large mass number, Ge is well suitable for direct searches of WIMPs with spin-independent interaction. The collaboration observed two candidate events in the signal region of the detector, against an expected background of misidentified surface events of 0.8 ± 0.1(stat) ± 0.2(syst); in addition, the expected nuclear recoil background counts 0.04 +0.04 -0.03 (stat) events from cosmogenic neutrons and 0.03 -0.06 events from radiogenic neutrons. We consider therefore a total of N bkg = 0.9 background events given by the sum of these independent contributions.

The recoil rate due to DM scatterings in the Cdms-Ge detector can be straightforwardly obtained by integrating the theoretical differential rate over the energy window of the detector, taking into account the nuclear recoil efficiency (E R ) [13] and summing over germanium's isotopes:

R = E max E min dE R (E R ) T dR T dE R , (29)
with E min = 10 keV nr and E max = 100 keV nr . In this case the integrated form factor is given by

F(N,N ) i,j (m χ ) Cdms-Ge = w T ξ T E max E min dE R (E R ) F (N,N ) i,j (E R , T ) . ( 30
)

Coupp

The Chicagoland Observatory for Underground Particle Physics (Coupp) is a 4 kg CF 3 I bubble chamber operated at SNOLAB in Ontario. Because of its unpaired proton, fluorine provides excellent sensitivity to p-WIMP spin-dependent interactions, while iodine enhances the sensitivity to spin-independent interactions due to its large mass number. Bubble nucleations, triggered by particles entering the liquid in the superheated phase, are recorded both photographically and by pressure and acoustic measurements. With appropriate chamber pressure and temperature, electron recoils due to the abundant gammaray and beta-decay backgrounds do not nucleate bubbles; therefore, the main background is constituted by neutrons and alpha decays. Alpha events are discriminated by an acoustic cut, while photographic and pressure measurements are used to identify and reject events occurring outside the fiducial volume as well as events with multiple bubble nucleation.

Even though the experiment is not able to tell the scattering energy on an event-byevent basis, the superheated fluid can be operated at different temperatures, corresponding to different energy thresholds for bubble nucleation. During the last run [16], data were collected in three contiguous data collection periods at temperatures of 39.0 • C, 36.2 • C and 33.5 • C, corresponding to thresholds of 7.8, 11.0 and 15.5 keV nr , respectively. Unfortunately there is still no precise understanding of how the probability of nucleating a bubble as a function of the energy deposition depends on the threshold. To estimate this uncertainty the collaboration uses two different parametrizations for the nucleation efficiency, fitting the data equally well, thus drawing the uncertainty band visible in figures 6 and 7 of [16]. 7In [16], data obtained for an effective exposure to single recoil events of 437.4 kg • days (taking into account the 79.1% detection efficiency) were presented. Twenty single nuclear recoil events passing all the analysis cuts were observed over the three energy bins. However, a more accurate analysis of the data led to think that not all these events are genuine single nuclear recoils. Upon implementation of a further time isolation cut, the number of signal events is reduced to thirteen as reported in Table 2. Due to uncertainties in the neutron background estimation, the collaboration does not attempt any background subtraction and instead treats all the twenty originally observed nuclear recoils as DM events. However, our Test Statistic cannot be applied to the case of zero background, and therefore we adopt the background estimates made by the collaboration, also reported in Table 2. We enforce the time isolation cut made by the collaboration, for this yields a bound closer to the experimental one.

The rate of DM events in the k-th energy bin, identified by its low energy threshold E thr k , is determined by convolving the theoretical differential rate with the target-dependent nucleation efficiency P T (E R , E thr k ):

R k = ∞ E thr k dE R T =C,F,I P T (E R , E thr k ) dR T dE R . ( 31
)
Scattering of iodine in CF 3 I is known to have a good efficiency, which can be thus taken to be 1 above threshold, [16]. We only consider the events passing the time isolation cut.

P I (E R , E thr k ) = θ(E R -E thr k ).
is uncertain, therefore the collaboration models it as

P C,F (E R , E thr k ) = 1 -exp a C,F 1 - E R E thr k , (32)
where a fit to the data yields a C,F = 0.15. The parameter a determines the rise in the nucleation efficiency with energy (higher a means steeper efficiency). This low value of a F therefore greatly decreases the sensitivity of the detector, especially at low energies.

The other parametrization of P(E R , E thr k ) proposed by the collaboration, consisting of a step function centered at the threshold, yields probably tighter constraints on the DM parameter space. Conservatively, we prefer to use formula (32), that is also the one used by the Picasso collaboration (which utilizes fluorine nuclei as well).

The integrated form factors will be then

F(N,N ) i,j (m χ ) Coupp = k w k T =C,F,I ξ T ∞ E thr k dE R P T (E R , E thr k ) F (N,N ) i,j (E R , T ) . (33)
Actually, Ref. [3] does not provide the form factors for carbon, and therefore we omit its contribution in Eq. (33). We note however that this is an excellent approximation, in that spin-dependent interactions are mainly probed by fluorine, while spin-independent interactions are more sensitive to heavier nuclei such as iodine and fluorine.

Picasso

The Picasso experiment, located at SNOLAB, searches for Dark Matter using superheated liquid droplets, a variant of the bubble chamber technique. The abundance of fluorine in the target liquid C 4 F 10 yields an excellent sensitivity to spin-dependent WIMP interactions with protons, especially at low DM masses due to small fluorine mass and the very low recoil detection threshold. Similarly to Coupp, particles which produce low ionization densities such as cosmic muons, gamma rays and beta radiation, become detectable only at high temperatures. They are therefore well separated from strongly ionizing WIMP induced recoils, thus allowing efficient suppression of such backgrounds. Alpha particles feature instead a flat spectrum up to high energy thresholds and constitute therefore the main background.

Of the 32 detector modules accommodated by Picasso, only 10 were used in the latest data analysis [17]. By varying the temperature, the collaboration constructed the spectrum of the particle induced energy depositions exploiting the strict correspondence of detector temperature to threshold energy. After correcting for cut acceptances and dead time, the rates recorded by the modules at each temperature are normalized with respect to the active fluorine mass; as for Coupp, the contribution of carbon is negligible and we disregard it in the following (the collaborations estimates it as 10% of the total spin-independent cross section [17], while it plays presumably no role in spin-dependent interactions). The counting rates of all detector modules result to be flat in the considered range of temperatures 28 • C < T < 48 • C, corresponding to low threshold energies between 1.7 and 55 keV, in contrast with a decreasing expected DM signal. The flat rate can be explained by the presence of alpha emitters in the droplets, whose origin is still uncertain. For this reason, the collaboration does not try to predict the background, rather it fits it to a detector-dependent constant using the data (different modules display different levels of alpha contamination, due to differences in the fabrication processes). Since the collaboration does not disclose the data, but only provides averaged quantities, we proceed as follows. We 'virtually' subtract the average background R j /8 from the data points R j (T k ) of the relative module (indicated by the index j) for the eight temperature values T k probed by the experiment. These rate fluctuations can then be directly compared, for a certain temperature T k , among different modules: averaging over different detectors we then get

∆R(T k ) ≡ j w j R j (T k ) -R j /8 j w j , (34)
where w j is the exposure of the j th module. What one should have used here are actually the exposures relative to the temperature T k , reflecting the time spent by each module at that temperature; since we do not have this information, we assume these exposures to be fixed fractions of the w j 's, independent on the detector. The ∆R(T k ) are depicted in figure 5 of [17]; the R j and the w j are given respectively in tables 3 and 1 of [17]. All these quantities are reported here in Tab. 3. From Eq. (34) one recovers the following formula for the total rate R obs k measured by the experiment for a given energy threshold E thr k (corresponding to a temperature T k ):

R obs k ≡ 1 w j w j R j (E thr k ) = ∆R(E thr k ) + R , (35)
where w ≡ j w j = 114.3 kg • days and R ≡ 1 8w j w j R j = 14.9 cpd/kg as from Tab. 3. The total rates R obs k can be computed in any given model of DM interactions, adding to the constant background R the DM contribution in the k-th energy bin R th k , and can then be compared with the experimental outcome Eq. ( 35) to a chosen confidence level. In analogy with Coupp,

R th k = ∞ E thr k dE R P F (E R , E thr k ) dR F dE R , (36)
where dR F /dE R is the DM scattering rate off fluorine nuclei. P F (E R , E thr k ) is defined in Eq. ( 32), whereas the Picasso collaboration chooses a F = 5 while allowing it to vary within ±2.5 around this central value.

Notice that, since the collaboration does not provide the exposures w k for the single energy bins, we are not able to determine the total number of events from the rate, which we used in Eq. ( 19) to explain the statistical procedure. However, the same procedure can be applied to rates rather than number of events. In this regard we can formally define the numbers

N th k , N bkg k and N obs k as N th k ≡ wR th k , N bkg k ≡ wR , N obs k ≡ wR obs k , (37)
where w is a dummy exposure that can take any value. The integrated form factors will then be

F(N,N ) i,j (m χ ) Picasso = w ξ F k ∞ E thr k dE R P F (E R , E thr k ) F (N,N ) i,j (E R , F) . (38)
The dependence on w will cancel in the functions Y (N,N ) i,j (m χ ) by their definition in Eq. ( 23).

5 A dictionary: from quark/gluon-level relativistic effective operators to NR operators

In this Section we review how to decompose the high-energy effective operators commonly used in model building into the non-relativistic bricks O NR k of Eq. ( 1). The starting point are the interaction operators of the DM with quarks and gluons; we review how the step up to the interaction at the nucleon level is performed and then how these are expressed in terms of NR operators, to which the results of Sec. 2 can be applied straightforwardly.

Incidentally, we will also see explicitly that different high-energy effective operators might have the same non-relativistic form, corresponding thus to the same O NR k . This can also give rise to interference effects that might significantly lower or enhance the scattering cross section, thus generating phenomenologies that are usually not taken into account by the effective operators analyses when they consider one operator at a time.

Notice that, while long-distance QCD effects induce energy-dependent corrections to the scattering amplitude, we will only present the matching from quark and gluon level to the nucleon level at lowest order. Next to leading order effects, including two-nucleon interactions [18], have been studied in the case of scalar interactions in Ref. [19], and for spin-dependent (axial-vector) interactions in Ref. [20,21].

Effective operators for fermion Dark Matter

At dimension six, the effective operators one can construct with a Dirac neutral DM field χ and quark fields q are O q 1 = χχ qq , O q 2 = χ iγ 5 χ qq , O q 3 = χχ q iγ 5 q , O q 4 = χ iγ 5 χ q iγ 5 q , O q 5 = χγ µ χ qγ µ q , O q 6 = χγ µ γ 5 χ qγ µ q , O q 7 = χγ µ χ qγ µ γ 5 q , O q 8 = χγ µ γ 5 χ qγ µ γ 5 q , O q 9 = χ σ µν χ q σ µν q , O q 10 = χ i σ µν γ 5 χ q σ µν q ,

where we do not take into account here flavor-violating interactions. Notice that the operators χ σ µν χ q i σ µν γ 5 q , χ i σ µν γ 5 χ q i σ µν γ 5 q (40)

are equal to O q 10 and -O q 9 , respectively, by virtue of the identity i σ µν γ 5 = -1 2 ε µνρτ σ ρτ . For a Majorana DM, only the bilinears χχ, χγ 5 χ and χγ µ γ 5 χ are non-zero.

Gauge-invariant interaction operators with gluons arise at dimension seven, and are

O g 1 = α s 12π χχ G a µν G a µν , O g 2 = α s 12π χ iγ 5 χ G a µν G a µν , O g 3 = α s 8π χχ G a µν Ga µν , O g 4 = α s 8π χ iγ 5 χ G a µν Ga µν , (41)
where Ga µν ≡ ε µνρσ G a ρσ , and the numerical overall factors have been chosen for later convenience.

The effective Lagrangian at the quark-gluon level is

L eff = 10 k=1 q c q k O q k + 4 k=1 c g k O g k , (42)
where the c q k and c g k are real dimensionful coefficients: 8 c q k will have dimensions of [mass] -2 and c g k of [mass] -3 . As briefly reviewed in Appendix B, these operators induce an effective Lagrangian at the nucleon level

L eff = 10 k=1 N =p,n c N k O N k , (43)
where the

O N k (N = p, n) are O N 1 = χχ N N , O N 2 = χ iγ 5 χ N N , O N 3 = χχ N iγ 5 N , O N 4 = χ iγ 5 χ N iγ 5 N , O N 5 = χγ µ χ N γ µ N , O N 6 = χγ µ γ 5 χ N γ µ N , O N 7 = χγ µ χ N γ µ γ 5 N , O N 8 = χγ µ γ 5 χ N γ µ γ 5 N , O N 9 = χ σ µν χ N σ µν N , O N 10 = χ i σ µν γ 5 χ N σ µν N , (44)
8 The c q k coefficients are not to be confused with c N k coefficients defined below nor with the c N k coefficients introduced in Eq. ( 2): the first ones are the coefficients in the expansion in terms of quark/gluon level effective operators, the second ones of the expansion in nucleon level operators, the last ones in NR operators.

and we denoted with N the nucleon field. Notice that the gluon operators contribute to the scalar operators

O N 1 , O N 2 , O N 3 and O N 4 . The couplings are c N 1,2 = q=u,d,s c q 1,2 m N m q f (N ) T q + 2 27 f (N ) T G q=c,b,t c q 1,2 m N m q -c g 1,2 m N , (45a)
c N 3,4 = q=u,d,s m N m q (c q 3,4 -C 3,4 ) + c g 3,4 m ∆ (N ) q , (45b)
c N 7,8 = q c q 7,8 ∆ (N ) q , ( (45c)
) c N 9,10 = q c q 9,10 δ (N ) q , ( 45d
) 45e
where C 3,4 ≡ q c q 3,4 m/m q with m ≡ (1/m u + 1/m d + 1/m s ) -1 ; the factors f T q , f T G , ∆

(N ) q and δ (N ) q are given in Appendix B, to which the reader should refer for the derivation of the equations (45) and where reference to the relevant literature is provided. For the quark scalar couplings c q 1 to c q 4 it is usually assumed c q k ∝ m q , as it would be the case for a DMquark interaction mediated by a higgs-like particle coupling to the quark masses (we will go back to this point in Sec. 6.1).

To compute the DM-nucleus amplitude, we have now to coherently sum up the interaction amplitude over all nucleons in the nucleus, taking also into account the physics of the bound state. The form factor formalism of Ref. [3] allows to make this for any type of interaction, only we have first to evaluate the matrix elements out χ, N |O N k |χ, N in and to express them in terms of the non-relativistic operators O NR i introduced in Eq. (1). To do so, we can expand the solution of the Dirac equation in its non-relativistic limit: in the Weyl or chiral representation for the spinors,

u s (p) = √ p µ σ µ ξ s √ p µ σµ ξ s = 1 2(p 0 + m) (p µ σ µ + m) ξ s (p µ σµ + m) ξ s = 1 √ 4m (2m -p • σ) ξ s (2m + p • σ) ξ s + O( p 2 ) ( 46
)
where σ µ = (1, σ), σµ = (1, -σ) and we approximated p µ = (m, p) + O( p 2 ) in the nonrelativistic limit. In this limit we can study the velocity, momentum and spin dependence of the fermion bilinears, when both fermions are on-shell. Up to and including the first order in the three-momenta,

ū(p )u(p) 2m , (47a) ū(p )i γ 5 u(p) 2i q • s , (47b) ū(p )γ µ u(p) 2m P + 2i q × s , (47c) ū
(p )γ µ γ 5 u(p) 2 P • s 4m s , (47d)
ū(p )σ µν u(p)

0 i q -2 P × s -i q + 2 P × s 4m ε ijk s k , (47e) ū
(p )i σ µν γ 5 u(p) 0 -4m s 4m s i ε ijk q k -2P i s j + 2P j s i , (47f)
where q = p -p is the exchanged momentum, and P = p + p . The spin operator is defined as s ≡ ξ † σ 2 ξ, where in its absence a ξ † ξ is understood. Finally, when contracting fermionic DM and nucleon bilinears, 9 the following expressions can be derived for the (matrix elements of the) O N k , at leading order in the non-relativistic expansion:

O

N 1 = O N 5 = 4m χ m N O NR 1 , O N 2 = -4m N O NR 11 , O N 3 = 4m χ O NR 10 , O N 4 = 4O NR 6 , O N 6 = 8m χ +m N O NR 8 + O NR 9 , O N 7 = 8m N -m χ O NR 7 + O NR 9 , O N 8 = - 1 2 O N 9 = -16 m χ m N O NR 4 , O N 10 = 8 m χ O NR 11 -m N O NR 10 -4m χ m N O NR 12 . (49)
Here we took q to be the momentum transferred by the DM to the nucleus. With these substitutions, the effective Lagrangian (43) gives origin to the nucleonic matrix element in the form of Eq. ( 2), and the results of Sec. 2 can be applied straightforwardly. 10 From Eq. ( 49) one can now see clearly that O N 1 and O N 5 correspond to the same nonrelativistic operator, and so do O N 8 and O N 9 . They are therefore indistinguishable by direct detection experiments alone. This also means that a bound computed on one of these operators is identical to the bound computed on the other, a difference arising only if different coefficients are chosen for the two in the Lagrangian. Furthermore, if a model 9 The following expression can be useful in expressing the result as function of only q and v ⊥ :

P χ m χ - P N m N = 2 v ⊥ , (48)
where P χ ( P N ) is the sum of the initial and final DM (nucleon) momenta. 10 Notice that, while in Ref. [3,4] the matrix element is defined as in . . . out , we use out . . . in instead. For this reason we expect to get a minus sign for each power of q occurring in the NR operators. However, we use an opposite convention for q so that the final result will be the same.

features both operators, strong cancellations or enhancements of the scattering cross section might arise due to their quantum interference. The simplest case in which this could happen is a theory where the DM interacts with quarks via a scalar and a vector exchange, thus producing both O N 1 and O N 5 .

Effective operators for scalar Dark Matter

For a scalar DM, the following contact operators for interaction with quarks are possible:

O q 1 = φ * φ qq , O q 2 = φ * φ q iγ 5 q , O q 3 = i (φ * ← → ∂ µ φ) qγ µ q , O q 4 = i (φ * ← → ∂ µ φ) qγ µ γ 5 q . ( 50
)
The operators in the first line are dimension five, while the others have dimension six and vanish identically for a real field φ. The operator ∂ µ (φ * φ) qγ µ q does not contribute to processes where the quarks are external because it is proportional (after integration by parts) to the divergence of the conserved current qγ µ q, that vanishes by virtue of the equations of motion. The operator ∂ µ (φ * φ) qγ µ γ 5 q, instead, reduces to 2m q φ * φ q iγ 5 q upon integration by parts and application of the equations of motion; it is therefore equivalent to O q 2 , the only difference being the quark mass factor that can be however absorbed in the arbitrary Lagrangian coefficient. At dimension six also the terms

φ * φ q i ← → / D q , φ * φ q i ← → / D γ 5 q (51)
exist, where

← → D µ = 1 2 ← → ∂ µ -ieQ q A µ -ig s T a
q g a µ is the ('hermitianized') quark covariant derivative. We do not take into account these operators here.

The most relevant gauge-invariant interaction operators with gluons have dimension six and are

O g 1 = α s 12π φ * φ G a µν G a µν , O g 2 = α s 8π φ * φ G a µν Ga µν , (52)
where the numerical factors have been chosen for later convenience.

The effective Lagrangian at the quark-gluon level is therefore

L eff = 4 k=1 q c q k O q k + k=1,2 c g k O g k , (53)
where the c q n and c g n are real dimensionful coefficients: c q k for k = 1, 2 will have dimension of [mass] -1 , c q k for k = 3, 4 dimension of [mass] -2 and c g k of [mass] -3 . This induces the effective Lagrangian at the nucleon level

L eff = 4 k=1 N =p,n c N k O N k , (54)
where the O N k (N = p, n) are

O N 1 = φ * φ N N , O N 2 = φ * φ N iγ 5 N , O N 3 = i (φ * ← → ∂ µ φ) N γ µ N , O N 4 = i (φ * ← → ∂ µ φ) N γ µ γ 5 N , (55)
and the couplings are, in analogy with the fermion DM case,

c N 1 = q=u,d,s c q 1 m N m q f (N ) T q + 2 27 f (N ) T G q=c,b,t c q 1 m N m q -c g 1 m N , (56a)
c N 2 = q=u,d,s m N m q [(c q 2 -C 2 ) + c g 2 m] ∆ (N ) q , (56b)
c p 3 = 2 c u 3 + c d 3 , c n 3 = c u 3 + 2 c d 3 , (56c)
c N 4 = q c q 4 ∆ (N ) q . ( 56d
)
We again point to Appendix B for the derivation of these formulas and for reference to the relevant literature.

As we did in Sec. 5.1 we can now match the (matrix element of the) relativistic operators (55) to the non-relativistic ones (1), at the leading order in the non-relativistic expansion, obtaining

O N 1 = 2m N O NR 1 , O N 2 = 2O NR 10 , O N 3 = 4m χ m N O NR 1 , O N 4 = -8m χ m N O NR 7 . (57)

Long range interactions

While contact interactions as the ones described above are either independent on the exchanged momentum or momentum-suppressed, long range interactions are enhanced at small momentum transfer. These arise from the exchange of a massless mediator, whose propagator is responsible for the enhancement. Since these interactions feature quite different spectra with respect to the ones seen so far, it is worth discussing also this case.

For definiteness we take into account here interactions of the DM with the Standard Model photon, that have already been studied in the context of DM direct searches (see e.g. [22][23][24][25][26][27][28][29] and references therein). Notice that interactions with gluons are not of long range type, as at the low energies relevant for direct detection the interaction is effectively with the entire nucleon (see above). The most relevant gauge-invariant interactions a fermionic Dirac DM χ can have with photons are

O C = Q χ e χγ µ χ A µ , (millicharged DM) ( 58
)
O M = µ χ 2 χ σ µν χ F µν , (anomalous DM magnetic moment) ( 59
)
O E = d χ 2 i χ σ µν γ 5 χ F µν . (DM electric dipole moment) ( 60
)
The dimensionless constant Q χ is the DM electric charge in units of e, whereas the parameters µ χ and d χ , with dimension [mass] -1 , are the magnetic and electric dipole moments (in units of e • cm) of χ, respectively. All these interactions vanish identically for Majorana DM particles; the lowest electromagnetic moment allowed for a Majorana particle is the anapole moment, χ γ µ γ 5 χ ∂ ν F µν , but this does not lead to a long range interaction and therefore we will not consider it here [27]. These interactions generate the following non-relativistic nucleonic matrix element (2), when the photon is to mediate the interaction with nucleons N = p, n (see also [4]):

M N C = 4 e 2 Q χ Q N m χ m N O lr 1 , (61)
M N M = 2 eµ χ Q N m N O NR 1 + 4 Q N m χ m N O lr 5 + 2 g N m χ O NR 4 -O lr 6 , (62)
M N E = 8 ed χ Q N m χ m N O lr 11 , (63)
where Q p = 1, Q n = 0 are the nucleons electric charge while g p = 5.59, g n = -3.83 are their g-factors. The operators O lr k ≡ O NR k /q 2 have been defined in Sec. 2. Consequently, we define appropriate form factors for the long range case, that take into account this dependence on q; from these, the relative rescaling functions Y lr(N,N ) i,j (m χ ) are then derived. In this way one is not concerned with any momentum dependence of the coefficients in (2); for example, for millicharged DM (61) we have c p 1 = 4 Q χ m χ m N and 0 for all the other coefficients.

For a scalar DM φ, the only relevant interaction operator with the photon displaying long range behavior is the one coming from the kinetic term, in case the DM has a small electric charge:

M N C = Q φ e i (φ * ← → ∂ µ φ)A µ . ( 64
)
For this operator, the effective interaction Lagrangian with nucleons is the same as Eq. ( 61), with the obvious substitution

Q χ → Q φ .
Other possible interaction operators of φ with photons, arising at dimension six, are

i (φ * ← → ∂ µ φ) ∂ ν F µν , φ * φ F µν F µν , φ * φ F µν Fµν , (65)
with F µν ≡ ε µνρσ F ρσ (see e.g. Ref. [30]). The terms

∂ µ ∂ ν (φ * φ) (∼)
F µν and i (φ * ← → ∂ µ φ) ∂ ν F µν , vanish identically because of the antisymmetricity of F µν and ε µνρσ . However, these interactions are not of long range type (see e.g. [31,32]) and therefore we don't include them here.

Summary and explicit examples

In this last Section we summarize the steps one can take without even reading the rest of the paper to be able to use our results and start setting bounds on his/her favorite DM candidate.

1a. Compute the amplitude for the relevant scattering process in your favorite DM model.

At the end of the day it should consist of a linear combination of quark-and gluonlevel high-energy operators. The combination will include, but will not generally be limited to, the operators of Eq. ( 39) and ( 41), for fermionic DM, or Eq. ( 50) and (52) for scalar DM (i.e. the amplitude should be in the form of Eq. ( 42) or ( 53)). The amplitude can also contain the long range operators of Eq. ( 58) to (60).

The coefficients of the operators (the c q,g k of the contact effective operators or Q χ , µ χ , d χ ) encapsulate the particle physics details of the model, including the unknown parameter(s) λ on which you wish to derive a bound. 1b. Dress up the quark-and gluon-level operators to the nucleon level in the standard way. For the effective operators, we review this process in Sec. 5 and Appendix B. One obtains an effective Lagrangian in the form of Eq. ( 43) or (54). The coefficients c N k now encapsulate the particle physics (including λ) and the internal nucleon dynamics of quarks and gluons.

1c. Reduce to the NR limit as illustrated in Sec. 5, in particular using the dictionary provided by Eq. ( 49), ( 57) or ( 61)-( 63) and (64). Matching the terms obtained in this expansion with the NR operators in Eq. ( 1) will allow to identify the c N i (λ, m χ ) as the respective coefficients, as in Eq. ( 2). 2a. For each experiment, plug in the c N i (λ, m χ ) in Eq. ( 22), using the interpolated functions Y (N,N ) i,j (m χ ) that we provide on the website. This provides you with an expression for λ B , implicit in the unknown parameter λ. We also present a few explicit examples. The first one is the simplest possible case one can imagine: considering one effective high-energy operator at a time. The second example aims to make contact between our formalism and the usual bounds in terms of spin-independent (SI) and spin-dependent (SD) cross sections. These, together with bounds on millicharged DM and DM with magnetic moment interaction, are also illustrated in the sample Mathematica notebook provided on the website. The third example is a more involved, but well defined, framework inspired by Minimal Dark Matter, in which the DM-ordinary matter scattering involves a combination of several (high-energy) operators with non trivial coefficients.

Bounds on individual relativistic effective operators

The simplest application of our formalism is the idealized case in which one considers one high-energy contact operator at a time. This is what actually is routinely done in papers employing the effective operator formalism for DM searches (see e.g. Ref. [33][34][35][36][37]). We use this simple case to explicitate the steps outlined above in a concrete example. For instance, if we study the case of the k = 1 quark operator for fermionic DM

L eff ≡ m q Λ 3 O q 1 =⇒ c q 1 = m q Λ 3 , c q k =1 = c g k ≡ 0 , ( 66
)
where Λ plays the role of the sole unknown parameter and, as usual, it represents the scale of the ultraviolet physics which has been integrated out. A possible overall numerical factor is conventionally assumed to be 1 or equivalently to be absorbed in the unknown scale Λ. Next, as dictated by Eq. (45a) Thus, following Eq. ( 49)

c N 1 = u,d,s f (N ) T q + 2 9 f (N ) T G m N Λ 3 . ( 67
)
c N 1 = 4m χ u,d,s f (N ) T q + 2 9 f (N ) T G m 2 N Λ 3 . ( 68
)
It is now straightforward to apply Eq. ( 22) and obtain

λ 2 B = 16 m 2 χ Λ 6 N,N =p,n m 2 N u,d,s f (N ) T q + 2 9 f (N ) T G m 2 N u,d,s f (N ) T q + 2 9 f (N ) T G Y (N,N ) 1,1 (m χ ) .
(69) Now by explicitly plugging this expression in the TS(λ B , m χ ) function for each experiment, one has the TS(Λ, m χ ) function which is needed to impose a bound on Λ. In Fig. 7, top left panel, we show the result using the values of the f (N ) T q quoted in Ref. [38] (see Table 4 in Appendix B).

We repeat the exercise for all the high-energy contact operators in Eq. ( 39) and Eq. ( 50), showing the results in Fig. 7. For O q 1,2,3,4 for fermionic DM we assume coefficients proportional to m q /Λ 3 , like in (66) as it is customary (since this applies in the case in which the mediator particle couples to quarks like the higgs boson does). Analogously, we assume coefficients proportional to m q /Λ 2 for scalar DM operators O q 1,2 .

'Usual' SI and SD interactions

As another example, we discuss here how to recover, within our formalism, the usual spinindependent and spin-dependent bounds presented by the experiments. As we saw in Sec. 2, there are different non-relativistic interactions that lead to a detector response which is independent on the nucleus spin or dependent on the nucleus spin. When experiments present results in terms of 'spin-independent' or 'spin-dependent' scattering cross section, however, they implicitly assume the lowest order interactions in the non-relativistic expansion, namely O NR , respectively, and could be used for the SI and SD analyses in place of O N 1 and O N 8 , the former ones have been historically considered first (probably as they are the operators relevant in describing elastic scattering of the supersymmetric neutralino). The effective DM-nucleon Lagrangians are

L N SI = λ SI O N 1 , (70)
L N SD = λ SD O N 8 . (71)
From this, using Eq. ( 49) to convert the operators to their non-relativistic limit expression, one immediately obtains

c N 1 (λ SI , m χ ) = 4 λ SI m χ m N , (72)
c N 4 (λ SD , m χ ) = -16 λ SD m χ m N . (73)
At this point a bound can be set on λ SI and λ SD using Eq. (22). Rather than presenting this kind of bounds, let us make contact with the usual physical cross sections. For SI scattering, DM-proton and DM-neutron couplings are customarily assumed to be equal, and so are the two cross sections (neglecting the proton-neutron mass difference). These read

σ SI p ≡ σ SI n = λ 2 SI µ 2 N /π
, where µ N is the DM-nucleon reduced mass. Therefore, to determine a bound on σ p from the usual SI interaction within the framework proposed in this work, one just needs to use λ 2 SI = πσ SI p /µ 2 N in Eq. ( 22), together with the value of c N 1 (m χ ) given in Eq. ( 72) and the Y (N,N )

(m χ ) functions provided in the website. That is, use

λ B (m χ ) 2 = σ SI p • 16πm 2 χ m 2 N µ 2 N Y (p,p) 1,1 (m χ ) + 2 Y (p,n) 1,1 (m χ ) + Y (n,n) 1,1 (m χ ) (74)
in the Mathematica interpolated function TS(λ B , m χ ) provided on the website. Notice that, exactly as for the original form factors of Ref. [3], one has Y

(N,N ) i,j = Y (N ,N ) j,i
.

In the SD case, contrarily to what happens for the SI interaction, the DM-p and DM-n cross sections are usually considered separately, in the assumption that only protons or neutrons contribute, respectively. The cross section on free nucleons is given by σ SD N = 3λ 2 SD µ 2 N /π, and as in the previous case we can invert this equation to express the bound on λ SD provided by Eq. ( 22) in terms of σ SD N , where N is either proton or neutron. Here we will use c N 4 (m χ ) as in Eq. ( 73), as well as the Mathematica interpolated functions Y (N,N )

(m χ ). Explicitly, using

λ B (m χ ) 2 = σ SD N • 256 3 πm 2 χ m 2 N µ 2 N Y (N,N ) 4,4 (m χ ) (75)
into TS(λ B , m χ ) will provide the desired bound for the nucleon N , either proton or neutron.

Our results for both SI and SD cross sections are plotted in Fig. 8, and they reproduce remarkably well the results given by the experimental collaborations in the references cited in Table 1. The only case where we find a notable difference of our bounds with the ones published by the experimental collaborations is for the spin-dependent interaction with protons in xenon. The reason for this difference has to be found in the uncertainty to which the xenon form factor for this interaction is known; in fact, as shown in Fig. 1 of Ref. [START_REF] Aprile | XENON100 Coll[END_REF], results found by different collaborations disagree even by orders of magnitude.

While the agreement is good even in the small DM mass region (m χ 10 GeV), we warn that that is the most critical one, where dedicated analyses have been performed by the experimental collaborations and by independent groups (see e.g. Ref. [39][40][41][42][43][44]).

DM DM DM

W q q h DM DM DM W W q q DM DM DM W W q q
Figure 9: One loop DM-quark scattering for fermionic MDM with Y = 0. Figure adapted from [45].

DM DM DM W Q g g h DM DM DM W W Q g g DM DM DM W W Q g g

Minimal Dark Matter-inspired model

As a final example, we consider a more involved model in which more than one operator contributes at the same time.

The Minimal Dark Matter (MDM) [45][46][47][48] construction consists in adding to the Standard Model (SM) the minimal amount of new physics (just one extra EW multiplet X) and searching for the minimal assignments of its quantum numbers (spin, isospin and hypercharge) that make it a good Dark Matter candidate without ruining the positive features of the SM. The theory univocally selects a fermionic 5-plet of SU(2) L with hypercharge Y = 0 as the best candidate [48].

Of course, the 'pure' MDM model as such is not of interest for our current analysis since it has no free parameters and therefore no parameter on which we can compute the bounds. We will therefore consider a modified version in which we assume that DM has a reduced SU(2) coupling g = g, where g is the ordinary SU(2) L gauge coupling. We will thus be able to compute the bounds from different experiments on the quantity . Moreover, while in the pure MDM model the mass of X is determined by the relic density requirement, here we leave it as a free parameter.

The MDM quintuplet scatters elastically on nuclei via the one loop diagrams in Fig. 9 and via the two loops diagrams in Fig. 10. Upon a careful computation 11 one arrives at determining the high-energy effective scattering Lagrangian

L eff = 3 32 π 2 g 4 M W q 1 m 2 h m q XX qq - 2 3m X Xγ µ γ 5 X qγ µ γ 5 q - 1 16 π 2 g 4 M 3 W m X X(i∂ µ γ ν )X 1 2 q i(D µ γ ν + D ν γ µ - 1 2 g µν / D)q - 1 16 π 2 g 4 M W Q κ Q m 2 h + 1 M 2 W α s 8 π XX G a µν G a µν , (76)
where Q = c, b, t in the last line, with κ c = 1.32, κ b = 1.19, κ t = 1 [49]. As apparent, the effective Lagrangian contains three dominant operators among those listed in Eq. ( 39)-( 41): O q 1 , O q 8 and O g 1 . It also contains the term with a twist-2 quark operator (let us denote the full operator as

O q A = X(i∂ µ γ ν )X O q,twist-2 µν
), which is not among the standard ones listed above and therefore has to be reduced to its NR limit explicitly. The coefficients explicitly read

c q 1 = 2 3 32 π g 4 M W m q 1 m 2 h , (77a)
c q 8 = -2 1 16 π g 4 M W m X , (77b)
c q A = -2 1 16 π g 4 M 3 W m X , (77c)
c g 1 = 2 1 16 π g 4 M W 3.51 m 2 h + 1 M 2 W . (77d)
O q A contributes O N 1 [49] and therefore, at the nucleon level, the effective operators O N 1 and O N 8 have coefficients (see Eq. ( 45))

c N 1 = 2 3 32 π g 4 M W m N 1 m 2 h q=u,d,s f (N ) T q + 2 27 f (N ) T G (1 -2 3 Q κ Q ) - 1 M 2 W 4 81 f (N ) T G + 1 2 q=u,d,s,c,b (q (N ) (2) + q(N) (2)) , (78a)
c N 8 = -2 1 16 π g 4 M W m X q ∆ (N ) q . ( 78b
)
in Ref. [49] (for the two loops diagrams). We always work in the limit of M W , m t m X , where M W is the W mass and m t is the top mass, relevant for the two loop diagrams. One can check that the term containing the twist-2 operator gives rise to the term proportional to 1/M 3 W in the original computation of [45], however with an opposite sign and with an o(1) difference in the coefficient of which we are not able to pin down the origin: we adopt the more recent result of [49]). The different sign causes a cancellation between the term in 1/M W m 2 h and the term in 1/M 3 W , which significantly suppresses the cross section [49]. Also, for the coefficient of the O q 8 operator we have a discrepancy of a factor of 2 with the corresponding computation in [49] of which we are not able to pin down the origin: we stick to the result in [45]. Finally, in principle one should also consider the twist-2 gluon operator G a µρ G a ρν + (g µν /4)G a αβ G a αβ . This has been done e.g. in Ref. [50]. Due to the aforementioned cancellation between terms in the scattering amplitude, its contribution might be relevant while dominated by uncertainties. Our estimate shows it to be subdominant anyway with respect to the gluon operator that we include: following [51], we do not include it in this analysis.

Here q (N ) (2) and q(N) (2) are the second moment of the parton distribution functions for a quark q or antiquark q in the nucleon N , whose values are given e.g. in [49].

c N 1 ( , m X ) = 2 3 8 π g 4 M W m 2 N m X 1 m 2 h q=u,d,s f (N ) T q + 2 27 f (N ) T G (1 -2 3 Q κ Q ) - 1 M 2 W 4 81 f (N ) T G + 1 2 q=u,d,s,c,b (q (N ) (2) + q(N) (2)) , (79a)
c N 4 ( , m X ) = 2 g 4 π m N M W q ∆ (N ) q , ( 79b
)
where the notation makes apparent that the role of the generic variable λ (employed in Sec. 2) is here played by the parameter. At this point we can make explicit use of Eq. ( 22) and write an expression for λ 2 B involving the above coefficients and the Y

(N,N ) 1,1 (m X ) and Y (N,N ) 4,4
(m X ) rescaling functions provided in the website. We just plug this explicit expression for λ 2 B into the TS so that we have a TS( , m X ) and draw a contour plot in order to get the resulting bound on . In Fig. 11 we show the result. This figure illustrates, for instance, that fairly large values of are still allowed by the current experiments. 12

Conclusions

In this paper we have described a method and provided a self consistent set of tools for obtaining bounds from some of the current leading direct detection experiments on virtually any DM model. The method is based on the use of non-relativistic operators to describe the interaction between DM and ordinary matter. It builds on the results obtained by Ref. [3] and incorporates into them the necessary detector and astrophysical ingredients. Indeed, our main result consists in building up a set of integrated form factors F which allow to compute the detector response for any combination of NR operators. Our main output, in turn, consists in providing (on this website) a set of Test Statistic functions (one for each experiment we consider) and a set of functions Y, which are just a rescaling of the integrated form factors. Armed with these, one can obtain a constraint on any parameter of a DM model by using the recipe spelled out in Sec. 6. The TS functions that we provide for each experiment reflect the current status of results: they will change whenever new results are released and we will update them accordingly. On the other hand, the integrated form factor (and, a fortiori, the scaling functions) change only if major changes occur in the experimental set-ups, e.g. in efficiencies, cuts, and thresholds. We do not foresee frequent updates for those, but we will consider them when necessary. As the NR operators describe any possible DM interaction with ordinary matter, the method q : from DarkSUSY [38], EOS [54], CC [55] and micrOMEGAs3.1 [56] (see Table 4 and Appendix B). The green dot marks the location of the 'pure' MDM candidate with = 1 and mass fixed to m X = 9.6 TeV by the relic density requirement [48].

is fully model-independent and can be applied to any DM high energy construction. Indeed, somewhat as an aside, we also reviewed the procedure to pass from high-energy effective operators (into which a large class of models can be cast) to NR operators.

We now conclude with a brief discussion of some features of the method and some future perspectives. We stress that the method is exact for 'counting experiments', i.e. those which can be considered as having a single energy bin (like Xenon100 and Cdms-Ge), but the approximation for experiments with more bins is very well verified (see the discussion in Sec. 3). We remind that we considered only 'null result' experiments, and only four of them, but clearly the method can be extended to any experiment. Also, the form factors for other nuclei potentially of interest for DD should be computed. Moreover, different statistical tools from the ones we employed could be used. We encourage the experimental collaborations to release, when possible, their own likelihood and to present their computation of the integrated form factors, enriched by a knowledge of their detector deeper than that available to us.

Striving to devise model-independent methods for DM phenomenology (in direct detection but also more generally) is a crucial direction to pursue in order to make sense of the current exciting experimental panorama.

A Velocity integral

The DM velocity distribution f E ( v) in the Earth frame is related to the velocity distribution in the galactic frame f G ( v ) by the Galilean velocity transformation

f E ( v) = f G ( v + v E ),
where v E is the Earth velocity with respect to the galactic frame. One has that v E = v + v ⊕ , where v ⊕ is the Earth rotational velocity, and v the sum of the galactic rotational velocity of our local system and the Sun's proper motion. More details can be found for instance in Ref. [52,53].

In this paper we consider for f G ( v) a truncated Maxwell-Boltzmann, corresponding to an isothermal sphere density profile for the DM:

f G ( v) = exp(-v 2 /v 2 0 ) (v 0 √ π) 3 erf(v esc /v 0 ) -2v 3 0 π(v esc /v 0 ) exp(-v 2 esc /v 2 0 ) . ( 80
)
The denominator is responsible for the normalization of the velocity distribution,

v vesc d 3 v f G ( v) = 1 , (81)
where v esc denotes the escape velocity of DM particles from our galaxy. For definiteness, in this work we assume v esc = 544 km/s [8]. At zeroth order in the non-relativistic expansion, the scattering amplitude does not depend on the DM velocity, and therefore from Eq. ( 5) one has dσ T /dE R ∝ v -2 . This means that, when computing the rate in Eq. ( 6), one has to compute the integral

I 0 (v min ) = v v min (E R ) d 3 v f E ( v) v . (82)
When going at first order, the differential cross section might feature also a v-independent term, thus giving rise to

I 1 (v min ) = v v min (E R ) d 3 v v f E ( v) . (83)
For a truncated Maxwell-Boltzmann velocity distribution, these integrals can be computed analytically and result in [25] and

I 0 (v min ) = 1 2 v 0 η E [erf(η + ) -erf(η -)] - 1 √ π v 0 η E (η + -η -) e -η 2 esc (84) Ref. f (p) T u f (n) T u f (p) T d f (n) T d f (N ) T s ∆ (p) u ∆ (p) d ∆ (
I 1 (v min ) = v 0 η - 2 √ π η E + 1 √ π e -η 2 -- η + 2 √ π η E - 1 √ π e -η 2 + + v 0 4 η E 1 + 2η 2 E [erf(η + ) -erf(η -)] ( 85
) - v 0 √ π 2 + 1 3η E (η min + η esc -η -) 3 -(η min + η esc -η + ) 3 e -η 2 esc ,
where we defined the normalized velocities

η E ≡ v E /v 0 , (86)
η esc ≡ v esc /v 0 , (87)
η min (E R ) ≡ v min (E R )/v 0 , (88) and η
± (E R ) = min(η min (E R ) ± η E , η esc ) . (89)

B From quarks and gluons to nucleons

Given the very low energies involved in DM scattering with matter, the relevant degrees of freedom to be considered are not quarks and gluons, but rather nucleons and nuclei. In Sec. 2 we recalled the formalism introduced in Ref. [3], allowing one to write the DM-nucleus cross section from the interaction with the nucleons. However, many models are formulated in terms of interactions with the fundamental degrees of freedom, namely quarks and gluons.

In this Appendix we summarize the formulas connecting amplitudes at the quark and gluon level with the matrix element with the nucleons, and provide connections to the relevant literature. In particular, we will review the gluon matrix elements N |G a µν G a µν |N and N |G a µν Ga µν |N , and the matrix elements of quark bilinears N |qΓq|N with Γ either 1, γ 5 , γ µ , γ µ γ 5 , σ µν or σ µν γ 5 . Whereas these quantities can be extracted from experimental data, an effort has been recently undertaken to compute them within lattice simulations. However, still no general agreement has been reached in some cases, and quantities such as the strange quark contribution to the nucleon mass and spin are still plagued by sizeable uncertainties. For a comparison, we collect some of the numerical values reported by few different collaborations in Table 4. We consider four standard references, Ref. [38,54,57,58], and the two more recent Ref. [55,56] (notice that Ref. [56] is the 'update' of [58]).

B.1 Scalar couplings

The matrix element of the scalar quark bilinear N |qq|N , at zero momentum transfer, can be computed in the following way [59]. At zero momentum m N = N |Θ µµ |N , where

Θ µµ = q m q qq + β(α s ) 4α s G a µν G a µν (90)
is the trace of the QCD energy-momentum tensor after applying the equations of motion; the gluon contribution comes from the trace anomaly. Now one can integrate out the heavy quarks h = c, b, t via the heavy quarks expansion, yielding at the lowest order a result that is reproduced by the substitution

m h hh → - α s 12π G a µν G a µν . (91)
Expanding the beta function in powers of α s we get finally, at the lowest order,

Θ µµ = q=u,d,s m q qq - 9α s 8π G a µν G a µν . (92)
The gluon contribution can then be expressed in terms of light quarks via

- 1 m N 9α s 8π N |G a µν G a µν |N = 1 - q=u,d,s f (N ) T q ≡ f (N ) T G , (93)
where the quantities

f (N ) T q ≡ N |m q qq|N m N (94)
express the light quark contributions to the nucleon mass. Values of the f (N )

T q coefficients reported by a few different collaborations are collected in Table 4.

We are now interested in the coupling of a scalar operator S to the nucleon, its interaction with the quarks being dictated by the Lagrangian L Sq = S q c q qq (95) with c q the interaction coefficients. This interaction induces a coupling of S to gluons via a quark loop, so to be as generic as possible we can also add an independent S-gluons scalar interaction with the same form [59],

L Sg = c g Λ α s 12π S G a µν G a µν , (96)
where the numerical factors have been chosen for later convenience; this coupling is generated again at the loop level, Λ being connected to the mass of the particles running in the loop, and c g being some coefficient whose value reflects the underlying process. The S-nucleon interaction can be written now as

L SN = N |L Sq + L Sg |N N N ≡ c N S N N , (97)
with N the nucleon field, and where using Eq. ( 91) and (93) we can express c N as

c N = q=u,d,s c q m N m q f (N ) T q + 2 27 f (N ) T G q=c,b,t c q m N m q -c g m N Λ . (98)

B.2 Pseudoscalar couplings

The pseudoscalar coupling of an operator P with quarks is L P q = P q c q q iγ 5 q .

(99)

Integrating out the heavy quarks h = c, b, t, one has a loop-induced coupling of P with gluons as in the previous case, reproduced by the substitution valid at lowest order [59]

m h h iγ 5 h → - α s 16π G a µν Ga µν , (100)
where Ga µν ≡ ε µνρσ G a ρσ . We can also add, as above, an independent P -gluons coupling, generated for instance by other particles in the loop. This can be written effectively as

L P g = c g Λ α s 8π P G a µν Ga µν , (101)
where again Λ is connected to the mass of the particles running in the loop, and the numerical factors have been chosen for later convenience.

Evaluating the gluonic operator between nucleon states is a problematic task; we rely on the analysis performed in [55,60], based on the relation

N |ū iγ 5 u + d iγ 5 d + s iγ 5 s|N = 0 , (102)
derived from the large-N c and chiral limits. This leads to

N α s 8π G a µν Ga µν N = m N m q=u,d,s ∆ (N ) q m q , (103)
where m ≡ (1/m u + 1/m d + 1/m s ) -1 . The coefficients ∆ (N ) q , defined by 2∆

(N ) q s µ = N |qγ µ γ 5 q|N with s µ the nucleon spin four-vector, parametrize the quark spin content of the nucleon N . They involve therefore an axial-vector quark operator rather than a pseudoscalar one, but the two are related by PCAC [60]. These coefficients are argued to be negligible for heavy quarks [61], while for light quarks they satisfy the following relations: ∆

(p) u = ∆ (n) d , ∆ (p) d = ∆ (n) u , ∆ (p) s = ∆ (n)
s . As for the scalar charges f (N ) T q , we collect in Table 4 values of the coefficients ∆ (N ) q reported by few different collaborations. The P -nucleon interaction can be now written as

L P N = N |L P q + L P g |N N iγ 5 N ≡ c N P N iγ 5 N , (104)
yielding [55] c

N = q=u,d,s m N m q (c q -C) + c g m Λ ∆ (N ) q , (105)
where we defined C ≡ q c q m/m q .

B.3 Vector couplings

The coupling of a vector operator V µ with the quark current qγ µ q, L V q = c q V µ qγ µ q , (106) leads to the same vector interaction with the nucleon current,

L V N = c N V µ N γ µ N . (107)
Since the quark current is conserved, the nucleon charge is obtained by the quark charges by merely summing over the valence quarks of the nucleons. The interaction coefficients are therefore c p = 2c u + c d and c n = c u + 2c d , for proton and neutron respectively.

B.4 Axial-vector couplings

An axial-vector interaction of quarks with an operator A µ ,

L Aq = c q A µ qγ µ γ 5 q , (108)
leads to an axial-vector interaction with nucleons

L AN = c N A µ N γ µ γ 5 N , (109)
with c N = q c q ∆

(N )

q . The quantities ∆ (N ) q have been introduced in Sec. B.2, and values quoted by a few collaborations can be found in Table 4.

B.5 Tensor couplings

The quark tensor interaction with an operator T µν , L T q = c q T µν q σ µν q , (110) leads to a nucleon tensor interaction

L T N = c N T µν N σ µν N , (111)
with c N = q c q δ (N ) q [62]. Using the relation σ µν γ 5 = (i/2) µνρτ σ ρτ one can also find the relative expression for an axial-tensor coupling. The tensor charges δ (N ) q are interpreted as the difference between the spin of quarks and the spin of anti-quarks in nucleons. They have been recently measured by Ref. [63] as δ

(p) u = 0.54, δ (p) d = -0.23 at Q 2 = 0.8 GeV 2 ;
a newer global analysis with more data has been presented by the same collaboration in Ref. [64], with results compatible with the previous ones within the quoted uncertainties. Ref. [56,58] quote instead the values δ = -0.05. See also [65] for other determinations of these values. Supposedly, the proton tensor charges are related to the neutron ones in the same way as the axial charges ∆ (N ) q are.

1 Addendum: LUX 2013 data On October 30 th , 2013, the Lux collaboration announced their first DM search results in Ref. [1A]. With the data collected in about 85 live-days, they were able to set a stronger bound on the spin-independent interaction cross section, with respect to the preexisting limits. Given the relevance of this result, we include it in our set of numerical tools (Release 2) to derive bounds from direct DM searches.

Description of Lux

The Large Underground Xenon (Lux) experiment, operated at the Sanford Underground Research Facility in South Dakota, is a dual-phase xenon time-projection chamber. As the Xenon100 experiment, it uses the prompt scintillation (S 1 ) and ionization (S 2 ) signals to reconstruct the deposited energy and to discriminate nuclear recoils from electron recoils.

In Ref. [1A], a non-blind analysis was conducted on data collected with an exposure w of 85.3 live-days × 118.3 kg of fiducial volume. After cuts, 160 events were found within the S 1 energy region 2 -30 photoelectrons. The collaboration found that all the events are compatible with the expected electron recoil background distribution.

Due to lack of detailed information on expected background and event distribution in the S 1 -S 2 space, it is very difficult to devise a coherent and spectrum independent analysis of the data. In Ref. [2A], bounds were computed by assuming that either 0, 1, 3, 5 or all the 24 events below the electron recoil band in Fig. 4 of Ref. [1A] are indistinguishable from the expected signal, while the remaining are background events. The statistical analysis adopted in Ref. [2A] does not need information on the expected background, however our analysis (described in Sec. 3) requires that information and therefore we can not use the same approach. We proceed therefore as follows. Motivated by the fact that calibration neutrons are expected to mimic the S 1 -S 2 signal of heavy DM particles in the detector, we assume that DM events on the S 1 -log 10 (S 2 /S 1 ) plane distribute evenly above and below the mean of the nuclear recoil event distribution (solid red line in Fig. 3 and 4 of Ref. [1A]). We restrict ourselves to the region below the line, where N bkg = 0.64 electron recoil background events are expected while the neutron background is negligible. Of the 160 observed events only one is found below the nuclear recoil mean.

The DM recoil rate below the nuclear recoil mean in the Lux detector, assumed to be half of the total rate as discussed above, is computed as

R = 1 2 E max E min dE R (E R ) T dR T dE R , (112)
with E min = 3 keV nr and E max = 18 keV nr the average lower and upper nuclear recoil energy threshold as quoted on page 41 of Ref. [3A]. We obtain the efficiency (E R ) by interpolating the black crosses in Fig. 9 of Ref. [1A]. The integrated form factor is therefore given by F(N,N ) i,j

(m χ ) Lux = 1 2 w T ξ T E max E min dE R (E R ) F (N,N ) i,j (E R , T ) . ( 113
)
The rescaling functions defined in Eq. ( 23), Y 2 Addendum: SuperCDMS 2014 data On February 28 th , 2014, the SuperCdms collaboration reported the result of a first search for DM particles using the SuperCdms detectors [1B]. An earlier analysis [2B] had employed only one of the detectors, operated in a new mode (Cdmslite, for Cdms low ionization threshold experiment) that allowed to obtain an exquisitely low threshold thus setting a limit on DM-nucleon scattering cross sections for very light DM particles. In comparison, the present analysis, which employs data collected between October 2012 and June 2013, allows to set more stringent bounds in the high DM mass region. While the Lux bound on heavy DM particles [1A] is at present the strongest one for usual spin-independent scattering, the SuperCdms limit could become relevant when other interactions are taken into account, e.g. isospin-violating DM [7]. Therefore, we include the SuperCdms result in our set of numerical tools (Release 3) to derive bounds from direct DM searches.

Description of SuperCdms

SuperCdms is an upgrade of the Cryogenic Dark Matter Search (Cdms II) experiment (see Sec. 4.2), and is operated at the Soudan Underground Laboratory as its predecessor. The experiment consists of 15 germanium target crystals, each instrumented with ionization and phonon detectors. The measured ionization and phonon energies can be used to derive the recoil energy and the 'ionization yield', i.e. the ionization to recoil energy ratio which is used to distinguish signal from background.

A blind analysis with an exposure w = 577 kg • days, using only the seven detectors with the lowest trigger thresholds, revealed 11 candidate events in the range 1.6 -10 keV nr [1B]. The background prediction is N bkg = 6.1 +1.1 -0.8 (stat + syst) events, with a negligible additional 0.098 ± 0.015 (stat) events from radiogenic and cosmogenic neutrons. The collaboration has however reason to believe that the adopted background model does not correctly account for a feature in one of the detectors, thus misestimating the background. The decision was made, prior to unblinding the data, to report an upper limit on the DM-nucleon scattering cross section.

The DM recoil rate is

R = E max E min dE R (E R ) T dR T dE R , (114)
with E min = 1.6 keV nr and E max = 10 keV nr , and the efficiency (E R ) taken to be the red line in Fig. 1 of [1B]. The integrated form factor is given by

F(N,N ) i,j (m χ ) SuperCdms = w T ξ T E max E min dE R (E R ) F (N,N ) i,j (E R , T ) . ( 115
)
The rescaling functions defined in Eq. ( 23), Y (N,N ) i,j (m χ ) for contact operators and Y lr(N,N ) i,j (m χ ) for long range operators, are plotted for SuperCdms in Fig. 13.

Results and examples adding SuperCdms

We have modified Fig. 1, Fig. 7, Fig. 8 and Fig. 11 in the main text to add the SuperCdms bound. (m χ ) for the SuperCdms experiment.

Figure 1 :

 1 Figure 1: Our benchmark constraints, i.e. the constraints at 90% CL on the benchmark DM constant λ B entering Eq. (20) as a function of the DM mass from Xenon100 (solid black), Cdms-Ge (dashed blue), Coupp (dashed-dotted orange), Picasso (dotted red), Lux (long-dashed green) and SuperCdms(short-dashed dark cyan).

Figure 2 :

 2 Figure 2: Absolute value of the rescaling functions Y (N,N ) i,j

Figure 3 :

 3 Figure 3: Absolute value of the rescaling functions Y (N,N ) i,j

Figure 4 :

 4 Figure 4: Absolute value of the rescaling functions Y (N,N ) i,j

Figure 5 :

 5 Figure 5: Absolute value of the rescaling functions Y lr(N,N ) i,j

Figure 6 :

 6 Figure 6: Absolute value of the rescaling functions Y lr(N,N ) i,j (m χ ), in the case of interference among the NR operator O NR 4 and the NR long range operators O lr 5 and O lr 6 .

c p 5 , 6 = 2 c u 5 , 6 + c d 5 , 6 , c n 5 , 6 = c u 5 , 6 + 2 c d 5

 56565656565

  2b. Plug the above expression for λ B in the Mathematica interpolated function TS(λ B , m χ ) of the benchmark model, that we provide on the website. At this point you possess a TS(λ, m χ ) function and you can derive a bound on λ at the desired confidence level, e.g. draw a contour plot of TS = 2.71 for a 90% CL.

Figure 7 :

 7 Figure 7: Bounds on individual contact effective operators, all at 90% CL.

1 9 lead to O NR 1 and O NR 4

 914 for the SI case and O NR 4 for the SD. As specified in Eq. (49), these interactions correspond to the relativistic operators O N 1 = χχ N N and O N 8 = χγ µ γ 5 χ N γ µ γ 5 N for fermion DM in Eq. (44), respectively. While also O N 5 and O N

Figure 8 :

 8 Figure 8: The bounds on standard SI and standard SD cross sections obtained within our formalism.

Figure 10 :

 10 Figure 10: Two loop DM-gluon scattering for fermionic MDM with Y = 0. Figure adapted from [49].

  For the proton one has u(p) (2) = 0.22,d (p) (2) = 0.11, ū(p) (2) = 0.034, d(p) (2) = 0.036, s (p) (2) = s(p) (2) = 0.026, c (p) (2) = c(p) (2) = 0.019, b (p) (2) = b(p) (2) = 0.012 (for the neutron one needs to exchange the values for the up quark with those for the down quark). Making use of Eq. (49), one determines that the NR operators for the model are O NR 1 and O NR 4 , with coefficients

Figure 11 :

 11 Figure 11: Bounds on in MDM-inspired model. For each experiment we use four different sets of values of the scalar f (N ) T q , f (N ) T G and axial charges ∆ (N )

  Grant (agreement n. 278234 -'NewDark' project). M.C. is supported in part by the French national research agency Anr under contract Anr 2010 Blanc 041301 and by the EU ITN network Unilhc. E.D.N. is supported in part by DOE grant DE-FG02-13ER42022. M.C. acknowledges the hospitality of the Institut d'Astrophysique de Paris (Iap) and the Theory Division of CERN where a part of this work was done. E.D.N. wishes also to thank M.C. and P.P.. But M.C. wishes to thank E.D.N. and P.P.. In turn P.P., guess what, thanks E.D.N. and M.C..

(

  N,N ) i,j (m χ ) for contact operators and Y lr(N,N ) i,j (m χ ) for long range operators, are plotted for Lux in Fig.12

  .

Figure 12 :

 12 Figure 12: Absolute value of the rescaling functions Y (N,N ) i,j

Figure 13 :

 13 Figure 13: Absolute value of the rescaling functions Y (N,N ) i,j

Table 1 :

 1 .2014 List of the experiments that we consider, the section where they are discussed and the corresponding references.

	Cdms-Ge	4.2	[13]
	Coupp	4.3	[16]
	Picasso	4.4	[17]
	Lux	Add. 1	[1A]
	SuperCdms Add. 2	[1B]

Table 2 :

 2 Effective exposures, observed counts and expected background for each data set of the Coupp experiment

	For fluorine and carbon the efficiency

Table 3 :

 3 Technical data of the Picasso analysis. Left: Exposures and average rates for each of the detectors used in the study. Masses are normalized to the mass of fluorine, i.e. kg = kg(F). From tables 1 and 3 of Ref.[17]. Right: Low energy thresholds defining the bins used in the analysis and respective rate fluctuations. Adapted from figure 5 of Ref.[17].

	Module	Average rate	Exposure w j	Threshold	Rate fluctuations
	number j	R j (cpd/kg)	(kg • days)	energy E thr k	∆R(E thr k )
	71	327.6	16.09	1.7	-6.0 ± 7.1
	72	134.2	17.69	2.9	-0.3 ± 1.8
	131	31.5	10.89	4.1	1.6 ± 9.0
	134	209.6	15.94	5.8	-0.2 ± 9.2
	137	69.9	16.33	6.9	0.05 ± 1.3
	141	25.2	13.37	16.3	1.4 ± 1.7
	144	60.8	6.18	38.8	-0.2 ± 1.7
	145	31.5	7.83	54.8	1.3 ± 4.7
	147	20.6	6.55
	148	20.0	3.43

Table 4 :

 4 A compilation of values of the scalar and axial charges from the indicated references.

	p)
	s

Note that this quantity coincides with what is denoted as a Lagrangian L in[3,4], e.g. in Eq. (55) of[3].

As we shall see in Sec. 5, gluons behave differently and need separate treatment.

ξ T = 10 3 N A m T ζ T /kg Ā, where N A = 6.022×10 23 is Avogadro's number, ζ T are the numeric abundances and Ā ≡ T ζ T A T .

From now on we will consider λ as representing exclusively those parameters whose value is unknown, and that one wishes to constrain using direct detection searches.

12,12 11,11 8,8 5,5 3,3 1,1 p,n n,n p,p p,n n,n p,p p,n n,n p,p p,n n,n p,p p,n n,n p,p p,n n,n p,p

10,10 9,9 7,7 6,6 4,4 p,n n,n p,p p,n n,n p,p p,n n,n p,p p,n n,n p,p p,n n,n p,p

The unknown details of L eff (E R ) at low energies, while being crucial for an understanding of the detector performances at low DM masses m χ ∼ 10 GeV, do not affect the Xenon100 bounds on heavier DM[START_REF] Aprile | XENON100[END_REF].

11,12 8,9 4,6 4,5 1,3 n,p p,n n,n p,p n,p p,n n,n p,p n,p p,n n,n p,p n,p p,n n,n p,p n,p p,n n,n p,p

Recently the collaboration also presented results based on the Si analysis[14,15]. We will not consider these here.

Fluorine gives the main contribution to the spin-dependent scattering cross section, while being also important in low energy spin-independent scattering. The large uncertainties on the probability of nucleating a bubble for scattering of a F nucleus in CF 3 I are therefore responsible for the width of the exclusion bands in figure6and (for low DM masses) in figure7of[16]. At large energies (and therefore for higher DM masses), the bound on the spin-independent interaction is driven instead by iodine, whose behavior is known with better accuracy thus making the band shrink.

These results are based on the computations in Ref.[45] and[49] (for the one loop diagrams) and

Or that, in turn, the pure MDM case is still far from the experimental sensitivity. This conclusion is different from that in[45] because of the reduction in the cross section due to the cancellations discussed above.

Acknowledgments. E.D.N. heartily thanks Liam Fitzpatrick and Ami Katz for useful discussions and explanations about their works [3,4]. We also thank Vincenzo Cirigliano and Andrea De Simone for discussions, and Giorgio Busoni for checking our code and pointing out two problems with three plots in Fig. 7 (now fixed). This work is supported by the European Research Council (Erc) under the EU Seventh Framework Programme (FP7/2007-2013) / Erc Starting

Results and examples adding Lux

We have modified Fig. 1, Fig. 7, Fig. 8 and Fig. 11 in the main text to add the Lux bound.

Acknowledgments. E.D.N. thanks Emilija Pantic for useful discussions.