Global limits and interference patterns in dark matter direct detection

We compare the general effective theory of one-body dark matter nucleon interactions to current direct detection experiments in a global multidimensional statistical analysis. We derive exclusion limits on the 28 isoscalar and isovector coupling constants of the theory, and show that current data place interesting constraints on dark matter-nucleon interaction operators usually neglected in this context. We characterize the interference patterns that can arise in dark matter direct detection from pairs of dark matter-nucleon interaction operators, or from isoscalar and isovector components of the same operator. We find that commonly neglected destructive interference effects weaken standard direct detection exclusion limits by up to one order of magnitude in the coupling constants.

The currently favored strategy in the analysis of dark matter direct detection experiments relies on drastic simplifying assumptions regarding the underlying dark matter-nucleon interaction [10,11]. For instance, the dark matter-nucleon interaction is often assumed to be independent of the momentum transferred in the scattering, and of the dark matter-nucleus relative velocity. Moreover, the possible interference of different dark matter-nucleon interactions is commonly neglected in the data analysis. At the same time, distinct experiments are usually interpreted separately, even when they are compatible and, most importantly, complementary.
Though this approach is a well motivated first approximation, the substantial progress recently made in the field, and the efforts planned for the next decade, motivate the exploration of more sophisticated strategies.
In this work we explore two extensions of the standard approach to dark matter direct detection. Firstly, we interpret current dark matter direct detection experiments without assuming the knowledge of the dark matter-nucleon interaction a priori. In other words, we do not artificially restrict the analysis to constant spin-independent or spin-dependent interactions. In contrast, we interpret current observations within the general effective theory of one-body dark matter nucleon interactions, which predicts 28 independent isoscalar and isovector dark matter-nucleon interaction types. Secondly, we interpret current direct detection experiments globally [12,13], i.e. varying all correlated coupling constants simultaneously, and considering several direct detection experiments in a single multidimensional Likelihood analysis.
Our approach allows to characterize the interference patterns arising in dark matter direct detection from pairs of dark matter-nucleon interaction operators, or from isoscalar and isovector components of a given operator. We refer to these interference patterns as multi-interaction interference effects. A key result of this work is to quantify the impact of multi-interaction interference effects on the calculation of dark matter direct detection exclusion limits.
Effective theories for dark matter-nucleon interactions were already explored in the context of dark matter direct detection in , in the analysis of neutrino telescope observations in [34][35][36][37][38], and in the context of helioseismology in [39][40][41][42]. The prospects for dark matter direct detection in general effective theories were first investigated in [26][27][28], and later in [32]. Ref. [26] also quantifies the bias in the dark matter particle mass and crosssection reconstruction induced by fitting observations assuming constant dark matter-nucleon interactions, when dark matter interacts with nucleons differently. The importance of multiinteraction interference effects in dark matter direct detection was first noticed in [13,43,44], and recently also in [32].
The paper is organized as follows. In Sec. 2 we review the non-relativistic effective theory of dark matter scattering from nucleons and nuclei. Sec. 3 is devoted to the statistical framework adopted in the analysis of the dark matter direct detection experiments introduced in Sec. 4. We present our results in Sec. 5, and conclude in Sec. 6. Important equations are listed in Appendix A.

Dark matter-nucleus scattering in effective theories
In this section we review the non-relativistic effective theory of dark matter scattering from nucleons and nuclei.

Interaction Hamiltonian
Assuming one-body dark matter-nucleon interactions mediated by a heavy spin-1 or spin-0 particle, the most general Hamiltonian density for dark matter-nucleus scattering is given by [17] It is the sum of A terms, one for each nucleon in the target nucleus. The A terms in the sum are linear combinations of 14 Galilean invariant quantum mechanical interaction operatorŝ O

The 14 operatorsÔ
(i) k are constructed from the momentum transfer operatorq, the relative transverse velocity operatorv ⊥ , and the dark matter particle and nucleon spin operators,Ŝ χ andŜ N , respectively [17]. In Eq. (2.1) t 0 = 1 is the identity in isospin space, t 1 = τ 3 is the third Pauli matrix, and m N is the nucleon mass. The isoscalar and isovector coupling constants, respectively c 0 k and c 1 k , are related to the coupling constants for protons and neutrons by c p k = (c 0 k + c 1 k )/2, and c n k = (c 0 k − c 1 k )/2. They have dimension mass to the power of −2.
The Hamiltonian density in Eq. (2.1) admits the following coordinate space representation [17]Ĥ where σ i denotes the set of three Pauli matrices representing the spin operator of the ithnucleon in the target nucleus. r i is the ith-nucleon position in the nucleus center of mass frame, and In the last expressionv ⊥ T =v ⊥ −v ⊥ N , wherev ⊥ N is an operator acting on the ith-nucleon space coordinate. Explicit coordinate space representations for the operatorsq,v ⊥ T , andv ⊥ N can be found in [37]. Eq. (2.2) shows that dark matter couples to the constituent nucleons through the nuclear vector charge and spin current (first line), the nuclear convection current (second line), and the nuclear spin-velocity current (last line). In Eq. (2.2) we have omitted the nuclear axial charge, since it does not contribute to dark matter-nucleus scattering cross-sections when nuclear ground states are eigenstates of P and CP , as it is commonly assumed for target nuclei in dark matter direct detection. Table 1. Non-relativistic quantum mechanical operators defining the general effective theory of onebody dark matter-nucleon interactions. Because of the nucleon mass, m N , in the equations above all operators have the same mass dimension. For simplicity, here and in the next sections we omit the nucleon index (i) adopted in Sec. 2.

Transition amplitude
We use Eq. (2.2) to calculate |M N R | 2 spins , i.e. the square modulus of the amplitude for dark matter-nucleus scattering. Assuming that nuclear ground states are eigenstates of P and CP , and averaging (summing) over initial (final) spin configurations, one finds 1 (2.4) The index k extends over the nuclear response functions W τ τ k (y) defined below in Eq. (2.7). Assuming the harmonic oscillator basis for single-nucleon states, y = (bq/2) 2 with (2.5) The 8 dark matter response functions R τ τ k in Eq. (2.4) depend on matrix elements of the operators in (2.3), and are listed in Appendix A. They are functions of the dark matter particle spin, of q 2 /m 2 N , and of where µ T and v are the reduced dark matter-nucleus mass and the dark matter-nucleus relative velocity, respectively. The 8 nuclear response functions W τ τ k (y) in Eq. (2.4) are quadratic in nuclear matrix elements, and defined as follows where |J, T, M T represents a nuclear state of spin J, isospin T , and isospin magnetic quantum number M T . The reduction operation, i.e. || · ||, is done via the Wigner-Eckart theorem. In Eq. (2.7), A and B can each be one of the following nuclear response operators The vector spherical harmonics, Y M LL1 (Ω r i ), are defined in terms of Clebsch-Gordan coefficients and scalar spherical harmonics: where e λ is a spherical unit vector basis. In Eq. (2.4), we use the notation W τ τ A (y) ≡ W τ τ AB (y) for A = B.
The 6 nuclear response operators in (2.8) arise from a multipole expansion of the nuclear charge and currents in Eq. (2.2). The multipole index L must be less than 2J. In our analysis of current direct detection experiments, we adopt the nuclear response functions of Ref. [21]. Nuclear response functions for dark matter capture by the Sun can be found in [37] in analytic form.

Scattering rate
In a dark matter direct detection experiment, the expected differential rate of scattering events per unit time and per unit detector mass is given by where m χ is the dark matter particle mass, ξ T is the mass fraction of the nucleus T in the target material, and ρ χ is the dark matter density in the solar neighborhood. In Eq. (2.10) the integral denotes an average over the local dark matter velocity distribution, f , in the galactic rest frame boosted to the detector frame, v min (q) = q/2µ T is the minimum velocity required to transfer a momentum q from the target nucleus to the dark matter particle, and  v e (t) is the time-dependent Earth velocity in the galactic rest frame. Here we consider a Maxwell- ) truncated at the local escape velocity v esc = 554 km s −1 , and with v 0 = 220 km s −1 .

Statistical framework
We compare the effective theory of dark matter-nucleon isoscalar and isovector interactions reviewed in Sec. 2 to current dark matter direct detection experiments in a global multidimensional statistical analysis. If not otherwise specified, for a given dataset d we assume the Likelihood function is the number of predicted scattering events at a given point in parameter space, k is the number of observed recoils in a given experiment, and µ B is the corresponding number of expected background events. The arrays c and η represent the 28 coupling constants of the theory and the nuisance parameters in Tab. 2, respectively. When the error σ B on the number of expected background events µ B is known, we analytically The constant term in the last line of Eq. (3.2) is an arbitrary normalization, that we choose such that ln L eff = 0 for µ S + µ B = k. We construct approximate 2D frequentist confidence intervals for pairs of model parameters from an effective chi-square defined as ∆χ 2 eff ≡ −2 ln L prof /L max , where L max is the maximum Likelihood, and L prof is the 2D profile Likelihood: Wilks' theorem guarantees that under certain regularity conditions the distribution of ∆χ 2 eff converges to a chi-square distribution with 2 degrees of freedom [46].
Here we mainly restrict the analysis to a frequentist approach in order to avoid the volume effects found in [13] studying a smaller sample of dark matter-nucleon interaction operators. However, in computing exclusion limits that do not involve marginalization, we also adopt a Bayesian approach, since it requires a smaller number of Likelihood evaluations. Within the Bayesian approach, we express exclusion limits on the coupling constants in terms of x% credible regions. Credible regions are defined as portions of the parameter space containing x% of the total posterior probability, and such that the posterior probability density function P(Θ|d) ∝ L eff (d|Θ)π(Θ) at any point Θ inside the credible region is larger than at any external point.
The function π(Θ) is the prior probability density function. We assume a uniform prior probability density function for the logarithm of m χ and c τ i , with τ = 0, 1 and i = 1, 3, . . . , 15. This choice allows us to sample the Likelihood function covering several orders of magnitude in all directions in parameter space. For the nuisance parameters we assume Gaussian priors, as shown in Tab. 2.
We sample the multidimensional Likelihood function (3.2) using the Multinest program [47][48][49]. We run Multinest with parameters set to n live = 20000 and tol = 10 −4 . Finally, we use our own routines to evaluate scattering rates and the Likelihood function, Getplots [50] to calculate the profile Likelihood, and Matlab to produce the figures.

Datasets
Here we list the experiments considered in the analysis. For each experiment we provide information needed to evaluate the corresponding Likelihood function. We refer to [13] for further details. In the analysis, we convolve Eq. (2.10) with a Gaussian probability distribution function of standard deviation σ to account for finite energy resolution effects.

CDMS-Ge
We consider data from the final exposure of the CDMS II experiment, corresponding to an exposure of 612 kg-days [1]. In the data analysis the CDMS collaboration has observed k = 2 candidate events within the signal region 10 -100 keV, with an expected number of background events in the same energy interval equal to µ B = 0.9 ± 0.3. Here we assume a maximum experimental efficiency of 32% at 20 keV, linearly decreasing to 20% at 10 keV and at 100 keV. In addition, we consider a Gaussian energy resolution with standard deviation σ = 0.2.

CDMS low threshold
In the 2 keV threshold analysis performed by the CDMS collaboration, the T1Z5 germanium detector has observed k = 36 events in the 2 -20 keV region [2]. Zero-charge, surface, and bulk events, can explain 75% of the observed counts [2]. In computing the Likelihood function, we therefore consider µ B = 36 × 0.75, with σ B = 0. We assume the efficiency reported in the inset of Fig. 1 in [2], a Gaussian energy resolution with energy dependent standard deviation σ = 0.293 2 + (0.056E R ) 2 , and an exposure for T1Z5 of 241/8 kg-days.

SuperCDMS
The SuperCDMS experiment has collected data using 7 germanium detectors, with a total exposure of 577 kg-days [3]. Tab. 1 in [3] provides details regarding the number of recorded events, as well as the number of expected background events for each detector separately. In our analysis we consider data from 5 detectors only, hence neglecting T5Z2 and T5Z3, since for these detectors the number of observed counts is significantly larger than the number of expected background events. In the calculations we consider a Gaussian energy resolution with σ = 0.3, and the detector efficiency of Fig. 1 in Ref. [3]. We assume an average exposure per detector of 577/7 kg-days. In the Likelihood function we set µ B and σ B as in Tab. 1 of Ref. [3].

CDMSlite
The SuperCDMS experiment has also collected data in the low ionization threshold operating mode (i.e. CDMSlite) during 10 live days of dark matter search [4]. The nuclear recoil energy threshold for this experimental configuration is 170 eV ee , corresponding to 841 eV nr . The average rate of nuclear recoils in the CDMSlite detector is 5.2 ± 1 counts/keV ee /kg-day between 0.2 and 1 keV ee , and 2.9 ± 0.3 counts/keV ee /kg-day between 2 and 7 keV ee . For these data we assume the Likelihood function where R [0.2,1] and R [2,7] represent the average rates between 0.2 and 1 keV ee , and between 2 and 7 keV ee , respectively. We calculate the expected average count rate in the CDMSlite detector using the efficiency reported in the inset of Fig. 1 in [4], and assuming perfect energy resolution.

XENON100
We use the data presented by the XENON100 collaboration in [5], corresponding to an effective exposure of 34×224.6 kg-days. XENON100 has observed k = 2 candidate signal events in the pre-defined nuclear recoil energy range 6.6 -30.5 keV nr , corresponding to 3 -30 photoelectrons. The number of expected background events in the same interval is 1.0 ± 0.2. The differential spectrum of expected photoelectrons S1 in XENON100 is given by In Eq. (4.2), n is the number of photoelectrons actually produced from a recoil energy E R . For n we assume a Poisson distribution of mean ν(E R ) = E R L eff (E R )L y S nr /S ee , where L y = 2.28 ± 0.04 PE/keV ee is the light yield, S ee = 0.58 and S nr = 0.95 are the electric field scintillation quenching factors for electron and nuclear recoils, and L eff (E R ) is the energy dependent scintillation efficiency. We assume for L eff (E R ) the following form is the best-fit scintillation efficiency reported in Fig. 1 of Ref. [51], and ξ Xe is a nuisance parameter, first proposed in Ref. [45] to logarithmically extrapolate the scintillation efficiency below 3 keV. In Eq.

XENON10
In a second study the XENON collaboration uses the ionization signal S2 only to measure the nuclear recoil energy of the detected events [6]. The experimental recoil energy threshold for this analysis is of about 1.4 keV nr . The relation between S2 (measured in PE) and the observed nuclear recoil energy is S2 = Q y (E R )E R . We extract the function Q y (E R ) from Fig. 1 in [6]. The experimental exposure corresponding to these data is 12.5×1.2 kg-days. In this low energy threshold analysis XENON10 observed k = 23 candidate events in the signal region 1.4 -10 keV nr . A large background contamination cannot be excluded [6]. We calculate the expected nuclear recoil energy spectrum in XENON10 as in Eq. 5.5 of [10], assuming a constant efficiency E(E R ) = 0.94. We also allow for µ B = 20 background events of unknown origin, with σ B = 10.

LUX
We construct the LUX Likelihood function as for XENON100. The first data release of the LUX experiment consists of 85.3 live days of dark matter search data, corresponding to an exposure of 250×85.3 kg-days [7]. LUX observed 160 events with S1 between 2 and 30 PE. Only one event is below the mean of the Gaussian fit to the nuclear recoil calibration events in Fig. 4 of [7]. The expected number of background events in the same region is 0.64 ± 0.16. We assume σ PMT = 0.37, and the experimental efficiency in Fig. 9 of Ref. [7] divided by 2 to account for the 50% nuclear recoil acceptance quoted by the LUX collaboration.

COUPP
Using a 4.0 kg CF 3 I bubble chamber, the COUPP experiment has measured k = 2, k = 3 and k = 8 bubble nucleations above the threshold energies 7.8 keV nr , 11 keV nr and 15.5 keV nr , respectively [8]. We calculate the expected number of dark matter scattering events above a threshold energy E th as follows [8] where (E th ) is the threshold dependent experimental exposure times the bubble detection efficiency, and P T (E R , E th ) is the probability that an energy E R nucleates a bubble above E th . The latter is given by [8] Here we neglect scattering from Carbon (α C = 0), assume perfect efficiency for bubble nucleation for Iodine (α I → +∞), and treat α F ≡ a COUPP as a nuisance parameter. For the energy dependent exposure we assume (7.8 keV) = 55.8 kg-days, (11 keV) = 70 kg-days and (15.5 keV) = 311.7 kg-days [8]. For the three threshold configurations the expected number of background events is µ B = 0.8, µ B = 0.7 and µ B = 3, respectively.

PICASSO
The PICASSO experiment searches for dark matter using superheated liquid droplets made of C 4 F 10 [9]. In the last run PICASSO operated assuming eight different bubble nucleation threshold energies, namely 1.7, 2.9, 4.1, 5.8, 6.9, 16.3, 38.8 and 54.8 keV. For each energy, the observed count rate above thresholdR i , i = 1, . . . , 8, is reported in Fig. 5 of Ref. [9], including the associated errors σ i . We calculate the expected scattering rate R i in the energy range [E th , +∞] using Eq. (4.3) with = 1. As for COUPP, we neglect scattering from Carbon, and consider α F = a PICASSO as a nuisance parameter. For the eight data pointsR i we assume the Gaussian Likelihood function: (4.5)

Global limits and interference patterns
We now compare the general effective theory of one-body dark matter-nucleon interactions to current observations, in a global statistical analysis of the direct detection experiments in Sec. 4. We present our results in terms of exclusion limits on the 28 coupling constants of the theory. Our investigation focuses on multi-interaction interference effects in the exclusion limit calculation. Multi-interaction interference effects are of two types. A first type involves pairs of dark matter-nucleon interaction operators. It generates terms proportional to c τ i c τ j , with i = j, in Eq. (2.10), arising from R τ τ Σ , R τ τ Σ , R τ τ Φ M and R τ τ ∆Σ . Inspection of Eq. (A.1) in Appendix A shows that there are 7 pairs of interfering operators in the general effective theory of one-body dark matter-nucleon interactions, namely (Ô 1 ,Ô 3 ), (Ô 4 ,Ô 5 ), (Ô 4 ,Ô 6 ), (Ô 8 ,Ô 9 ), (Ô 11 ,Ô 12 ), (Ô 11 ,Ô 15 ), and (Ô 12 ,Ô 15 ). Other pairs of operators do not interfere, partially since nuclear ground states are eigenstates of P and CP , and partially because of theirŜ χ dependence. We refer to the remaining operators in Tab. 1 as non-interfering dark matter-nucleon interaction operators. c τ i c τ Table 3. A second type of multi-interaction interference arises when isoscalar and isovector components of the same operator combine in Eq. (2.10) forming terms proportional to c τ i c τ i , τ = τ . Any operator in Tab. 1 can produce isoscalar-isovector interference patterns of this type as long as off-diagonal nuclear response functions are generated in the dark matternucleus scattering.
Let us now consider a pair of operators (Ô i ,Ô j ) and an experiment characterized by the signal region S. Let us also assume c τ i = 0 and c τ j = 0 for fixed (τ, τ ) and (i, j) pairs, and all other coupling constants equal to zero. Then, contours of constant µ S (m χ , c, η) at a given m χ and η represent ellipses in the c τ i − c τ j plane of equation where the constants a τ τ ii , a τ τ jj , and a τ τ ij are given by and similarly for a τ τ jj and a τ τ ij . We can therefore introduce a correlation coefficient r τ τ ij to quantify the degree of interference between a pair of operators (Ô i ,Ô j ), or between the isoscalar and isovector components of a given operator: In Tab. 3 we list some of the r τ τ ij coefficients different from zero for the LUX experiment at m χ = 10 TeV.
Below we fit m χ , η, and c to observations globally, i.e. including all the available data, and simultaneously varying sets of coupling constants with r τ τ ij = 0.

Non-interfering operators
We start with an analysis of the non-interfering operators. Non-interfering interaction operators areÔ 7 =Ŝ N ·v ⊥ ,Ô 10 = iŜ N ·q/m N ,Ô 13 = i(Ŝ χ ·v ⊥ )(Ŝ N ·q/m N ) and O 14 = i(Ŝ χ ·q/m N )(Ŝ N ·v ⊥ ). In contrast toÔ 7 andÔ 10 , the interaction operatorsÔ 13  Figure 1. In all panels, cyan contours represent 2D 90% confidence intervals from a fit of m χ , η and of the coupling constants in the legends to the direct detection experiments indicated in parenthesis. In each panel, colored regions correspond to the 2D profile Likelihood associated with the cyan contours. Yellow contours denote 2D 90% credible regions obtained by fitting m χ and a single coupling constant to the LUX data. andÔ 14 depend on both the momentum transfer operator and the transverse relative velocity operator. All non-interfering operators contribute to the nuclear spin current through the nuclear response operators Σ LM ;τ and Σ LM ;τ , which also appear in the theory of electroweak scattering from nuclei, and characterize the familiar spin-dependent dark matter-nucleon interaction operatorÔ 4 . At the same time, the operatorÔ 13 induces a nuclear spin-velocity current through the nuclear response operatorΦ LM ;τ , which is specific to dark matter-nucleon interactions.
The top panels in Fig. 1 show the results that we obtain fitting c 0 7 , c 1 7 , η and m χ to current direct detection experiments. Results are presented in terms of 2D profile Likelihoods (colored regions) and associated 2D 90% confidence intervals (cyan contours) in the planes For reference, in the same panels we also report the 2D credible regions that we obtain fitting m χ , and a single coupling constant (c 0 7 in the left panel, and c 1 7 in the right panel) to the LUX data. For single coupling constant fits, where marginalization is not involved, we adopt a Bayesian approach, as it requires a relatively small number of Likelihood evaluation. Volume and prior effects in dark matter direct detection are discussed in [13,45,52,53]. The bottom panels in Fig. 1, and the four panels in Fig. 2 show the results of similar analyses which focus on the operatorsÔ 10 ,Ô 13 andÔ 14 .
Comparing different panels in Figs. 1 and 2, we find that confidence intervals and profile Likelihoods in the m χ − c 0 i and m χ − c 1 i planes are similar for a given operatorÔ i . This result is expected, in that the (0,0) and (1,1) components of the nuclear response functions W τ τ Σ , W τ τ

Isoscalar-isovector interference patterns
For any interaction operator in Tab. 1, including operators that do not interfere in pairs aŝ O 7 ,Ô 10 ,Ô 13 , andÔ 14 , we observe isoscalar-isovector interference patterns in the rate (2.10). As an example, let us focus on the operatorÔ 7 . The operatorÔ 7 contributes to the square  Figure 4. Same as for Fig. 3, but now for the operatorsÔ 11 ,Ô 12 , andÔ 15 .
modulus of the scattering amplitude as follows The last term in (5.4) describes a destructive interference between the isoscalar and isovector components of the operatorÔ 7 , as once integrated over a typical signal region, it gives a negative contribution to the total scattering rate. We observe similar interference patterns for all operators in Tab. 1.
In general destructive interference effects make direct detection exclusion limits weaker. For instance, an otherwise excluded value of c 0 7 can remain compatible with observations if compensated by an appropriately large value of c 1 7 , because of the negative ∝ c 0 7 c 1 7 term in Eq. (5.4).
For this reason, in all figures of this work single coupling constant fits to LUX data results in stronger exclusion limits for m χ 5 GeV, as in these fits destructive interference effects are neglected. For smaller values of m χ , SuperCDMS, and CDMSlite dominate the exclusion limit calculation.

Interfering operators
Interfering operators divide into 4 independent subsets. The first subset consists of the operatorsÔ 1 = 1 χN andÔ 3 = iŜ N · [(q/m N ) ×v ⊥ ]. The operatorÔ 1 contributes to the nuclear vector charge through the nuclear response operator M LM ;τ , whereas the operator O 3 contributes to the nuclear spin current, and to the nuclear spin-velocity current through Σ LM ;τ and Φ LM ;τ , respectively. The operatorsÔ 1 andÔ 3 generate a transition probability proportional to In Eq. (5.5) terms with τ = τ describe isoscalar-isovector interference effects similar to those discussed in Sec. 5.2. The term ∝ c τ 1 c τ 3 arises from the interference ofÔ 1 andÔ 3 . Integrating Eq. (5.5) over a typical signal region, cancellations between different terms occur, as the integrated nuclear response functions W τ τ M and W τ τ Σ , for τ = τ , and W τ τ Φ M , for τ = τ , are negative for many isotopes.
Because of cancellations in (5.5), numerical noise affects the exclusion limits derived simultaneously varying η, m χ , c 0 1 , c 1 1 , c 0 3 , and c 1 3 . To circumvent this problem, while exploring all interference patterns, here we present exclusion limits obtained in four complementary ways. In a first analysis we fit m χ , c 0 1 , c 0 3 , and η to current dark matter direct detection experiments. In a second analysis, we simultaneously place limits on the constants c 1 1 and c 1 3 while fitting the same data. Our third and fourth analysis respectively consider the constants c 0 1 − c 1 1 , and c 0 3 − c 1 3 as free parameters in the global fit. Fig. 3 shows the 2D 90% confidence intervals (colored contours) and profile Likelihoods (colored regions) resulting from the four analyses described above. To derive the contours in the panels we fit m χ , η, and the constants in the legends to observations. Results are presented in the m χ − c 0 1 , m χ − c 1 1 , m χ − c 0 3 and m χ − c 1 3 planes. For comparison, in each panel we also report the 2D 90% credible region that we obtain fitting m χ and a single coupling constant to the LUX data.
In the top panels of Fig 3 (corresponding to the m χ − c 0 1 and m χ − c 1 1 planes), we observe significant differences between distinct contours. As expected, limits obtained from single coupling constant fits to LUX data are generically stronger in the m χ 5 GeV mass region. Global fits of the c 0 1 and c 1 1 coupling constants to current observations (cyan contours) instead produce the weakest exclusion limits, as in Eq. (5.5) cancellations between the term proportional to c τ 1 c τ 1 , with τ = 0 or τ = 1, and the one proportional to c 0 1 c 1 1 can be large. In contrast, distinct contours in the bottom panels of to the nuclear vector charge through the nuclear response operator M LM ;τ , similarly toÔ 1 . The operatorsÔ 12 andÔ 15 induce a nuclear spin current through Σ LM ;τ and Σ LM ;τ , and a nuclear spin-velocity current through Φ LM ;τ .
The coupling constants c τ 11 , c τ 12 , and c τ 15 contribute to the scattering amplitude in a way similar to c τ 1 and c τ 3 in Eq. (5.5). Also in this case, we therefore expect cancellations in the scattering rate due to destructive interference between operator pairs, or between isoscalar and isovector components of the same operator. To study the impact of multi-interaction interference effects in the exclusion limit calculation, we proceed as for the operatorsÔ 1 and O 3 . We therefore divide the coupling constants c τ 11 , c τ 12 , and c τ 15 into 5 subsets: (c 0 11 , c 0 12 , c 0 15 ), (c 1 11 , c 1 12 , c 1 15 ) and (c 0 i , c 1 i ), with i = 11, 12, 15. We then fit m χ , η, and a single subset of coupling constants at the time to observations. Finally, from each fit we derive exclusion limits on the coupling constants in analysis. This procedure allows to separate interference effects due to operator pairs, from those related to isoscalar and isovector components of the same operator. Fig. 4 shows the 2D 90% confidence intervals and profile Likelihoods that we find in the five analyses described above. Results are presented in 6 m χ vs coupling constant planes. To derive the contours in the panels we fit m χ , η, and the constants in the legends to observations. For comparison, in each panel we also report the 2D 90% credible regions of a single coupling constant fit to the LUX data. In the m χ − c 1 12 , m χ − c 0 15 and m χ − c 1 15 planes distinct exclusion limits are comparable. In the other panels, contours differ by up to 1 order of magnitude in the coupling constants. As forÔ 1 andÔ 3 , such differences arise from multi-interaction interference effects. There are two additional subsets of interfering operators. One consists of the operatorŝ O 4 =Ŝ χ ·Ŝ N ,Ô 5 = iŜ χ · [(q/m N ) ×v ⊥ ], andÔ 6 = [Ŝ χ · (q/m N )][Ŝ N · (q/m N )]. The other one includesÔ 8 =Ŝ χ ·v ⊥ andÔ 9 = iŜ χ · [Ŝ N × (q/m N )]. Eq. (2.3) shows the contribution of these operators to the nuclear charges and currents. The operatorsÔ 5 andÔ 8 are the only operators contributing to the nuclear convection current (through ∆ LM ;τ ). Figs. 5 and 6 show the 2D 90% confidence intervals and profile Likelihoods that we find analyzing the operatorsÔ 4 ,Ô 5 ,Ô 6 ,Ô 8 , andÔ 9 . In these figures we do not observe strong destructive interference patterns.
We conclude this section with an analysis of selected pairs of coupling constants. We focus on largely correlated coupling constants (see Tab. 3). Figs. 7 and 8 show the 2D profile Likelihoods that we obtain fitting m χ , η and the parameters in the legends to current observations. Results are presented in the planes c 0 11 − c 0 12 , c 0 12 − c 0 15 , c 1 11 − c 1 12 , c 1 12 − c 1 15 , c 0 11 − c 1 11 , and c 0 12 − c 1 12 . In each panel, blue lines represent ellipses constructed from Eq. (5.1) with the right hand side fixed at a given reference value, and assuming LUX as experimental setup. As expected, the red elliptical regions in Figs. 7 and 8 are in general aligned with the blue ellipses. However, the 2D profile Likelihood in the c 0 12 − c 0 15 plane constitutes an exception. The mismatch observed in this figure is due to the fact that the blue ellipse in the top right panel of Fig. 7 neglects a large c 0 11 − c 0 12 correlation, which is instead taken into account deriving the corresponding red elliptical region.

Conclusions
We compared the general effective theory of one-body dark matter-nucleon interactions to current direct detection experiments in a global multidimensional statistical analysis. In this study, we included data from a variegated sample of direct detection experiments in a single Likelihood function. In the fits we simultaneously varied subsets of correlated model parameters, covering a large volume of the parameter space. We presented our results in terms of 90% confidence level exclusion limits on the 28 isoscalar and isovector coupling constants of the theory. We found that current direct detection experiments are able to probe all isoscalar and isovector coupling constants, including those measuring the strength of dark matter-nucleon interactions commonly neglected in this context. For instance, the interaction operatorŝ O 11 = iŜ χ · (q/m N ) andÔ 12 =Ŝ χ · (Ŝ N ×v ⊥ ), and the familiar spin-dependent interaction operatorÔ 4 =Ŝ χ ·Ŝ N are equally constrained by current data, though the former are by far less explored.
We characterized the interference patterns arising in dark matter direct detection from multi-interaction effects involving pairs of dark matter-nucleon interaction operators, or isoscalar and isovector components of the same operator. Most importantly, we quantified the impact of multi-interaction interference effects on the calculation of direct detection exclusion limits.
We found that destructive interference effects weaken direct detection exclusion limits derived neglecting the superposition of different interaction operators by up to one order of magnitude. Destructive interference effects largely affect exclusion limits on the strength of the familiar spin-independent interactionÔ 1 , and are less important for the standard spin-dependent operatorÔ 4 .

A Dark matter response functions
Below, we list the dark matter response functions that appear in Eq. (2.10). The notation is the same used in the body of the paper. In the figures, for definitiveness we assume j χ = 1/2.