Interplay of the LHC and non-LHC Dark Matter searches in the Effective Field Theory approach

We present accurate and up-to-date constraints on the complete set of dimension five and six operators with scalar, fermion and vector Dark Matter (DM). We find limits using LHC mono-jet data, spin inde- pendent and spin dependent direct searches, relic density and CMB, and show the interplay between high and low energy data in setting bounds on the parameter space. In order to properly compare data taken at different energies, we take into account the effect of the running and mixing of operators. We also take into account the local density uncertainties affecting direct detection data, and apply EFT validity criteria related to the cut on the invariant mass of DM pair production at the LHC, which turns out to be especially important for the case of vector DM. Finally, we estimate the potential of the future LHC runs to probe DM parameter space.


Introduction
Understanding the nature of Dark Matter (DM) is one of the greatest puzzles of modern particle physics and cosmology. Although overwhelming observational evidences from galactic to cosmological scales point to the existence of DM [1][2][3], after decades of experimental effort only its gravitational interaction has been experimentally confirmed. Currently, no information is available on the DM properties, such as its spin, mass, interactions other than gravitational, symmetry responsible for its stability, number of states associated to it, and possible particles that would mediate the interactions between DM and the standard model (SM) particles.
If DM is light enough and interacts with SM particles directly or via some mediators with a strength beyond the gravitational one, its elusive nature can be detected or constrained in different ways: • from direct production at colliders, resulting in a signature exhibiting an observed SM object, such as jet, Higgs, Z, or photon, that recoils against the missing energy from the DM pair [4][5][6][7]; • via the relic density constraint obtained through the observations of cosmic microwave background (CMB) anisotropies, such as those of WMAP and PLANCK collaborations [1,8]; • from DM direct detection (DD) experiments, which are sensitive to elastic spin independent (SI) or spin dependent (SD) DM scattering off nuclei [9][10][11][12]; • from DM indirect detection searches, that look for SM particles produced in the decay or annihilation of DM present in the cosmos, both with high energies observables (gamma-rays, neutrinos, charge cosmic rays) produced in the local Universe [13][14][15][16][17][18], and by studying the effects of energy produced by DM annihilation in the early universe on the properties of the CMB spectrum [1,19,20].
In this work we obtain the present constraints on three scenarios for the DM particles: complex scalars (φ), Dirac fermions (χ) and complex vectors (V µ ). In order to describe the interactions of the new states, we parametrize the DM interactions with the SM quarks and gluons as an effective field theory (EFT) that contains a complete set of operators of dimension six or less. Notice the presence of the coupling g * in the definition of the effective operators, which we insert according to the Naive Dimensional Analysis [21]. Moreover, for the vector DM case we choose the parametrisation suggested in Ref. [22] that takes into account the high energy behaviour of the scattering amplitudes that are enhanced by an energy factor (E/m DM ) for every longitudinal vector DM polarisation.
Here, our goal is to explore the complementarity of the collider and non-collider experiments mentioned above for all DM EFT operators and DM spin listed in Table 1 [22][23][24] for DM masses in the GeV-TeV range; for lighter DM masses see for example [25]. For the sake of consistency of our analyses, we obtain the present constraints from the LHC data taking into account the validity of the EFT [26][27][28] using the prescription of Ref. [29]; see Section 3 for details. Furthermore, the correct comparison between the LHC and non-collider bounds requires that we consider in our analyses the running of the EFT operators from the TeV scale down to the GeV one. This is important because the running of operators leads to mixing between them at low energy which can give rise to stronger DM DD limits [30][31][32][33][34][35][36].
As it is well known, the DM DD searches are plagued by the uncertainty on the local DM density [37][38][39][40][41][42] which propagates to the limits on the DM DD cross sections reported by the experiments, inducing variations up to one order of magnitude. Here we estimate the impact of this uncertainty on the DM bounds, presenting three scenarios that range from a conservative to a more optimistic one. On the other hand, indirect constraints from CMB are unaffected by the usual unknowns related to the DM density profile within structures. For this reason, we also include the CMB data in our analyses that lead to more robust limits. This paper is organized as follows: in Section 2 we study the running of the EFT DM operators from the TeV to the GeV scale and the effect of their mixing, while we present our analyses framework and available constraints in Section 3. Section 4 contains our main results for all operators in Table 1 that show the complementarity of the collider and non-collider constraints. Finally, we draw our conclusions in Section 5.

Direct detection and Running effect of the EFT operators
In this section we demonstrate the importance of the running of the operators for the DM DD constraints. As we know, in Quantum Field Theory radiative corrections may be important to properly assess the phenomenological implications of a model. In the case of DM EFT this is even more so, because (i) the bounds on the parameter space involve experiments with very different energy scales, and (ii) the Wilson coefficients can vary substantially between the typical LHC energies (a few TeV) and the typical energies of DM DD experiments (below the GeV).  The first radiative effects to be considered were the QCD ones [30][31][32], while only later it has been realized that EW loops could be important as well [33][34][35][36]43]. In our discussion we will make use of the renormalisation of currents [32], which we found fits best the purpose of our study. As can be seen from Table 1, all the operators are written as a product where A denotes some combination of Lorentz indexes and J χ and J SM are currents constructed out of DM and SM fields only. Since we suppose χ (and hence J χ ) to be a gauge singlet, only the renormalisation of the SM currents has to be computed. The lowest dimensional currents considered in this work are the dimension 3 and 4 operators 1 where q denotes any of the SM quarks and G is the gluon field. At the low energies involved in the DM DD experiments, the relevant degrees of freedom are nucleons and nuclei, rather than quarks and gluons. For this reason, we should match the amplitudes involving quarks and gluons with matrix elements involving nucleons (see for example [44,45] and Table 2 for some of the operators). Let us start with the scalar quark current. The matrix elements at zero momentum transfer are [46,47] for light (q = u, d, s) and heavy quarks (Q = c, b, t), respectively. In the above equation, N stands either for proton p or neutron n. The quantity f (N ) T q amounts to the light quark contribution to the nucleon mass m N , while f (N ) T q . Notice that both the quark condensate and the running mass depend on the scale µ, in such a way that their product is scale independent. As a consequence, the form factors are also scale independent. The light quarks form factors are known from hadron spectroscopy and lattice calculations. In our numerical analysis we use the quark scalar form factors presented in micrOMEGAs [48]. Using the matrix elements of Eq. (3) the DM-nucleon spin independent (SI) cross section can be written as where λ is an operator dependent coefficient that will be defined below.
To discuss the importance of the running, let us discuss the example of the C1 operator. We will suppose that such operator is generated at the scale Λ ∼ µ LHC ∼ TeV with Wilson coefficient g 2 * . As already explained, the qq operator runs as the inverse of the quark mass, in such a way that the C1 operator at an arbitrary scale µ looks like The same is true for the D1 and V1 operators, and the λ coefficient of Eq. (4) reads where N = n, p.
A similar effect is present for the operators involving the G µν G µν current (C5 and V11). In this case, the matrix element of the gluon current is [46,47] and the combination α s (µ)G µν G µν (and hence the form factor) is scale independent. Suppose now that the C5 operator is generated at the Λ scale with Wilson coefficient g 2 * . At an arbitrary scale µ, the operator looks like The same result applies to the V11 operator. The λ coefficients of Eq. (4) reads In principle, in addition to the running of the individual Wilson coefficients, a mixing between operators is generated [30][31][32][33]. However, for operators involving the {m q qq, G µν G µν } currents the mixing was found to be numerically unimportant [33].
Let us now turn to the analysis of the {qγ µ q, qγ µ γ 5 q} currents. In this case, the mixing between operators can be numerically important. Let us consider for instance the operators involving the qγ µ γ 5 q axial vector current, i.e. the operators C4, D7-D8, V4. These operators are responsible for spin dependent DM scattering at DD experiments, with bounds much weaker than those of spin independent experiments. Nonetheless, the qγ µ γ 5 q axial vector current mixes with the qγ µ q vector current during the running, and a SI cross section is generated [33][34][35]. Suppose for instance the C4 operator is generated at the Λ scale above the top mass where the sum is taken over all light (q) and heavy (Q) quark flavors. Then, using the results of References [33][34][35], the operators present in the Lagrangian at the DD scale are where we show only the most relevant contribution, coming from top loops. In the previous expression, α t ≡ y 2 t /4π, with y t the top Yukawa coupling. Notice that the top contribution is present only down to the top scale, where the top quark is integrated out. As already stressed, in addition to the SD operators of the first two lines, the running has generated the SI operators of the last two lines. Even though the Wilson coefficients of the SI operators are smaller than those of the SD ones, it has been shown in Reference [35] that they are sufficient to put bounds on g 2 * /Λ 2 which are up to a factor of 100 stronger with respect to the typical bounds that can be obtained from SD experiments (i.e. considering only the operators in the first two lines). This shows clearly the importance of the running in setting consistent bounds on the parameter space of the DM EFT. Let us point out however that for the effect to be numerically relevant, it is instrumental for the coupling between the DM and the top currents to be switched on in Eq. (11). If this is not the case, the SI operators are still generated in the running, but with much smaller Wilson coefficients and weaker bounds (see Reference [35] for more details). In the analysis of Sec. 4 we will numerically implement the running and the mixing of the currents using the runDM code [33][34][35]49], which takes into account all the contributions.

Analysis setup and constraints
In this section we describe the analysis setup and constraints used in this study. In particular we delineate the limits originating from CMB, direct detection experiments and collider searches.
Direct and indirect detection constraints are affected by uncertainties of astrophysical nature. On one hand, the scattering of DM off nuclei on the Earth depends on the DM local density and velocity distributions around Earth. On the other hand, the DM self-annihilation rate in our galaxy depends on its particle density distribution therein. For what is concern of this paper, whenever possible we make the conservative choice to select targets that can reduce as much as possible the uncertainties, and thoroughly account for the remaining ones. In practice, this means that: (i) for indirect searches we adopt CMB limits, as the energy injection of DM annihilation is unaffected by the usual unknowns related to DM density profile within structures; (ii) for direct searches we explicitly take into account the systematic effects generated by the astrophysical uncertainties in the determination of the local DM density ρ 0 .

CMB constraints
The observation of byproducts of DM annihilation (or decay) in astrophysical targets can be used to determine (or constrain, in case of missing observations) relevant DM properties such as its mass and annihilation cross section. Such bounds depend on the unknown DM distribution within the astrophysical objects chosen as targets. A detailed analysis has shown that, choosing the CMB as target, the leading signal of DM annihilation is produced around redshift z ∼ 600 [50]. This makes the CMB a quantitatively competitive channel for indirect searches [1,19,20], since at z ∼ 600 the DM has not fallen into structures yet, and the observation is free from the usual astrophysical sources of uncertainties (density profile within a halo, distribution and density of subhalos, mass of the smallest bound halo). Moreover, additional sources of systematics affecting the CMB constraints have also been thoroughly examined [51,52], and shown to affect the results below the sensitivity needed for this paper. In the following, we will neglect them leaving our conclusions unaltered.
In order to set the CMB bounds on the quantity of our interest Λ, we first obtain the observational bound on the thermally averaged annihilation cross-section at redshift 600: where σv j (600) is the thermally averaged partial annihilation cross section for the j-th channel at redshift 600 and f j (z, m DM ) is the fraction of annihilation energy that is absorbed by the plasma at redshift z. The quantity p ann is constrained by Planck TT, TE, EE and lowP data [1]: p ann < 4.1 × 10 −28 cm 3 s GeV at 95% C.L., (13) and the values of the variable f j (z, m DM ) are taken form Ref. [53]. In order to bound the new physics scale Λ, we numerically compute the velocity dependent annihilation cross section with micrOMEGAs [44], obtaining the relationship between σv j (600) and Λ for each effective operator and each final state. It may be noted that for those operators for which the s-wave process is dominating, the thermally averaged cross section is constant 2 and σv j (600)= σv j (0). On the other hand, for those operators in which the p-wave contribution dominates the annihilation cross section, the CMB bound is almost ineffective, since it is suppressed by the low T CM B temperature. A possible bound in this case can be obtained from the Big Bang Nucleosynthesis (BBN), since this process happens at T BBN ∼ MeV. As shown in References [54,55], the bounds on σv are generically weaker than those obtained from the CMB, and we have explicitly checked that in all the cases in which the CMB bound is ineffective, the BBN bound is also not relevant.

DM direct detection constraints
The determination of the DM mass and elastic scattering cross section in a DM DD experiments is affected by uncertainties associated to the flux of DM particles crossing the Earth at any given time. Although the uncertainties on the Sun's relative motion with respect to the Galactic Center, or the exact morphology of the Galactic bulge, do not affect the conclusion of the presence of a sizable component of DM at the Sun's location [37], they impact the reconstruction of the DM profile throughout the Milky Way [38]. More specifically, two sources of uncertainties are relevant for DM DD experiments: the local DM distribution 3 and its velocity structure [40][41][42]56]. Although relevant, their effect is often overlooked in putting DD bounds on the parameter space of a model. Recently, Ref. [57] has shown that the effect of known Galactic uncertainties on the reconstruction of particle physics parameters may overcome those on the velocity structure. We thus decide to follow closely the analysis therein and take into account the different determinations of the local DM density due to the variation of the galactic parameters. In particular, we consider three possible values ρ 0 = 0.06 GeV/cm 3 , ρ 0 = 0.3 GeV/cm 3 and ρ 0 = 1.8 GeV/cm 3 (corresponding to the lowest and largest possible ρ 0 , and to the one used by the experimental collaborations), and apply the experimental bounds stemming from Xenon1T SI searches [10] and from PandaX-II SD analyses [58].
In addition, when we obtain DD bounds, we also take into account the relic density constraint. Direct detection experiments constrain DM particles assuming that their relic density matches the one of the cold DM component. The simplest way to compute the detection rates is to rescale the DM density distribution according to the prescription where ρ 0 is the local DM density, Ω χ is the theoretical relic density obtained via micrOMEGAs for every operator listed in Table 1 and Ω DM 0.12 from Ref. [1].

Collider constraints
In the upcoming sections, we also present bounds coming from collider searches. In order to perform the required simulations, we implemented the different effective operators independently in FeynRules [59] and LanHep [60] and generated the signal using MadGraph5 aMC@NLO [61]. The hadronization and parton showering was done using PYTHIA 6.4 [62], with subsequent detector simulation performed using MadAnalysis5 [63] and Delphes [64]. In order to obtain the limits on the scale Λ we consider the CMS analysis of final state presenting jets and missing transverse energy [5] based on data obtained at 13 TeV with an integrated luminosity of 35.9 fb −1 . This analysis was performed as a counting experiment in 22 independent signal regions characterized by (i) E miss T > 250 GeV, (ii) one jet with p jet T > 100 GeV and (iii) |η j | < 4.5. In order to simulate the DM contribution to this process we studied for all operators in Table 1. As it is well known, higher dimensional operators such as the ones in our EFT can lead to perturbative partial wave unitarity violation at high energies, signaling a maximum value of the center of mass energy for its applicability. Therefore, in order to guarantee the validity of the EFT we impose in our simulation that the invariant mass of the DM pair M χχ,φφ,V V satisfies [29]: 2 We remind the reader that the thermally averaged cross section is given by with T the temperature at which the process is computed. 3 See Ref. [39] for a recent and thorough review.

8
In our statistical analysis we use the simplified likelihood approach given in Ref. [65], taking into account the full correlation and covariance matrix provided in [5]. More specifically, we defined the likelihood function where the s i (Λ, g * ) is the expected number of events of the DM signal in i th bin, b i is the respective number of background events and n i is the number of observed events. In our case the signal cross section for each bin is a function of the coupling g * and the scale Λ. For practical purposes we consider three different benchmark values for g * = 4π, 6 and 1. The systematic uncertainties of the SM backgrounds and the DM signal are treated as nuisance parameters and they are approximated by zero-mean Gaussian variables θ i and a covariance matrix V . We define our test statistic function as whereθ Λ is a θ vector that minimizes the logarithm of Eq. (17) for a given value of Λ. On the other hand,Λ andθ are the values of Λ and of the θ vector that globally minimize the logarithm of Eq. (17). We find the upper limit on the scale of the mediator varying Λ until T S(Λ) = 4. 4 In addition to the present limits we also perform a projection of the CMS reach for an integrated luminosity of 300 fb −1 . For this projection we assume a gaussian likelihood, that the number of background events scales with the luminosity and that the uncertainty on the background scales as the square root of the luminosity. We set a lower background limit to be 1% of the background, based on post-fit numbers with respective background error provided by ATLAS and CMS [5,66].

Non-collider constraints
In this section we present combined non-collider constraints for all operators under study. In particular, we obtain the bounds originating from DM DD searches, indirect DM searches from the CMB and relic density assuming the freeze-out mechanism and standard cosmology. The results in this section are obtained with three different tools: micrOMEGAs [48], a modified version of the code released with [45] and runDM [33][34][35]49].

Complex scalar DM
We start presenting in Fig. 1 the results for complex scalar DM. In all the panels, the area below the three blue curves represents the region excluded by SI DM DD experiments, taking into account the uncertainty on the local DM density as discussed in Section 3.2. The shaded blue region is the conservative exclusion, while the middle and upper contours correspond to the center (namely the one used by the experimental collaboration) and the most optimistic exclusions, respectively. The dashed yellow line corresponds instead to the region of parameter space in which the predicted relic density matches the DM density observed by the Planck collaboration. Above the dashed yellow curve, the predicted relic density is larger than the experimental one. The computation is done assuming that the observed DM relic density is associated with only one particle, interacting with the SM via only one of the effective operators listed in Table 1. The region excluded by the CMB measurement is indicated by a shaded green area below the solid green line. Finally, the grey-shaded region represents the region Λ < 2m DM , excluded by the requirement of the validity of the EFT.
From Fig. 1(a) we can learn that the C1 operator is strongly constrained by SI DM DD searches since this operator leads to a cross section that is neither velocity nor transferred-momentum suppressed. In fact, we can see from this panel that any scale Λ below 1000 TeV is excluded even for m DM = 1 TeV. Moreover, the other two blue contours associated with the astrophysical uncertainties show that the bounds can be tightened by almost an order of magnitude. Therefore, it is important to be conservative to make a robust exclusion of the parameter space, as we do in the present work. We can also see from Fig. 1(a) that the relic density data exclude Λ above 20 TeV for all DM mass range under consideration, since in this region the φ particle would over-close the universe. Notwithstanding, one should note that this bound is quite model dependent, e.g. adding an      additional particle that could co-annihilate with the DM could change the relic density dramatically. Finally, the CMB measurement leads to quite strong bounds that exclude Λ below 1-10 TeV for g = 1, depending on the DM mass. This is due to the fact that the C1 operator gives rise to s-wave annihilation cross section, which is strongly constrained by the CMB bound.
We show in Fig. 1(b) the non-collider constraints on the operator C2. As expected, the CMB and relic density bounds on C2 are identical to those on C1. However, the SI DM DD constraints are absent, since this Ωχ ∼ ΩDM CMB operator violates parity. Therefore, the limit from CMB plays a crucial role for this operator, excluding the parameter space below the solid green line for g = 1. The unconstrained window can be tightened for DM masses above 20 GeV using the relic density information if we assume that there is no process that leads to DM co-annihilation.
We display in Fig. 1(c) the limits for the operator C3, that contains a vector quark current. The bounds coming from SI DM DD are weaker than those on the operator C1 since C3 is a dimension six operator, while C1 is dimension five. The uncertainty on this constraint is roughly a factor of 3 due to the higher dimensionality of the operator. Moreover, the annihilation cross section for the C3 operator is velocity suppressed (p-wave), therefore the effect on the CMB spectrum is negligible. The complementarity of the relic density and DD DM bounds completely rule out WIMP models that give rise at low energy to this operators only.
The bounds on the C4 operator, that contains a pseudo-vector quark current, are shown in Fig. 1(d). At variance with the operator C2, the SI DM DD limits are non-vanishing for operator C4 since the running and mixing effects play an important role; for further details see Section 2. Notice that the sudden drop of the SI DM DD constraint around m DM 200 GeV is due to the rescaling of the DM density distribution given in Eq. (14). Similarly to the operator C3, C4 is also velocity suppressed, and consequently it is not constrained by the CMB.
Finally, Figs. 1(e) and 1(f) present the non-collider bounds for the C5 and C6 operators, involving the gluon field strength. As we can see, there is no DD bound on C6 due to its parity violating nature. The bound on C5 is instead important, given the large gluon content of the proton. For DM masses below 10 GeV, the bound is comparable to the one on the C3 operator. For larger masses, the bound weakens because the SI cross section is suppressed by an additional factor of m 2 DM with respect to the C3 case (see Eq. (9)). Still, the SI bound is always the dominant one, in the region in which the EFT is valid. The constraints of the operator C6 are instead much weaker, since it does not contribute to SI DM DD, being dominated by the CMB data and relic density.
From Fig. 1 we can notice that the slope of the CMB constraint is negative for the C1 and C2 operators, in   contrast with it being positive for C5 and C6. This happens because the leading term for the annihilation cross section is independent of the DM mass for C1 and C2, while it is proportional to the DM mass squared for C5 and C6. This picture is generic -if the amplitude for DM pair annihilation is not velocity suppressed, then using dimensional analysis one finds that: a) for dimension five operators, the annihilation cross section is constant with respect to the DM mass, leading to the slight negative slope because the "effective" photon spectrum for the CMB constraint drops with the invariant mass; b) for dimension six operators the annihilation cross section is proportional to M 2 DM and one observes positive slope in (M DM , Λ) plane for the CMB constraint.

Dirac fermion DM
We present in Figure 2 the non-collider constraints for the operators D1-D4, that exhibit scalar or pseudoscalar interactions of fermion Dirac DM with quarks. From this figure we can see the operator D1 is strongly constrained by the SI DM DD searches while the operators D2, D3 and D4 are not bound by these searches since they give rise either to spin-dependent interactions (D3 and D4) or the spin-independent cross section is momentum suppressed (D2). Notice that the uncertainty of the SI DM DD limits on D1 is of the order of 3 since this operator is of dimension six; this is true for all operators D1-D10. On the other hand, the CMB data lead to constraints on the operators D2 and D4 since they exhibit an s-wave annihilation cross section, while D1 and D3 are not bound by CMB due to their p-wave cross sections. The relic density constraint set a strong upper bound on the operators D2, D3 and D4 and rule out the operator D1 unless there is a mechanism to avoid the overproduction of DM in the early universe. As before, the theoretical consistency of our framework, i.e. the EFT validity region, becomes important just for large DM masses. We show in Fig. 3 the non-collider limits on the fermion operators D5-D8, which contain vector and axial vector currents. The SI DM DD bound on the D5 operator is very strong, since its SI cross section is unsup-   Figure 1 for the D9-D10 fermion operators involving tensor currents.
pressed. On the other hand, the constraints on the operator D6 are milder, since its SI cross section is both velocity and nuclear recoil momentum suppressed. The operator D7 has in principle only SD cross sections. However, as discussed in Section 2, the running from the higher scale mixes the operators D7 and D5, leading to SI limits on D7. The same is true for the D8 operator, which in the running to low energy mixes into D6. Notwithstanding, the strongest limits on D8 come from the PandaX-II SD searches [12], since the D6 SI cross section is momentum suppressed in contrast with the D8 SD one. We can see from the Fig. 3 that the CMB data put stronger constraints on the D5 and D7 operators since they lead to unsuppressed s-wave annihilation cross sections. The limits on the operators D8 are much weaker since its s-wave cross section is dumped by a factor m 2 q /m 2 DM , while the operator D6 is not limited due to its p-wave annihilation cross section. Finally, we show in Figure 4 the non-collider limits on the operators D9 and D10 which exhibit tensor currents. The DD cross section due to operator D9 is spin dependent, however, it is not suppressed. Notice that the conservative limit on D9 is restricted to a very narrow DM mass range. To understand this result, it is interesting to compare the operators D8 and D9, since they have the same low energy limit. Despite the numerical coefficient of D9 being twice the one of D8 [45], the constraints on D8 are stronger than the ones on D9 because the bounds on the latter are weakened by the rescaling in the relic density given in Eq. (14). On the other hand, the SI bounds on D10 are the most stringent ones up to DM masses of the order of 400 GeV despite momentum suppressed cross section. The CMB bound on both operators D9 and D10 is important in part of the parameter space, since these operators exhibit s-wave annihilation cross sections. Moreover, for both operators, the relic density leaves only a small window below Λ ∼ 10 TeV unconstrained.
Before concluding, let us observe that also in the fermion DM case the CMB bound has positive slope, as was happening in the case of the C5 and C6 operators. Again, this is due to the fact that when the annihilation cross section is not velocity suppressed, it grows as m 2 DM .

Vector DM
Before we present the non-collider bounds on vector DM we would like to stress that the operators V1-V12 exhibit high energy asymptotic behaviors that correspond to dimension-seven and -eight operators [22], and therefore we included additional powers of m DM /Λ in their parametrisation; see Table 1. With this parametrisation the limits on Λ from the LHC and non-collider experiments are of the same order as the one for scalar and fermion DM operators; for further details see Ref. [22]. We present the non-collider limits on the effective operators containing vector DM in Fig. 5. When the exclusion plots are very similar, e.g. as is the case for the V2, V6, V9M and V10M operators, for instance, we show only one representative plot. Let us discuss these results starting with DD bounds. The most stringent SI constraints are on the V1 and V3 operators, since these operators lead to non-suppressed SI cross sections. We can see from panels 5(a) and 5(c) that for g = 1, the constraint on Λ can reach 10 TeV for large DM masses around 1 TeV. Another operator for which the SI constraints are relevant is V4, although it exhibits a pseudo-vector quark current. Notwithstanding, it develops a small mixing with the V3 operator due to the running from the scale Λ to the 1 GeV scale. As can be seen from Fig. 5(d), the SI DD limits on V4 are weaker than those on V1 or V3; still, they dominate the constraints for DM masses in the range 20-300 GeV. Also the V11 operator is quite constrained by SI searches, which are the dominant bound for m DM 500 GeV (notice that, although the bounds on V11 and V12 are presented in the same panel, the SI bound applies only to the V11 operator, since V12 is parity-violating and do not contribute to the SI cross section). The other twelve operators involving vector DM either do not contribute to SI processes or their contributions are velocity/momentum suppressed, and consequently, are not bound by SI DM DD searches. It is interesting to observe that the astrophysical uncertainties amount to a factor of 2.
The SD DM DD data can be used to constrain a few operators that lead to SD cross sections that are unsuppressed. We can see from Fig. 5(g) that for the V9P and V10P operators, the most constraining limits in the DM mass range 20-200 GeV stem from SD DM DD searches. On the other hand, the SD bounds on the operator V5, see Fig. 5(e), are rather weak. Nevertheless, it can be the most stringent one for DM masses between 10 and 50 GeV depending on the local DM density.
The CMB data constrain operators exhibiting s-wave DM annihilation channels that, in the case of vector DM, are V1, V2, V5, V6, V9M, V10M, V11, and V12. In the case of the operator V1, the DM bounds are looser than the ones coming from SI DM DD, as we can see from Fig. 5(a). On the other hand, for the operator V5, the CMB limits are tighter than the conservative SD DM DD ones; see Fig. 5(e). For the remaining operators constrained by CMB data, i.e. V2, V6, V9M, V10M and V12 the CMB bounds are the more robust ones on these operators as we can learn from Figs. 5(b) and 5(h). Nevertheless, the analyses of the relic density predicted by the operators in the last class indicate that vector operators are either ruled out (V1, V3, V4, V9P, and V10P) or very strongly constrained (V2, V5, V6, V7M, V7P, V8M, V8P, V9M, V10M, V11, and V12). As remarked before, these last limits can be evaded with simple modifications of the model.
It is interesting to notice that the operators V7P, V7M, V8P, and V8M are only bound by the relic density data and EFT validity region, as we can see from Fig. 5 PandaX-II PandaX-II

LHC sensitivity and its complementarity to non-collider constraints
In this section we present the LHC bounds on all the operators listed in Table 1 in the (m DM , Λ) plane, analogously to what has been done for the non-collider constraints in the previous section. Applying the analysis described in Section 3.3 we obtain limits on the DM EFT operators using the CMS monojet and mono W/Z (hadronically decaying) searches [5] at LHC Run 2 with an integrated luminosity of 35.9 fb −1 . Moreover, we also assess the LHC potential to probe the DM EFT for a projected integrated luminosity of 300 fb −1 .
Our results are shown in Figs. 6-10. In these figures, the red shaded area represents the LHC exclusion region at 95% CL for g = 1. Furthermore, the area inside the solid orange (blue) curve is excluded at 95% CL by the presently available CMS mono-jet data for g = 6 (g = 4π). For the sake of comparison, we also display in these figures the excluded region for g = 1 due to non-collider searches, represented by the light purple shaded area. We do not include the relic density constraints into the "non-collider" excluded area since these bounds are model-dependent and can easily be evaded, e.g. by the addition of co-annihilating DM partners. The region where the EFT approach is not valid (Λ > 2m DM ) is represented by the grey triangle in the right bottom side of the figures. Finally, the regions inside the dashed red, orange and blue curves represent our estimates for the future 95%CL exclusion by the LHC with 300 fb −1 , taking g = 1, 6 or 4π, respectively.
A common feature of the monojet excluded regions is that it exhibits upper and lower limits in Λ for a fixed DM mass. The lower limit is caused by the cut given in Eq. (16), that guarantees that the EFT is applied only within its validity limit. Furthermore, this cut has an important effect when the bound on Λ is low. On the other hand, for higher values of Λ (above few TeV) the impact of this cut on the cross sections is negligible for scalar operators. We will see below that there are also important differences in the effects of this cut for different DM spin cases since the DM pair invariant mass distributions depend strongly on the DM spin, as discussed in Ref. [22].

Complex scalar DM
We present in Fig. 6 the CMS monojet bounds on the EFT operators C1-C6 that contain scalar DM states. From the top left panel of this figure, we can learn that the LHC data is able to exclude up to DM masses of 70 GeV for Λ 600 GeV and g = 1. We can also see how these bounds evolve with the coupling g : as g varies, the upper limit on Λ does not scale as g 2 as naively expected. The reason for this behavior is the cut in Eq. (16), that reduces more significantly the available phase space for smaller g since this is associated to smaller Λ upper limits. This behavior is true for all operators containing vector and pseudo-vector quark currents. Moreover, this panel shows that the conservative non-collider bounds are the strongest ones for g = 1 This situation persists even when a larger integrated luminosity is accumulated.
The top right panel of Fig. 6 shows the collider bounds on the operator C2. As we can see, the collider limits on C1 and C2 are the same. Moreover, the most stringent limits on this operators originates from non-collider data, irrespective of the integrated luminosity. In addition, from the middle left panel of this figure we see that results for the operator C3 are similar to the ones for C1 and C2, except that the collider limit on Λ is reduced.
The scale of Λ which the LHC can probe for scalar DM EFT operators strongly depends on the effective operator: it is about 0.6 TeV for the C1 and C2 operators, it is only about 0.3 TeV for the C3 and C4 operators, and it is about 3 TeV for the C5 and C6 operators. The LHC searches are the least sensitive for the operators C3 and C4 for two reasons: first of all, these operators contain explicit momentum dependence through the derivative, and LHC cuts on the monojet are not hard enough to enhance this operator. Secondly the invariant DM pair mass distribution is shifted to higher values than for other scalar DM operators (see detailed discussion in Ref. [22]), therefore the cut given by the Eq. (16) reduces their signal more then for other EFT operators with scalar DM.
From Fig. 6 we can also see that LHC plays an important complementary role in probing the DM parameter space for the C4 and C6 operators. Notice that the LHC searches are especially important to test the operator C6 that involve gluons. In fact, the sensitivity of non-collider experiments to the EFT parameter space for this operator is very poor, as one can see from Figure 6(f). Moreover, even for the operators C1, C2 and C3, the LHC could help elucidate the nature of DM by independently probing the DM parameter space, especially if some DM signals would take place at collider and non-collider experiments.
In the case that the DM interactions are stronger, e.g. g = 4π, we can scale the CMB and DD limits shown in Fig. 6 by g 2 (g ) for dimension 5 (6) operators to compare with the LHC results. The conclusions for g = 4π are the same as the ones above for g = 1, except for the range of masses where the DM DD is dominant that shrinks to the interval 20-140 GeV. Moreover, we project that the LHC limits on C4 will be the strongest ones for an integrated luminosity of 300 fb −1 .
non-coll. (e) (f) non-coll.   Figure 6: LHC monojet constraints on EFT operators with scalar DM, as indicated in each panel. The area inside the red, orange and blue solid curves is excluded by current LHC data at 95% CL for g = 1, 6 and 4π, respectively. The projected LHC limits for 300 fb −1 are indicated by dashed thin lines. The combined exclusion regions from CMB and DM DD searches for g = 1 are given by the light-purple area below the purple curve. The grey triangular region marks the region where the EFT approach is not valid.

Dirac fermion DM
From Fig. 7 we can see that the LHC monojet data exclude DM masses up to 150-200 GeV and Λ 900 GeV for the operators containing scalar and pseudo-scalar quark currents and g = 1. Furthermore, in the case of effective operators containing vector and pseudo-vector quark currents, the LHC excluded region is slightly larger for g = 1: m DM 200 GeV and Λ 1.1 TeV. In the case of couplings approaching the strongly interacting regime (g = 4π) the exclusion region is extended to m DM 2 TeV and Λ 20 TeV for all operators D1-D8. Fig. 7 also allows us to see the complementarity between the collider and non-collider searches. First of all, the DM SI DD bounds are more stringent than the LHC ones for operators that possess unsuppressed contributions to DM SI DD, i.e. D1 and D5; see panels (a) and (e) of this figure. On the other hand, the LHC bounds are stronger than the non-collider ones for the operators D6 and D8 since these operators only exhibit velocity suppressed contributions to DM SI DD searches. The same conclusion applies to the operator D3 that is not limited by neither the CMB data nor the DD searches. Moreover, the CMB data and the LHC monojet searches complement nicely each other for the operators D2, D4 and D7 because the LHC searches dominates the bounds for DM masses smaller than 50 GeV while the CMB limits are stronger above this mass. Fig. 8 depicts the LHC limits on the operators D9 and D10 that contain a tensor quark current. The LHC monojet searches excludes the region m DM 600 GeV and Λ 3 TeV for these operators. For larger couplings g = 4π the region is expanded to m DM 2 TeV and Λ 35 TeV. For the operator D9 and g = 1, the CMB bounds dominates in a small region for m DM 600 GeV while the LHC limits are stronger below this mass. The most stringent bounds on the operator D10 have three different sources depending on the DM mass for g = 1. At low DM masses 15 GeV the bounds are dominated by the LHC searches while for heavier masses 300 GeV the most important constraints come from the CMB data. On the other hand, the SI limits are more stringent in the mass window 15 m DM 300 GeV.
As for the case of scalar DM, DD and CMB limits on the Dirac fermion DM rescale as g 2 (g ) for dimension 5 (6) operators. In the case of g = 4π the region where the CMB constraints are more stringent than the collider ones on D2, D4, D7, and D10 changes to DM masses in excess of 350, 350, 400, and 620 GeV respectively. In the case of the other effective operators for Dirac fermion DM the conclusions remain the approximately the same for the value of g .

Vector DM
We display in Figures 9 and 10 the present constraints on vector DM operators that stem from the LHC monojet searches. Notice that the LHC cross sections do not change significantly when we replace a current by the corresponding pseudo one. Consequently, we do not present the results for the pseudo currents unless their non-collider constraints are different.
From the top panels of Fig. 9 we learn that the LHC monojet constraints on the operators V1 and V2 are equal: for g = 1, the excluded region is given by m DM 40 GeV and 100 ≤ Λ ≤ 400 GeV. As g increases to 4π, the excluded region becomes m DM 2 TeV and 100 ≤ Λ ≤ 9000 GeV. Notice that the lower limit on Λ (a) non-coll.
non-coll.  does not depend on the coupling g . Since the V1 contribution to the SI DM DD is unsuppressed the strongest bound on this operator comes from the DD searches. On the other hand, for the operator V2, the collider and non-collider limits are clearly complementary: the collider limits dominate for DM masses smaller than 40 GeV while the CMB data give rise to the most stringent limits at DM masses larger than 40 GeV. The operators V3 and V4 do not lead to any collider limit for g = 1 even for a larger integrated luminosity (300 fb −1 ), therefore, the non-collider limits are the most important ones, showing again the synergy between low and high energy data. This lack of sensitivity on these operators for g = 1 originates from the shape of the invariant mass distribution of the DM pairs that is shifted towards high values. Consequently, the EFT validity cut in Eq. (16) discards a large fraction of the events in this case. This reduction of number of signal events is overcome only for larger values of coupling g = 6 (4π) for which the LHC reach can be m DM 600 (900) GeV and Λ 2.5 (4) TeV.
The top left panel of Fig. 10 displays the LHC monojet constraints on the operators V5 and V6, which are identical. For g = 1 the LHC bounds are the most stringent on these operators for DM masses up to 180 GeV and λ 600 GeV while the CMB limits play an important role just in a small DM mass range around 200 GeV. On the other hand, we can see from the right top panel of this figure, that the operators V7M and V8M are not bound at all for g = 1. In the collider case the reason is the same behind the null results for the operators V3 and V4. These operators can only be probed by the LHC monojet searches for higher values of g .
The left middle panel of Fig. 10 contains the bounds on the operators V7P and V8P that are only constrained by the LHC monojet searches that exclude the narrow region m DM 100 GeV and 250 Λ 450 GeV for g = 1. For higher values of the coupling, e.g. g = 4π this region is substantially expanded to m DM 1.3 TeV and 250 Λ 6500 GeV. Once again the lower limit on Λ depends on g and on the integrated luminosity. The limits on the operators V9M and V10M are equal and are shown in the right middle panel of Fig. 10. In this case we witness a nice complementarity between high and low energy data for these operators for g = 1.
non-coll. We can see from the left lower panel of Fig. 10 that the limits on the operators V9P and V10P are rather loose, however, the non-collider and collider ones are complementary. The LHC monojet data exclude the area m DM 40 GeV and 160 Λ 280 GeV, while the SD DM DD searches are more relevant for DM masses in the range 40-100 GeV. Moreover, the excluded area is increase for larger values of g or for larger integrated luminosities. Last, but not least, the limits on the operators V11 and V12 are presented in the lower right panel of this figure. In this case, the monojet constraints are more stringent than the direct detection (valid only for the operator V11) and CMB ones, excluding a large fraction of the parameter space for g = 1: m DM 400 GeV and 200 Λ 1200 GeV. Notice that these are the only operators that the Λ lower limits changes for higher couplings or integrated luminosities.
In brief, from Figs 9 and 10 we learn that LHC is sensitive to a substantial numbers of operators containing vector DM, such as V2, V5 V6, V7P, V8P, V9-V12 which can not be probed at all or are poorly bound with non-collider searches. The LHC limits strongly depend on the type of the operator ranging from few hundred GeV (V1, V2, V5,V6, V7P, V8P, V9M, V9P, V10M, V10P) to about 1.3 TeV (V11,V12). Moreover, we also assess the impact of the EFT validity cut in Eq. (16) for vector DM, a fact that was not explored in previous studies.
In this scenario, the DD and CMB constraints rescale at a fixed DM mass as g 2/3 and g 1/2 for operators suppressed by Λ 3 and Λ 4 respectively. For g = 4π, the dominance of the CMB bounds on the operators V2, V5/V6, V9M/V10M, and V11/V12 occurs for DM masses larger than 920, 1000, 850, and 1500 GeV respectively. Moreover, the collider limits on V3 are more stringent than the DD ones for m DM ≤ 50 GeV. For the other vector DM operators, the conclusions remain the same for the upper bounds on Λ. Nevertheless, there is a beautiful complementarity between the non-collider and LHC searches since the former exclude the region of small Λ that the collider searches do not cover due to the validity of the EFT.

Conclusions
In this work we have presented accurate and up-to-date constraints on the complete set of dimension five and six operators connecting SM quarks and gluons with a DM candidate, which can be a complex scalar, a Dirac fermion or a complex vector.
We have performed a comprehensive analyses of the complementarity between collider and non-collider searches to probe DM parameter space, including LHC mono-jet data, bounds from SI and SD direct searches, relic density limits and CMB indirect constraints due to the injection of energy produced by DM annihilation in the early universe. Since the characteristic energy scale for LHC and direct searches differs by about six orders of magnitude, to correctly evaluate the experimental sensitivity to the DM EFT operators we have taken into account their running and mixing from the TeV scale to the GeV one. This effect is especially important for operators with pseudo-vector SM quark current (C4, D6, D7 and V4 in Table 1) which, mixing into operators with vector SM quark current, develop a SI cross section. Another important point has been to take into account the realistic uncertainty in the DM DD searches limits due to the uncertainty on the local DM density. In some cases, the uncertainty can quantitatively shift the bounds by one order of magnitude, an important result that must be taken into account to properly address the excluded parameter space. As for indirect searches, we have chosen to present only the bounds coming from CMB data, since they are not plagued by uncertainties on the DM distribution. These bounds are particularly important when the direct detection rates are suppressed. We have also found that the EFT validity criteria M DM,DM < Λ plays an important role for the case of vector DM, when the invariant mass of the DM pair M DM,DM is typically larger than in the case of scalar and fermion DM, and there is no collider bound for relatively small values of Λ.
Our updated bounds are summarized in Figs. 6-10 for all the operators listed in Table 1. In the case of scalar DM we have found an important synergy between the different experiments: while for the operators C1 and C3 (with scalar and vector quark currents, respectively) the most stringent bounds stem from SI searches, the best limits on the C5 and C6 operators (describing DM-DM-Gluon-Gluon interactions) come from the LHC mono-jet data, since DD rates of DM scattering on gluon component of the nucleons are suppressed. At the same time, the most important constraints on the C2 operator come from the CMB data, while the limits on C4 are dominated by the LHC searches at low m DM and by SI bounds for large DM mass values.
The overall picture for the Dirac fermion DM is similar to the one of scalar DM: there is a synergy between the collider and non-collider searches. The strongest limits on the D1 and D5 operators are due to the SI searches, while the LHC mono-jet data put the strongest constraints on the D3, D6 and D8 operators. On the other hand, the operators D2, D4 and D9 are bounded by the LHC mono-jet searches at small DM masses, while at larger masses the CMB data dominate the limits. At the same time, the best limits on the D7 and D10 operators are defined by the interplay of the LHC searches, direct detection data and CMB constraints, depending on the DM mass.
In the case of vector DM we can also see an important complementarity between the different experiments to probe DM parameter space. The most stringent limits on the V1, V3 and V4 Wilson coefficients are set by SI searches. On the other hand, the LHC mono-jet data provides the strongest bounds on the V5, V6, V7P, V8P, V11, and V12 operators for any value of the DM mass. Furthermore, there is a synergy between the LHC and CMB data in probing the V2, V9M and V10M operators in different DM mass regions. One should also note the interplay between the LHC mono-jet searches and the SD bounds in bounding the V9P and V10P operators. We have also found that notwithstanding of combination of all data there are no limits 22 on the operators V7M and V8M for g = 1. There is one more important point to stress about the interplay of different data in probing vector DM operators: the lack of the LHC sensitivity to small values of Λ discussed above in section 4.2 is complemented by the potential of CMB data to probe this region. This complementarity is manifest for many operators (V1, V2, V3, V4, V9M, V10M, V9P, V10P, V11 and V12), for which CMB data partly or completely cover the lower Λ region which LHC is unable to probe.
As a general remark, in our analysis we have assumed that just one operator is non-vanishing at a time, which may or may not be the case depending on the underlying theory realized in nature. Nevertheless, our studies can be easily adapted to some specific scenario where the integration of heavy mediators can lead to more than one non-vanishing Wilson coefficient. For instance, a t-channel scalar mediator coupled to fermion DM would generate a (χq)(qχ) effective operator. After a Fierz transformation, this operator can be written as a combination of the fermion operators D1, D4, D5, D8, and D9 [22]. According to Fig. 7, it is clear that in this case the most stringent bounds would be set by SI searches, so that at least in first approximation the bound on (χq)(qχ) coincides with the bound on D1. Of course, the validity of such a procedure must be analyzed case by case.