Constraints on dimension-seven operators with a derivative in effective field theory for Dirac dark matter
CConstraints on dimension-seven operators with a derivative in effective field theory for Dirac darkmatter
Tong Li ∗ and Yi Liao
1, 2, † School of Physics, Nankai University, Tianjin 300071, China Center for High Energy Physics, Peking University, Beijing 100871, China
The effective field theory (EFT) for dark matter (DM) has been widely used to investigate dark matter detec-tion in both theoretical prediction and experimental analysis. To form a complete basis of effective operators forDirac DM EFT at dimension seven, eight new four-fermion operators with a derivative in DM currents have re-cently been introduced. We discuss the experimental observables and constraints for the theoretical predictionsof these new operators to constrain the DM mass and relevant energy scale. The observables from thermal relicabundance, indirect and direct detection, and LHC constraints are presented.
I. INTRODUCTION
Despite of substantial efforts in theory and experiment, themicroscopic properties of dark matter (DM) particles are stillunknown. Due to the plethora of competing theoretical mod-els in the current status, it is not feasible to extract the funda-mental properties of dark matter particles by contrasting the-oretical predictions with observation. To avoid this problem,the effective field theory (EFT) offers a good approach basedon general and minimal theoretical assumptions regarding thephysics underlying dark matter particles [1–4]. The EFT prin-ciple formulates specific dark matter models as a quantumfield theory, under the assumption that the DM candidate is asingle particle beyond the Standard Model (SM) and all otherdegrees of freedom are either heavy enough to be integratedout or have negligible strength for the observable spectrum.The simplest way to build an EFT for DM is to introducea new SM gauge singlet field, χ . We assume it to be a Diracfermion in this article. We assume further it is odd under anew parity while all SM fields are even, so that the χ parti-cle is guaranteed to be stable and can only be created or an-nihilated in pair. The interaction Lagrangian containing allLorentz and gauge invariant operators involving a pair of the χ field is schematically written as an expansion in the dimen-sion d of operators L χ = (cid:88) d,i,f C di O di,f , (1)where the Wilson coefficients are generically parameterizedas C di = Λ − di with Λ i being an effective cutoff scale for EFT.We will restrict ourselves in this work to the interactions ofDM with the quarks, and assume the coefficients are universalin quark flavor f = q . The framework includes two elec-tromagnetic dipole operators at dimension-5 (dim-5) and fourdim-6 operators formed from products of vector and axial-vector currents [3]. Those dim-6 operators have been adoptedto analyze the data of DM searches in indirect (ID) and directdetection (DD) experiments. For the unsuppressed vector op-erators, the DD measurements result in more stringent limits ∗ [email protected] † [email protected] than what can be expected from relic abundance [5]. More-over, the spin-1 mediator scenario via a vector or axial-vectorinteraction is highly constrained by the Z (cid:48) → dijet searchesat the Large Hadron Collider (LHC), and the mass of the(axial-)vector mediators has been excluded up to a scale of4-5 TeV [6–11]. The largely pushed energy scale for dim-6operators motivates us to consider the effects of higher dimen-sional operators. The higher dimensional operators could beinduced by underlying theories at perhaps a different energyscale from dim-6 operators, and thus may likely dominate theDM relic abundance when only one operator is switched on ata time. In fact, the detailed phenomenology of the six dim-7operators involving a scalar or tensor fermion current has beenwidely investigated [3, 12, 13]. For instance, for the scalaroperator scaled by quark mass, the DD experiments and theLHC have excluded the energy scale below about 1 TeV and100 GeV, respectively.The EFT Lagrangian with complete and independent dim-7operators describing a pair of the Dirac DM χ field interactingwith the quarks, gluons and photon is given by L dim − χ = (cid:88) i,q C i O i,q + L ( G aµν ) + L ( F µν ) , (2)where the first terms are four-fermion interactions whose Wil-son coefficients are assumed to be quark flavor universal, C i = Λ − i . The last two terms consist of eight operators thatcouple DM to a pair of gluon or photon field strength ten-sors. Among the four-fermion operators there are six com-monly considered ones, i.e., m q ¯ χ O χ χ ¯ q O q q with O χ,q ∈{ , iγ , σ µν } , which are suppressed for light quarks. Re-cently, it was pointed out that there exist additional eight four-fermion operators with a derivative acting on the DM fields asshown in Table I [14], where ¯ χi ←→ ∂ µ χ = ¯ χi∂ µ χ − ¯ χi ←− ∂ µ χ . Sincethe phenomenological studies of these new operators are stilllacking, we will fill this gap.In this work we confront the theoretical predictions of thesenew dim-7 operators with the experimental observables andconstraints to infer the most probable mass of DM and its in-teraction strengths with ordinary matter. We will also spec-ify distinct features of these operators against other opera-tors in DM searches. This paper is organized as follows. InSec. II, we study the experimental observables and constraintsfor these operators, and their allowed regions of the cutoff a r X i v : . [ h e p - ph ] A ug Operator O i,q Coefficient C i D : ¯ χi ←→ ∂ µ χ ¯ qγ µ q / Λ D : ¯ χiγ i ←→ ∂ µ χ ¯ qγ µ q / Λ D : ¯ χi ←→ ∂ µ χ ¯ qγ µ γ q / Λ D : ¯ χiγ i ←→ ∂ µ χ ¯ qγ µ γ q / Λ D : ∂ µ ( ¯ χσ µν χ )¯ qγ ν q / Λ D : ∂ µ ( ¯ χσ µν iγ χ )¯ qγ ν q / Λ D : ∂ µ ( ¯ χσ µν χ )¯ qγ ν γ q / Λ D : ∂ µ ( ¯ χσ µν iγ χ )¯ qγ ν γ q / Λ TABLE I. New four-fermion dim-7 operators and associated energyscales. scale and DM mass are given in Sec. III. Finally, we sum-marize our main results in Sec. IV.
II. OBSERVABLES AND CONSTRAINTS
In this section we discuss the experimental observables andconstraints for the new dim-7 operators with a derivative inDM currents, including thermal relic abundance, indirect anddirect detection, and LHC bounds.
A. Relic Density
The thermally averaged cross sections of DM pair annihi-lation into a quark pair through each new dim-7 operator at atime are, to the leading order in DM velocity v , (cid:104) σv (cid:105) D = N C m χ π Λ v (cid:88) q Θ( m χ − m q )(3 − β q ) β q , (3) (cid:104) σv (cid:105) D = N C m χ π Λ v (cid:88) q Θ( m χ − m q )(3 − β q ) β q = (cid:104) σv (cid:105) D (Λ → Λ ) , (4) (cid:104) σv (cid:105) D = N C m χ π Λ v (cid:88) q Θ( m χ − m q ) β q , (5) (cid:104) σv (cid:105) D = N C m χ π Λ v (cid:88) q Θ( m χ − m q ) β q = (cid:104) σv (cid:105) D (Λ → Λ ) , (6) (cid:104) σv (cid:105) D = 2 N C m χ π Λ (cid:88) q Θ( m χ − m q )(3 − β q ) β q , (7) (cid:104) σv (cid:105) D = 4 N C m χ π Λ (cid:88) q Θ( m χ − m q ) β q , (8)where β q = (cid:113) − m q /m χ , N C = 3 for quarks, and Θ is theHeaviside function for only taking into account on-shell two-body annihilations. We assume all kinematically accessiblequarks in the final state of the DM annihilation. Note that the annihilation rates for operators D and D are equivalentto those of D and D , respectively, when all of the fourfermions are on-shell by making repeated use of the equationsof motion. The velocity scaling for these annihilation crosssections is collected in the second column of Table II. Onecan see that they display various types of velocity dependencein the annihilation rate. The annihilation rates are d -wave foroperators D and D , and so are proportional to the fourthpower of DM velocity v . Operators D , D , D and D however have a p -wave term. Operator O i,q Anni. (cid:104) σv (cid:105) NR operator O Ni D O ( v ) 1D O ( v ) s χ · q D O ( v ) s N · v ⊥ D O ( v ) ( s χ · q )( s N · v ⊥ )D O (1) q , ( s χ · q )( s N · q ) , s χ · ( v ⊥ × q ) , q s χ · s N D O ( v ) s χ · q D O (1) s χ · ( s N × q )D O ( v ) ( s χ · q )( s N · v ⊥ ) TABLE II. Velocity scaling of annihilation cross sections (cid:104) σv (cid:105) andNR DM-nucleon operators O Ni for dim-7 operators considered inthis work. s χ ( s N ) is the DM (target nucleon) spin, and q and v ⊥ are scattering exchange momentum and velocity defined in Ref. [14]. The thermal DM relic abundance is determined by the equa-tion [15] Ω χ h = 1 . × GeV − M Pl x F √ g ∗ a + 3 b/x F + 20 c/x F , (9)for the expansion of annihilation cross section (cid:104) σv (cid:105) ∼ a + bv + cv . Here, M Pl ≈ . × GeV is the Planckmass, h is the Hubble parameter, g ∗ is the number of relativis-tic degrees of freedom, and T F is the freeze-out temperatureappearing in x F = m χ /T F . We vary x F and g ∗ in the rangeof < x F < [16, 17] and < g ∗ < [18], respec-tively, and adopt the relic abundance measured by Planck, i.e. Ω χ h = 0 . ± . [19]. Note that these choices arerather simplistic but are sufficient to estimate the relic densityin the context of EFT [5]. B. Indirect Detection
Dwarf galaxies are bright targets to search for DM anni-hilation through gamma rays. The Fermi Large Area Tele-scope (LAT) has searched for gamma ray emission from thedwarf spheroidal satellite galaxies (dSphs) of the Milky Waybut detected no excess. The Fermi-LAT thus set an upper limiton the DM annihilation cross section from a combined anal-ysis of multiple Milky Way dSphs [20, 21]. For individualdwarf galaxy targets, the Fermi-LAT collaboration tabulatedthe delta-log-likelihoods as a function of the energy flux bin-by-bin. The gamma ray energy flux from DM annihilation inthe j th energy bin is given by Φ Ej,k ( m χ , (cid:104) σv (cid:105) , J k ) = (cid:104) σv (cid:105) πm χ J k (cid:90) E max j E min j E dN γ dE dE, (10)where J k is the J factor for the k th dwarf and dN γ /dE de-scribes the gamma-ray spectrum from DM annihilation. Theenergy flux only depends on m χ , (cid:104) σv (cid:105) and J k , and is calcu-lable for DM annihilation processes given by the above EFToperators. The likelihood for the k th dwarf is L k ( m χ , (cid:104) σv (cid:105) , J k ) = L J ( J k | ¯ J k , σ k ) (cid:89) j L j,k (Φ Ej,k ( m χ , (cid:104) σv (cid:105) , J k )) , (11)where L j,k is the likelihood tabulated by the Fermi-LAT foreach dwarf and calculated gamma-ray flux and the uncertaintyof the J factors is taken into account by profiling over J k inthe likelihood below [20] L J ( J k | ¯ J k , σ k ) = 1ln(10) J k √ πσ k × exp (cid:20) − σ k (cid:16) log ( J k ) − log ( ¯ J k ) (cid:17) (cid:21) , (12)with the measured ¯ J k and error σ k . Then one can perform ajoint likelihood for all dwarfs L ( m χ , (cid:104) σv (cid:105) , J ) = (cid:89) k L k ( m χ , (cid:104) σv (cid:105) , J k ) , (13)where J is the set of J k factors. In our implementation weadopt the corresponding values of L j,k and ¯ J k , σ k for 19dwarf galaxies considered in Ref. [21].As Fermi-LAT found no gamma ray excess from the dSphs,for a given m χ , one can set an upper limit on the DM annihila-tion cross section by taking J factors as nuisance parametersin the maximum likelihood analysis. We follow Fermi’s ap-proach and take the delta-log-likelihood as below −
2∆ ln L ( m χ , (cid:104) σv (cid:105) ) = − L ( m χ , (cid:104) σv (cid:105) , (cid:98)(cid:98) J ) L ( m χ , (cid:100) (cid:104) σv (cid:105) , (cid:98) J ) , (14)where (cid:100) (cid:104) σv (cid:105) and (cid:98) J maximize the likelihood while (cid:98)(cid:98) J maximizesthe likelihood for given m χ and (cid:104) σv (cid:105) . The 95% C.L. upperlimit on the annihilation cross section for a given m χ is de-termined by demanding −
2∆ ln L ( m χ , (cid:104) σv (cid:105) ) ≤ . . Weobtain the spectrum of photons induced by annihilation intoquarks using the PPPC4DMID code [22] and perform the like-lihood analysis using Minuit [23]. Once the annihilation crosssection obtained from a certain set of m χ and Λ is larger thanthe limit, we claim the corresponding parameter values areexcluded by the Fermi-LAT dSphs measurement. Due to thesuppression by the extremely non-relativistic DM velocity, weexpect ID constraints to be relatively weaker for d- and p-waveoperators. C. Direct Detection
We show in the third column of Table II the non-relativistic(NR) operators of DM scattering off the nucleon induced fromthe considered dim-7 operators at quark level [14]. One cansee that, for operators D -D , the scattering rates are ei-ther suppressed by the spin of the target nucleus s N or thescattering momentum exchange q or both, rendering weakDD constraints. Only operator D leads to non-momentum-suppressed spin-independent (SI) DM-nucleon scattering, andis thus highly constrained by the direct DM detection. TheNR reduction of the operator D to the DM-nucleon levelis [14, 24] C D O D → C N D O N = 2 m χ (2 C D + C D ) O N = 6 m χ Λ O N , with O N = 1 χ N , (15)for the interaction of DM with the nucleon. The SI DM-nucleon scattering cross section is thus given by σ SI χN (D ) = µ χN π (cid:0) C N D (cid:1) = µ χN π (cid:18) m χ Λ (cid:19) = 4 . × − cm (cid:16) µ χN (cid:17) (cid:16) m χ
10 GeV (cid:17) (cid:18) (cid:19) , (16)where µ χN = m χ m N / ( m χ + m N ) is the reduced mass with m N being the nucleon mass. This prediction can then be com-pared directly to the limits set by DD experiments to yield alower bound on the cutoff scale Λ for a given m χ . In Fig. 1,we show the SI DM-proton scattering cross section versus DMmass for different values of Λ . For instance, for Λ = 1 TeV the whole range of m χ > GeV is excluded by Xenon1T [25, 26], while DM with m χ < GeV can evade theDD limit when Λ = 10 TeV.
FIG. 1. SI DM-proton scattering cross section versus DM mass foroperator D at Λ = 10 , TeV is compared with Xenon1T lim-its at 90% CL (solid black [25] and dashed black [26]) and neutrinobackground (green).
D. LHC Constraints
The LHC constraints on DM EFT stem from searches forlarge missing energy events produced alongside with a vis-ible object such as a jet, lepton, or photon, i.e. the so-called mono-X searches. The most stringent constraint forthe operators we consider comes from the mono-jet searchcorresponding to an integrated luminosity of 36.1 fb − at acentre-of-mass energy of 13 TeV [27]. In order to estimatethe mono-jet constraint on our EFT setups, we create UFOmodel files using FeynRules [28] and interface them withMadGraph5 aMC@NLO [29] to generate signal events com-posed of DM pairs with a jet from initial-state radiation. Thesignal events are then passed to Pythia [30] and Delphes [31]for parton shower and detector simulation, respectively. Fol-lowing the event selection in Ref. [27], we require the leadingjet satisfying p T > GeV and | η | < . and the missingtransverse momentum with E miss T > GeV. Ref. [27] pro-vides the observed 95% confidence level (CL) upper limit onthe visible cross section, defined as the product of cross sec-tion and efficiency corresponding to the above selection cuts.Once the visible cross section obtained from a certain value of m χ and Λ is larger than the limit, we claim the correspondingparameters are excluded by the mono-jet search at 95% CL.In Fig. 2 we compare the normalized distributions in thetransverse momentum of the leading jet for the D operator(black), a dim-6 operator ¯ χγ µ χ ¯ qγ µ q (red), and a dim-7 oper-ator m q ¯ χχ ¯ qq (green), assuming m χ = 100 GeV at 13 TeVLHC. The signal distribution of the new dim-7 operator witha derivative does not decrease as fast as the other two in thehigh energy region due to the derivative enhancement [32]. (GeV) T leading jet p400 600 800 1000 1200 1400 -3 -2 -1
10 D15 q m g q c m gc qq cc q mD15 q m g q c m gc qq cc q mD15 q m g q c m gc qq cc q m FIG. 2. Normalized leading jet p T for D operator (black solid),dimension-6 operator ¯ χγ µ χ ¯ qγ µ q (red dashed) and dimension-7 op-erator m q ¯ χχ ¯ qq (green dotted), assuming m χ = 100 GeV.
III. RESULTS
In this section, for each operator in Table I we show itsallowed region in the cut-off scale Λ i versus m χ by the ob-servables discussed in the previous section. Figs. 3, 4, and5 correspond to the operators resulting in d -, p - and s -waveannihilation rates, respectively. The correct thermal DM relic abundance with correspond-ing Λ and m χ is given by the red band. The band is derivedfrom Eq. (9) and reflects the assumed ranges of values for x F and g ∗ . As expected before, the more the annihilation rate issuppressed, the weaker the ID constraint such as the exclusionby Fermi-LAT dSphs becomes. Most severely, as indicated bythe blue squares, Fermi-LAT excludes a majority of space be-low Λ i (cid:39) TeV for operators D and D .The mono-jet search at 13 TeV LHC excludes the param-eter space to the left of the orange solid line that essentiallyamounts to a lower limit on the EFT scale Λ i (cid:46) TeV. ThisLHC constraint is more sensitive to the low m χ region, thuscomplementary to the indirect detection. Besides, the limit ofscattering cross section from Xenon 1T at 90% CL severelyconstrains the operator D such that only the blank band inthe left panel of Fig. 3 remains to be explored by future directdetection experiments.Finally, the EFT approximation is valid above the black dot-ted lines, i.e. roughly for Λ i > m χ / (2 π ) . The region yieldinga correct relic density is compatible with EFT validity. Largercouplings will violate perturbative unitarity whence the EFTexpansion breaks down and cannot give a reliable descriptionof an underlying theory.In summary, to avoid overproduction of DM and ensure thevalidity of the EFT approximation, the viable Λ i – m χ regionmust fall between the red band and the dotted line. Further-more, the energy scale Λ i has to be greater than about 1 TeVto satisfy the LHC bound. In particular, this squeezed regionis entirely excluded by direct detection for D and mostlyexcluded by indirect detection for D and D . (a) (b) FIG. 3. The allowed region of cut-off scale Λ vs. m χ by Planck(red band) for operators D (a) and D (b). EFT is valid above thedashed line. The blue region is excluded by the null measurementof dwarf galaxies by Fermi-LAT. The excluded region by Xenon1T(purple) and the region below the neutrino background (green) arealso shown for operator D . The orange curve represents the LHCbound. IV. CONCLUSION
In this work we have investigated new dimension-7 opera-tors in effective field theory for Dirac fermionic dark matter.These operators involve a derivative in the DM currents and (a) (b)
FIG. 4. Results for operators D (D ) and D (D ), as labeled inFig. 3. (a) (b) FIG. 5. Results for operators D (a) and D (b), as labeled inFig. 3. their phenomenology has not yet been studied in the literature.We discussed the experimental observables and constraints forthese operators to confine the DM mass and relevant energyscale. We found that these operators induce various s -, p - and d -wave annihilation rates and are thus, to different extents,constrained by indirect DM detection such as the Fermi-LATdSphs. In spite of this, the correct thermal relic abundancecan be achieved in the parameter space allowed by indirectdetection. The mono-jet search at 13 TeV LHC excludes theparameter space with energy scale Λ i (cid:46) TeV and m χ (cid:46) TeV. And only one of the operators gives non-momentum-suppressed spin-independent DM-nucleon scattering, and isthus highly constrained by direct detection experiments.
ACKNOWLEDGMENTS
TL thanks Thomas Jacques for helpful discussions. Thiswork was supported in part by the Grants No. NSFC-11575089 and No. NSFC-11025525, by The National KeyResearch and Development Program of China under GrantNo. 2017YFA0402200, and by the CAS Center for Excel-lence in Particle Physics (CCEPP). T.L. is supported by “theFundamental Research Funds for the Central Universities”,Nankai University (Grant Number 63191522, 63196013). [1] M. Beltran, D. Hooper, E. W. Kolb and Z. C. Krusberg, Phys.Rev. D , 043509 (2009) doi:10.1103/PhysRevD.80.043509[arXiv:0808.3384 [hep-ph]].[2] J. Fan, M. Reece and L. T. Wang, JCAP , 042 (2010)doi:10.1088/1475-7516/2010/11/042 [arXiv:1008.1591 [hep-ph]].[3] J. Goodman, M. Ibe, A. Rajaraman, W. Shepherd,T. M. P. Tait and H. B. Yu, Nucl. Phys. B , 55 (2011)doi:10.1016/j.nuclphysb.2010.10.022 [arXiv:1009.0008 [hep-ph]].[4] C. Bal´azs, T. Li and J. L. Newstead, JHEP , 061 (2014)doi:10.1007/JHEP08(2014)061 [arXiv:1403.5829 [hep-ph]].[5] C. Bal´azs, J. Conrad, B. Farmer, T. Jacques, T. Li,M. Meyer, F. S. Queiroz and M. A. S´anchez-Conde, Phys.Rev. D , 083002 (2017) doi:10.1103/PhysRevD.96.083002[arXiv:1706.01505 [astro-ph.HE]].[6] V. Khachatryan et al. [CMS Collaboration], Phys. Rev. Lett. , 031802 (2016) doi:10.1103/PhysRevLett.117.031802[arXiv:1604.08907 [hep-ex]].[7] A. M. Sirunyan et al. [CMS Collaboration], Phys. Lett.B , 520 (2017) doi:10.1016/j.physletb.2017.02.012[arXiv:1611.03568 [hep-ex]].[8] A. M. Sirunyan et al. [CMS Collaboration], JHEP ,097 (2018) doi:10.1007/JHEP01(2018)097 [arXiv:1710.00159[hep-ex]]. [9] A. M. Sirunyan et al. [CMS Collaboration], arXiv:1803.08030[hep-ex].[10] M. Aaboud et al. [ATLAS Collaboration], Phys. Rev.D , 052004 (2017) doi:10.1103/PhysRevD.96.052004[arXiv:1703.09127 [hep-ex]].[11] M. Aaboud et al. [ATLAS Collaboration], arXiv:1801.08769[hep-ex].[12] G. D’Ambrosio, G. F. Giudice, G. Isidori and A. Strumia, Nucl.Phys. B , 155 (2002) doi:10.1016/S0550-3213(02)00836-2[hep-ph/0207036].[13] J. Goodman, M. Ibe, A. Rajaraman, W. Shepherd,T. M. P. Tait and H. B. Yu, Phys. Rev. D , 116010(2010) doi:10.1103/PhysRevD.82.116010 [arXiv:1008.1783[hep-ph]].[14] J. Brod, A. Gootjes-Dreesbach, M. Tammaro and J. Zupan,arXiv:1710.10218 [hep-ph].[15] G. B. Gelmini, P. Gondolo and E. Roulet, Nucl. Phys. B ,623 (1991). doi:10.1016/S0550-3213(05)80036-7[16] G. Jungman, M. Kamionkowski and K. Griest, Phys. Rept. , 195 (1996) doi:10.1016/0370-1573(95)00058-5 [hep-ph/9506380].[17] Q. H. Cao, C. R. Chen, C. S. Li and H. Zhang, JHEP , 018(2011) doi:10.1007/JHEP08(2011)018 [arXiv:0912.4511 [hep-ph]].[18] T. Schafer, hep-ph/0304281. [19] P. A. R. Ade et al. [Planck Collaboration], Astron. As-trophys. , A13 (2016) doi:10.1051/0004-6361/201525830[arXiv:1502.01589 [astro-ph.CO]].[20] M. Ackermann et al. [Fermi-LAT Collabora-tion], Phys. Rev. Lett. , 231301 (2015)doi:10.1103/PhysRevLett.115.231301 [arXiv:1503.02641[astro-ph.HE]].[21] A. Albert et al. [Fermi-LAT and DES Collaborations], As-trophys. J. , 110 (2017) doi:10.3847/1538-4357/834/2/110[arXiv:1611.03184 [astro-ph.HE]].[22] M. Cirelli et al. , JCAP , 051 (2011) Erratum: [JCAP , E01 (2012)] doi:10.1088/1475-7516/2012/10/E01,10.1088/1475-7516/2011/03/051 [arXiv:1012.4515 [hep-ph]].[23] F. James and M. Roos, Comput. Phys. Commun. , 343(1975). doi:10.1016/0010-4655(75)90039-9[24] A. De Simone and T. Jacques, Eur. Phys. J. C , 367 (2016)doi:10.1140/epjc/s10052-016-4208-4 [arXiv:1603.08002 [hep-ph]]. [25] E. Aprile et al. [XENON Collaboration], Phys. Rev. Lett. , 181301 (2017) doi:10.1103/PhysRevLett.119.181301[arXiv:1705.06655 [astro-ph.CO]].[26] E. Aprile et al. [XENON Collaboration], arXiv:1805.12562[astro-ph.CO].[27] M. Aaboud et al. [ATLAS Collaboration], JHEP ,126 (2018) doi:10.1007/JHEP01(2018)126 [arXiv:1711.03301[hep-ex]].[28] A. Alloul, N. D. Christensen, C. Degrande, C. Duhrand B. Fuks, Comput. Phys. Commun. , 2250 (2014)doi:10.1016/j.cpc.2014.04.012 [arXiv:1310.1921 [hep-ph]].[29] J. Alwall et al. , JHEP , 079 (2014)doi:10.1007/JHEP07(2014)079 [arXiv:1405.0301 [hep-ph]].[30] T. Sjostrand, S. Mrenna and P. Z. Skands, JHEP ,026 (2006) doi:10.1088/1126-6708/2006/05/026 [hep-ph/0603175].[31] J. de Favereau et al. [DELPHES 3 Collaboration], JHEP ,057 (2014) doi:10.1007/JHEP02(2014)057 [arXiv:1307.6346[hep-ex]].[32] Y. Cai, M. A. Schmidt and G. Valencia, JHEP1805