A combined analysis of PandaX, LUX, and XENON1T experiments within the framework of dark matter effective theory

Weakly interacting massive particles are a widely well-probed dark matter candidate by the dark matter direct detection experiments. Theoretically, there are a large number of ultraviolet completed models that consist of a weakly interacting massive particle dark matter. The variety of models makes the comparison with the direct detection data complicated and often non-trivial. To overcome this, in the non-relativistic limit, the effective theory was developed in the literature which works very well to significantly reduce the complexity of dark matter-nucleon interactions and to better study the nuclear response functions. In the effective theory framework for a spin-1/2 dark matter, we combine three independent likelihood functions from the latest PandaX, LUX, and XENON1T data, and give a joint limit on each effective coupling. The astrophysical uncertainties of the dark matter distribution are also included in the likelihood. We further discuss the isospin violating cases of the interactions. Finally, for both dimension-five and dimension-six effective theories above the electroweak scale, we give updated limits of the new physics mass scales.


Introduction
The search for particle dark matter (DM) is one of the most important topics in both particle physics and astrophysics. Yet, no clear evidence appears in direct, indirect, or collider searches of the weakly interacting massive particles (WIMPs), which are the most intriguing DM candidates as motivated by the thermal production of DM and its relic density. Instead, several stringent limits have been reported, which pushes the WIMP mass heavier or the interaction between the DM and standard model (SM) particles weaker. Among those experiments, the direct detection ones have made rapidly improved limits on probing the interaction between DM and the quark sector of SM. Such a limit from xenon-type detectors has been improved by more than one order of magnitude in recent years, from XENON100 with 34 kg target [1], LUX with 250 kg [2], PandaX with 500 kg [3], and finally to a ton-scale detector XENON1T [4]. For a WIMP mass m χ ∼ 30 GeV, the latest -1 -

JHEP11(2017)024
XENON1T experiment sets upper limits on the spin-independent WIMP-nucleon cross section to ∼ 10 −46 cm 2 , which is only two orders of magnitude higher than the neutrino floor.
These stringent limits constrain DM models severely. For some well-motivated DM models, the parameter space has significantly shrunk or the survival region needs to be fine-tuned, such as the blind-spot region for the neutralino DM in the supersymmetric models, (see e.g. refs. [5][6][7] for the current status of the blind-spot region and refs. [8,9] for the latest comprehensive global study). On the other hand, the effective theory approach begins to catch more attention, such as the spin-1 mediator models including the Anapole, magnetic dipole, and electric dipole DM. The DM dipole interaction with the SM photon can be generated by a new mediator that is kinetically mixed with the SM photon. Such a mechanism can give rise to a velocity-dependent cross section [10]. In such a velocitydependent framework, the published experimental limits for spin-independent and spindependent WIMP-nucleon cross section cannot be applied directly. Although there are many more WIMP candidate ultraviolet (UV) completed models than the given examples here, some of them may result in similar phenomena in direct detection experiments at low energies. Therefore, a model-independent limit from these experiments will be very useful to link the WIMP models with the direct detection experiments.
Regardless of the model complexity, the DM direct detection is to search for recoil events due to the WIMP scattering off nuclei inside the detector. Because of the very small momentum transfer q ∼ O(MeV) (for WIMPs) compared with mediator masses, such an interaction can be expressed using the effective field theory (EFT) whose heavier mediators can be integrated out and only a new physical scale Λ and the information of spin and initial velocity are left. The low-energy EFT of DM has been extensively studied (e.g., see [11][12][13][14][15][16][17][18][19][20][21][22]). To study the EFT in a model independent framework, with possibly the spin and velocity dependence, a classification of 14 operators based on the spin and velocity of DM and target nuclei has been widely adopted in literature [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38]. 1 The interaction between the DM and nuclei is generally expressed by their response functions [23,25]. Any UV completed model can then be described with a combination of these 14 model-independent operators. This framework has also been adopted in studying the DM captured inside the Sun [40][41][42][43][44].
In this work, we update the current limits on the EFT operators of DM from the most recent PandaX, LUX, and XENON1T results, adopting similar methodology as that of ref. [25]. Comparing with previous EFT works, we have the following three main improvements. Firstly, we reconstruct three likelihoods of PandaX, LUX, and XENON1T with proper consideration of the astrophysical uncertainties. Secondly, using the reconstructed likelihood functions, we obtain combined limits on the DM mass and coupling parameters for the 14 operators as given in ref. [25]. For a more generic purpose, we also discuss the situation of the relativistic Lagrangian which can consist of more than one operator in the non-relativistic (NR) limit. We then discuss the isospin conserving (ISC) and isospin violating (ISV) scenarios. The ratios of maximum cancellation for ISV at the event rate level are computed. Finally, we show the updated lower bounds of new physics energy scale (or mediator mass) Λ for dimension-five and dimension-six DM effective theories above JHEP11(2017)024 Table 1. Non-relativistic quantum mechanical operators defining the general effective theory of one-body DM-nucleon interactions. Taken from ref. [25].
electroweak symmetry breaking scale. The code and the likelihoods of experiments will be incorporated in the LikeDM tool [45]. An alternative package GAMBIT [46] also solves the similar problems. This paper is organized as follows. In section 2 we briefly introduce the theory of DM direct detection in terms of effective operators. In section 3, we describe the experimental results of PandaX, LUX, and XENON1T, and their likelihood functions, including a detailed discussion of our treatment of the DM astrophysical parameters. In section 4 we present the results of our scans. We summarize our findings in section 5.

DM direct detection theory
When the Earth sweeps through the local halo together with the sun, the DM from local halo may hit the underground target nuclei via the DM-nucleus scattering. The DM direct detection is designed to detect the nuclear recoil energy due to such an interaction. Unfortunately, no firm detection of such interactions has been reported in present leading xenon-target experiments, such as PandaX-II [3], LUX [2], and XENON1T [4]. Usually the null result is interpreted as an upper limit on the cross section of the DM-nucleon interaction, and presented separately on the spin-independent and spin-dependent components. To date, the most stringent limit on the spin-independent cross section σ SI p comes from the XENON1T experiment [4].
Generally the exclusion limits by the experiment groups are based on simplified assumptions and cannot be used directly for a number of interactions. For example, a possible enhancement at small momentum transfer for the velocity dependent scattering cross section is usually not properly investigated in the standard spin-independent or spin-dependent approach by the experimental groups. To study the constraints on a variety of different operators, we adopted the particle model independent method developed in refs. [24,25].  operators in table 1 are the mass dimension six operators, it is convenient  to introduce some mass-scale parameters into the coefficients to make them dimensionless. Hence, by including the physical scale Λ and expanding the isospin index of coefficients, one can rewrite the operator as where τ 3 is the 3rd Pauli matrix and i runs from 1 to 15 for different operators with O 2 omitted. Following the convention in the code DMFormFactor [25], we set Λ to be the Higgs vacuum expectation value Λ = 246 GeV. The new coefficients c p i and c n i are treated as input parameters, and the relations c p i = c n i and c p i = c n i correspond to ISC and ISV scenarios, respectively.
In the limits of small recoil energy and slow DM velocity, one can expand all possible effective interactions with the following four three-vectors in the non-relativistic limit, where S χ and S N are spins of the DM and nucleon, the vector q is the transfer momentum, and v ⊥ is the velocity that is perpendicular to the momentum transfer q. The velocity v ⊥ is defined as v ⊥ = v + q/(2µ χN ) where µ χN is the dark matter-nucleon reduced mass. In classical limits, these four three-vectors describe the states before and after a scattering. Generally speaking, if the operators are nucleon spin-dependent, the vector S N appears in the operators whose cross sections are not coherently enhanced by a factor of target atom number square as the nucleon spin-independent case. Regarding the velocity dependent operators, they may come from the anapole or dipole interaction. For example, in ref. [38], it is shown that the anapole interactionχγ 5 χ depends on q 2 ; the magnetic dipole interactionχσ µν χ depends on q 4 and (qv ⊥ ) 2 ((v ⊥ /q) 2 ) for the case of mediator mass heavier (lighter) than the momentum transfer; the electric dipole interactionχσ µν γ 5 χ depends on q 2 (q −2 ) for the case of mediator mass heavier (lighter) than the momentum transfer. Following the conventions of ref. [32], we write the differential event rate of scattering between DM and the target nuclues per unit detector mass as a function of the recoil energy Q as dR dQ = where m χ and m T are the DM and target masses. The parameter ξ T is defined as where η T can be found in the website. 2 The differential cross section is  and the averaged amplitude can be written as The amplitude includes the nuclear response functions W τ τ k (y) and the DM response functions R τ τ k , which were described in ref. [25]. The indices τ and τ run over proton p and neutron n. We also briefly describe them in appendix A. The symbols M , ∆, Σ , Σ ,Φ , and Φ indicate the type of response function. Note that the dimension of the averaged amplitude is mass −4 . One shall not confuse the transfer momentum q in eq. (2.6) with the recoil energy Q in eq. (2.5). These two quantities are related by q 2 = 2m T Q.
In the astrophysics part, the DM velocity distribution function f ( v + v e ) can be described by the Maxwell-Boltzmann distribution [47], or directly extracted from N -body simulation [48]. The Earth's velocity in the galactic rest frame, v e , is added in order to translate the reference frame from the galactic rest frame to the Earth rest frame. In this work, we adopt the soft-truncated Maxwell-Boltzmann distribution [49]. The local density ρ 0 indicates the DM mass density near the Sun. We analytically integrate the velocity distribution over the solid angle and present the result in appendix B. The minimum DM velocity can be written as a function of recoil energy, assuming that the scattering between DM and nucleus is elastic where M r is the reduced mass. Finally, the efficiency of the dark matter detector (Q) needs to be included in the calculation. The total event rate is then We will see in the following section that the event rate R is the best model independent quantity in our likelihood functions.

PandaX-II experiment
The PandaX-II is a half-ton xenon dual-phase detector at the China Jinping underground Laboratory. The PandaX collaboration has published results based on their Run 8 (19.1 live days) and Run 9 (79.6 live days) data with exposure of 5845 kg-day and 27155 kg-day, respectively. For the Run 8, the total observed event number after all the cuts is 2 and the expected background event number is 2.4 ± 0.8. For the Run 9, the observed event number is 1 with an expected background event number of 2.4 ± 0.7. The 90% confidence level -5 -JHEP11(2017)024 (CL) upper limit on the (m χ , σ SI p ) panel has been shown in figure 5 of ref. [3]. Usually, the 90% CL limit in the experimental results of the DM direct detection, is computed with the two-tail convention, which corresponds to the 95% CL limit in the one-tail convention. We will adopt the one-tail convention in the analysis.
Based on the data of PandaX-II, one can simply build the likelihood function as [50] L PandaX ∝ i=run8,run9 Due to the complexity of the experimental performance one usually cannot exactly reproduce the σ SI p,95% line given in ref. [3]. Therefore, we introduce an additional correction factor f (m χ ) to account for this discrepancy, where s exp 95 (m χ ) is the event number inferred from the experimental limits, the σ SI p,95% line from ref. [3], and s stat 95 is the 95% CL limit computed via our likelihood analysis. The default values for astrophysical parameters, quoted by ref. [3], are v 0 = 220 km/s, ρ 0 = 0.3 GeV/cm 3 , and v esc = 544 km/s. s stat 95 is computed by eq. (3.1), which is DM model independent. Taking into account the correction factor f (m χ ), we have which is then used to compute the likelihood function given by eq. (3.1), instead of s i = RE i as done previously. With this correction factor included in the likelihood function, we are able to reproduce the PandaX results. Thus, in our analysis, the total likelihood of PandaX-II, L PandaX , is a function of both R and m χ .

LUX experiment
The LUX experiment is also a dual-phase xenon (250 kg) time projection chamber. Combining two data sets WS2013 and WS2014-16, the LUX has a similar exposure (3.35 × 10 4 kg-day) to that of the PandaX, and gives comparable but slightly stronger limits (perhaps due to a different analysis method 3 ), on the DM-nucleon cross section. The LUX JHEP11(2017)024 experiment did not report the total observed and expected events after cuts, which makes the construction of the likelihood function much more difficult. We construct a different likelihood function here than that in the PandaX case. By assuming that null signal has been detected, we build the total likelihood function of LUX as where s i,95 (m χ ) is the number of events computed from the 95% CL limit curve in ref. [2].
Here, the 95% CL is equivalent to 1.64σ far from the central value in 1-dimensional Gaussian likelihood. One has to bear in mind that the theoretical prediction of the event number s i depends on the specific particle model of DM, because the efficiencies of WS2013 and WS2014-16 are different and it is hard to compute the likelihood in advance with a particle model-independent quantity. Hence, for the purpose of a particle model-independent study, we employ two independent R variables and two likelihood functions. Since the result of σ SI p,95% of ref. [2] has been used in the likelihood eq. (3.4), we do not need to consider the correction factor f (m χ ) as in the PandaX likelihood.

XENON1T experiment
The XENON1T [4] is the first ton-scale xenon-type detector. With a fiducial mass of ∼ 1042 kg and a running of 34.2 live days, the XENON1T data give the currently most stringent limit on the spin-independent cross section between DM and nucleon. Like PandaX and LUX, the XENON1T is also a dual-phase detector, and all the setups of the PandaX likelihood can be directly applied for XENON1T, with only the replacement of o = 1 and b = 0.36 +0. 11 −0.07 [4]. A piecewise function is used to incorporate both the positive and the negative background uncertainties. The efficiency can be found in figure 1 of ref. [4]. Similarly, a correction factor f (m χ ), which is different from the one of PandaX, is needed to compensate the discrepancy between the official Xenon1T limit and our analysis.
One has to bear in mind that both XENON1T and PandaX groups used the unbinned likelihood analysis. On the contrary, our likelihood is a binned method which relies on the event number. In the binned method, all the information of each event are folded which makes our constraints weaker. Hence, our likelihood is more conservative than the actual ones by XENON1T and PandaX data.

Combined likelihood
Before combining three experimental data sets, it is necessary to consider the shared systematical uncertainties which may be multiply counted. Since LUX, PandaX, and XENON1T are three independent experiments, the systematical uncertainties from instrumentation are expected to be independent. However, the astrophysical uncertainties of the DM distribution should be the same. They should only be considered once in the likelihood calculation.  The combined likelihood can be simply written as the product of likelihoods of the three experiments, PandaX-II run8+run9, 4 LUX WS2013+WS2014-16, and XENON1T. In the left panel of figure 1, we show the upper limits of R at the 95% (solid) and 99% (dashed) CL for PandaX, LUX, and XENON1T. In all the mass region, the XENON1T limits are lower than the other two experiments. The capabilities of PandaX and LUX in constraining the DM-nucleon cross section are similar with each other, with slightly different sensitive mass regions. By comparing between the 95% CL (solid) and 99% CL (dashed) lines, we find that the PandaX has larger gap between the two lines than LUX, which can explain why LUX can constrain the parameter space more stringent than PandaX, although their exposure is comparable. In the right panel of figure 1, we give the 95% CL upper limits of the spin-independent cross section based on individual likelihood and the combined one. Our combined constraint improves by a factor of ∼ 1.3 compared with the XENON1T result.

DM astrophysical nuisance parameters
As aforementioned, the astrophysical uncertainties have to be properly considered, particularly for the operators where the cross section is velocity dependent. The DM velocity can be simply described by a soft truncated Maxwell-Boltzmann distribution (see eq. (B.4)). The parameters of the velocity distribution, the DM local velocity, the escape velocity, and the local density can be determined by the kinematics of stars or gas in the Milky Way [52]. Such methods are, however, subject to systematic uncertainties of the DM halo JHEP11(2017)024 Table 2. The DM astrophysical nuisance parameters given in ref. [52] which considers the rotation curve of the Milky Way within ∼ 100 kpc.
profile. See ref. [53] for a compilation of the results on the local density measurements by different methods and/or data sets. In this work we adopt the velocity parameters derived in ref. [52], which employed an updated precise measurement of the rotation curve using the LAMOST data to determine the DM halo properties. In table 2, we list the parameter values and their 1σ uncertainties of the local DM density, the local DM velocity, and the escape velocity given in ref. [52]. Unlike GAMBIT DDCalc [54] whose nuisance distribution of DM local density is log-normal, we adopt normal distribution. The uncertainties of the local velocity and escape velocity will be included in the likelihood calculation, via a profile likelihood approach.
Note that the integration of the DM velocity distribution over the solid angle is computationally heavy; we thus provide three analytical formulae valid in different velocity regions in appendix B. The local DM density affects the normalization of the event rate. We also quote the results given in ref. [53], i.e., 0.20-0.56 GeV cm −3 , which is indicated by a band in the plots of our final results.

Numerical results
In this section, the combined limits from the PandaX, LUX and XENON1T are presented on (i) the non-relativistic operator coefficients, (ii) the relativistic effective operator couplings, and (iii) the high energy EFT energy scale Λ. We only focus on the DM masses between 5 GeV and 1 TeV. Note that our likelihood functions are particle model independent but only the spin-1/2 DM is considered in this work as illustration. For DM with different spins, the limits can be different.
We present the limits at the 95% CL with the CL b hypothesis, namely −2∆ ln L = −2(ln L − ln L 0 ) = 2.71. Here L 0 is the likelihood without DM signal (background only hypothesis), and L is the corresponding likelihood for given DM model parameters. L 0 is slightly smaller than the global maximal likelihood in the background-and-DM hypothesis, due to a less-than-1σ excess for the XENON1T data. The choice of the backgroundonly hypothesis instead of the global maximum L as reference makes our results more conservative.
The astrophysical parameters, including the local velocity and escape velocity, are treated as nuisance parameters following Gaussian distributions with parameters given in table 2. As we have discussed above, a larger uncertainty band of the local density measurements [53] is adopted to account for possible systematic uncertainties.

JHEP11(2017)024
level. Unlike such traditional definitions, we define our ISV coupling ratio for maximum cancellation at the event rate level, in which we take the experimental efficiencies and the DM velocity distribution into account. This new definition shall be more generic. As a comparison, the ISC scenario (c n = c p ) is also presented.

Non-relativistic operators
In figure 2, the 95% upper limits for the operators that consist of only two basis threevectors given in eq.  . This is because that the contribution of the nuclear response function Φ is larger than M , since Φ possess not only the scalar contribution but also the quasicoherent one [24,25]. For a similar reason, the limits for O 12 is two orders of magnitude stronger than that of O 9 , since the operator O 12 also has the quasicoherent contribution from Φ .
The  [25], the structure S N · v ⊥ , where the DM couples to nucleon via the axial charge, can have a vanishing intrinsic velocity contribution so that it leads to standard spin-dependent operators. In addition, for the 5-vector combination O 15 , we found that the difference between ISC and ISV are the smallest among the 14 operators. This is also expected because the interference terms in O 15 are from high dimensional operators.
Finally, we want to discuss the ISV coupling ratio for maximum cancellation. In order to study the interference between the χ-n and χ-p amplitudes, we introduce an interference parameter N np and rewrite the total event number in a matrix form, as  where N nn is the predicted event number with c n i = 1 and c p i = 0, and N pp is the one with c n i = 0 and c p i = 1. The interference N np can be easily obtained from the numerical computation. Applying the minimum condition for the variable c n i /c p i , the minimum value of N (maximum cancellation) is located at Such a coupling ratio is determined at the event level, which can slightly vary for different targets, cut efficiencies, and DM velocity distributions. In figure 5, we present the predicted ratio of the coefficient c n 13 to c p 13 by the minimum condition (blue dashed line), compared with the maximum cancellation at the 95% CL (red solid line) from figure 4. The difference between these two is very small, suggesting that the value −N np /N nn gives an excellent estimate of the ISV coupling ratio. One may note that the ratio becomes a constant of 0.225 for O 13 if the DM mass is much heavier than 100 GeV.    Figure 6 shows the ratios c n i /c p i for the other 13 operators except O 13 . We can find that such ratios for all the operators approach constant at high m χ regions. The ISV coupling ratio c n i /c p i is between −0.82 to −0.62 for the operators shown in the left panel of figure 6, and between −0.043 to −0.015 for the operators shown in the right panel. From eq. (4.2), we can see that the sign of c n i /c p i for the maximum cancellation is always determined by the interference term N np . Interestingly, except for the operator O 13 where the interference is negative, the interference of all the other operators is positive so that c n i and c p i have opposite signs for maximum ISV cancellations.

Relativistic effective Lagrangians
In previous sections, we have mentioned that all the four-point effective Lagrangians can be expanded by the NR effective operators shown in table 1. In fact, many relativistic Lagrangians correspond to one operator in the non-relativistic limit with a simple relation- ship (see table 1 of ref. [25]). For example, the ratio of d 3 (the coupling of the interaction iχγ 5 χN N ) to c 11 (the coefficient of the operator iS χ · q m N ) is −m χ /m N . Again here we use the notation following ref. [25]. However, for some relativistic Lagrangian L i the expansions of operators are not so trivial and may contain several operators, such as Here we omitted the energy scale Λ in front of the interactions. It will be included in the computation as we did for the case of non-relativistic operators. The effective Lagrangians L 9 and L 10 represent the magnetic dipole interaction with proton and neutron, respectively. The electric dipole moment interaction with proton is part of L 18 , and the Anapole interaction can be presented by L 13 . For more explicit  expressions of the Anapole and electromagnetic moments in the effective theory, one can refer to ref. [67].
To reuse the code developed for effective operators, we can use the relationship presented in eq. (4.3). Taking L 13 as an example, if one sets its coupling as d 13 , the coefficients of effective operators for the event rate computation are c 8 = 2d 13 , c 9 = 2d 13 , and the other coefficients are zero.
In the left panel of figure 7, we present the upper limits on coupling d p 12 for the ISC (red solid) and ISV (blue dashed) scenarios. For the ISV scenario for all Lagrangians in eq. (4.3), only the Lagrangian L 12 has a positive maximum cancellation ratio. Therefore -15 -

JHEP11(2017)024
we show its results separately from the others. The maximum cancellation ratio of neutron to proton couplings for L 12 as a function of m χ is shown in the right panel.
The limits on the couplings for other Lagrangians are shown in figure 8, and the maximum cancellation ratios of neutron to proton couplings are summarized in figure 9.

High energy scale theory
The new physics models beyond the standard model often appear at higher energy scale (i.e. greater than Z boson mass), while the non-relativistic operators apply at the low energy scale. The match between these two scales may not be trivial. For example, the longdistance corrections due to DM scattering with a pion exchanged between two nucleons can generate a coupling c i p,n proportional to q −2 , which would significantly change the results [36,67]. In this subsection, we adopt the Mathematica package DirectDM [36, 37] 5 to calculate the relationship between the non-relativistic operators and the dimension-five and dimension-six effective DM-quark interactions. All contributions of quarks whose mass is less than Z boson are included.
Using similar conventions with that of ref. [36], we consider two dimension-five effective interactions, where F µν is the electromagnetic field strength tensor, and Λ is the new physics energy scale. The CP even operator Q 1 and the CP odd operator Q 2 represent the magnetic dipole and the electric dipole interactions, respectively. Because protons and neutrons couple to the electromagnetic field differently, and there is no electric dipole interaction between DM and neutrons, the isospin is not conserved in these dimension-five interactions. The interaction Q and DM-neutron coefficients c n 4 , c n 6 , in which three coefficients c p 5 , c p 6 , c n 6 are proportional to the inverse square of transfer momentum, q −2 . However, the interaction Q 2 can only be expanded by operator O 11 with the DM-proton coefficient c p 11 which is also proportional to q −2 . See appendix A of ref. [36] for the exact relations.
In addition, we consider four dimension-six effective interactions between DM and quarks (q) For the sake of simplicity, we only consider the scenario where all the DM-quark couplings are unified. Therefore, isospin is conserved under such a simple assumption. Compared with dimension-five interactions, the mappings from these four interactions to the low energy DM-nucleon coefficients do not depend on the transfer momentum. The interaction Q 1,q only relates to c p,n 1 . In addition to c p,n 1 , Q 2,q contains two more 5 We thank F. Bishara, J. Brod, B. Grinstein and J. Zupan for providing us the DirectDM code [36,37].  coefficients c p,n 8 and c p,n 9 . The interaction Q 3,q has one more relevant coefficient c p,n 7 added to that of Q (6) 2,q , and the interaction Q (6) 4,q has one additional relevant coefficient c p,n 4 added to that of Q (6) 3,q . In figure 10, we present the 95% lower limits on Λ for dimension-five interactions (left) and dimension-six interactions (right). It is shown that the constaints on the electric dipole interaction Q (5) 2 are stronger than that of the magnetic dipole interaction Q 4,q has a significant positive contribution from O 4 which makes the energy scale of Q (6) 4,q larger than that of Q (6) 3,q .

Summary
As the first time to consider a combined analysis of three most recent and powerful experimental data sets from PandaX, LUX, and XENON1T, we have approximately reconstructed their likelihood in terms of event rates and reported a new combined limit based on a spin-1/2 DM. To consider the possible impact of astrophysical uncertainties from the DM local velocity, escape velocity, and local density, we introduce these parameters as nuisance parameters in our likelihood.
In the low energy effective operator framework, we apply our combined likelihood for each effective operator and derive their 95% CL upper limits on the plane (m χ , c p i ) for the ISC and ISV (with maximum cancellations) scenarios. As expected, our combined limits of the effective operators are more stringent than previous studies [32,67]. In addition to the low energy effective operator framework, we also study the high energy effective Lagrangians which can usually be expressed as the combination of several operators. The 95% CL limits (upper or lower) on the (m χ , d p i ) and (m χ , Λ) planes are presented with several representative high energy effective Lagrangians.
The inclusion of the uncertainties of the DM velocity distribution and the escape velocity of the solar location leads to 5% uncertainties for the velocity independent operators and 7.5% for the velocity dependent operators. The maximum uncertainties appear around m χ 10 GeV whose recoil energies fall into the small efficiency region. However, such astrophysical parameters depend on the data sets and modeling of the luminous matter distribution, as well as the prior assumption of the DM density profile, which may be subject to additional systematic uncertainties. Therefore we quote a relatively large band of the local DM density measured by different analyses [53] to show a potential uncertainty range of our results.
For the xenon target detectors, we report new ISV coupling ratios c n i /c p i for the maximum cancellation between the contributions of DM-proton and DM-neutron couplings. Firstly, we find that the ratio is not a constant with respect to the DM mass. However, as long as the DM mass is heavy, the ratio asymptotically approaches a constant. Secondly, only the operator O 13 has a positive ratio, and the rest operators have negative ratios for the maximum cancellation. Finally, the well-known number of c n i /c p i = −0.7 only agrees with the operators O i=1, 3,5,8,11,12,15 . For the high energy effective Lagrangian cases, only the transfer momentum dependent operators can have a large ratio change with respect to the DM mass. The ratios are between their compositions of lower energy operators. In this paper, we only studied one-nucleon contributions. If the two-nucleon correlations would be included, the result of ISV case can be changed [61,68].
Lastly, we would like to suggest the future experiments to publish their limits of event rates R togehter with the spin-independent and spin-dependent cross sections. There are two advantages to present the limits on the plane (m χ , R). Firstly, this can help theorists to go beyond the neutralino-like benchmark scenario which is only relevant to -19 -

JHEP11(2017)024
operators O 1 and O 4 . As performed by the SuperCDMS collaboration [30] and XENON collaboration [31], the EFT operators can provide a more general framework to explore the DM-nuclear interaction. With the likelihood information of R, one can simply obtain the experimental limits for their interested operators. Secondly, it also helps to unfold the astrophysical uncertainties. The systematic uncertainties for the next generation DM direct detection detectors will be more important because the statistical precision can be much improved for the future ton-scale detectors and more kinematic information can be obtained for example in the DM directional detection experiments.

A DM response function
Following the definition of ref. [25], we use the notations M , ∆, Σ , Σ ,Φ , and Φ to represent the DM currents by the vector charge, vector transverse magnetic, axial transverse electric, axial longitudinal, vector transverse electric, and vector longitudinal operators, respectively. In the DM response functions, one needs also to consider two interference terms, Φ M and ∆Σ . There are eight possible DM response functions: B Integration of the DM velocity distribution over solid angle where We adopt the soft truncated Maxwell-Boltzmann distribution of f : where Θ(x) is the step function, and "Norm" is the normalization factor In the above formulae, v 0 , v e and v esc are all in the Galactic center (GC) frame. The relationship between the DM velocity v + v e in the GC frame and the velocity magnitude v in the Earth frame is v + v e = (v sin θ cos ϕ, v sin θ sin ϕ, v cos θ + v e ).
Note that we set v z to be parallel to v e . Hence, the integration of angle ϕ in eq. (B.3) can be done prior to θ, which gives a prefactor 2π. Finally, we can integrate I( v + v e ) over the angle θ in terms of variables a = v esc /v 0 , b = v e /v 0 , and x = v/v 0 . We find that I( v + v e ) has 3 types of solutions, depending on the relationships between a, b and x.

JHEP11(2017)024 C Confidence limits
To determine the exclusion limits, a confidence level (CL) is usually used, defined as where the desired CL value is the integral likelihood for the range 0 < χ 2 < χ 2 c . If the likelihood L is Gaussian distribution (∝ exp(−χ 2 /2)), then we can simply obtain χ 2 c = 2.71 at CL = 0.95 (one sided confidence limit). Using such a definition, conventionally we can describe the consistency of the data with the background-only hypothesis CL b and the signal-plus-background hypothesis CL s+b .
However, in many real cases one can introduce nuisance parameters and model parameters to the likelihood function as systematic uncertainties which makes it difficult to find CL b and CL s+b by integrating the multi-dimensional parameter space in eq. (C.1). Therefore, it is convenient to introduce a test-statistic χ 2 (s, b) = −2 ln L(s, b) to numerically pin down the confidence level for two hypotheses of background-only and signal-plus-background.
In the m χ − R plane, we divide m χ into various bins. For each m χ bin, the χ 2 varies as a function of R. Then, we solve the following equation to get a value of R 95 corresponding to the 95% upper limits based on the CL b method. Repeating the process for different m χ bins, we can get a curve on the m χ − R 95 plane corresponding to the 95% exclusion limit. Regarding to CL s+b method, the same procedure as described above can be performed, but solving a different equation to get the R 95 limits. Since we are using the null hypothesis (background only), we adopt the CL b method to evaluate the limits. It is worthy mentioning that our limits will be hence slightly weaker than the CL s+b method.

D Limit including PandaX run10 data
In this appendix, we will show the update limit including the latest PandaX run9 +run10 data [51]. For the new Run9 , the total observed event number after all the cuts is 1 and the expected background event number is 3.2 ± 0.9. For the Run10, there is none observed event number but the expected background event number is 1.8 ± 0.5. With the updated exposure E 9 = 26180.44 kg-day and E 10 = 27871.65 kg-day, we replace the old PandaX likelihood (run8+run9) to new one (run9 +run10) in the figure 11. One can see that the improvement is around 20% in the spin independent cross section but it is not significant in the (m χ , c p i ) plane.