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.


I. 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 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. Ref. [5][6][7] for the current status of the blind-spot region and Ref. [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 velocity-dependent framework, the published experimental limits for spin-independent and spin-dependent 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 modelindependent limit from these experiments will be very useful to link the WIMP models with the direct detection experiments. results of our scans. We summarize our findings in Sec. V.

II. 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].
The non-relativistic quantum mechanical operators are listed in the Table I Table I 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 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. (6) with the recoil energy Q in Eq. (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.  [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]  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. (9), 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. (9), 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 χ .

B. 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 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 3 The PandaX collaboration adopted the boosted-decision-tree (BDT) method to optimize the rejection of the nuclear recoil background, while the LUX group did not use the BDT. Furthermore, a time-dependent mapping between the true recoil position and the observed S2 coordinate is required to interpret the LUX data.
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. (12), we do not need to consider the correction factor f (m χ ) as in the PandaX likelihood.

C. 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 Fig. 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.

D. 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. Here, we adopted exactly the same astrophysical setup as used in the experimental computation, CL for PandaX, LUX, and XENON1T. In all the mass region, the XENON1T limits are 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 Fig. 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.

E. 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. (B4)).
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 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 II, 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.

IV. NUMERICAL RESULTS
In this section, the combined limits from the PandaX, LUX and XENON1T are pre- The astrophysical parameters, including the local velocity and escape velocity, are treated as nuisance parameters following Gaussian distributions with parameters given in Table II.
As we have discussed above, a larger uncertainty band of the local density measurements [53] is adopted to account for possible systematic uncertainties.
In many ISV studies [55][56][57][58][59][60][61][62][63][64][65][66], the ISV coupling ratio for maximum cancellation between the DM-proton contribution and the DM-neutron contribution is defined at the amplitude 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.  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 Φ ′′ .   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 Fig. 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 Fig. 4. The difference between these two is very small, suggesting that the value −N np /N nn gives an excellent    In previous sections, we have mentioned that all the four-point effective Lagrangians can be expanded by the NR effective operators shown in Table I. In fact, many relativistic Lagrangians correspond to one operator in the non-relativistic limit with a simple relationship (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. (15). 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 Fig. 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. (15), only the Lagrangian L 12 has a positive maximum cancellation ratio. Therefore 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 Fig. 8, and the maximum cancellation ratios of neutron to proton couplings are summarized in Fig. 9.

C. 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 long-distance 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 dimensionsix 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  In addition, we consider four dimension-six effective interactions between DM and quarks (q) 3,q = (χγ µ χ)(qγ µ γ 5 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 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 2,q , and the interaction Q 4,q has one additional relevant coefficient c p,n 4 added to that of Q 3,q . In Fig. 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 2 are stronger than that of the magnetic dipole interaction Q 1 . The dominant operators among the six operators in Q 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 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. 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. (B3) 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.
1. Region x + b < a or v + v e < v esc : 2. Region (x + b > a, |x − b| < a) or (v + v e > v esc , |v − v e | < v esc ): 3. Region |x − b| > a or |v − v e | > v esc : 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 Fig. 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.