Thermal WIMPs and the Scale of New Physics: Global Fits of Dirac Dark Matter Effective Field Theories

We assess the status of a wide class of WIMP dark matter (DM) models in light of the latest experimental results using the global fitting framework $\textsf{GAMBIT}$. We perform a global analysis of effective field theory (EFT) operators describing the interactions between a gauge-singlet Dirac fermion and the Standard Model quarks, the gluons and the photon. In this bottom-up approach, we simultaneously vary the coefficients of 14 such operators up to dimension 7, along with the DM mass, the scale of new physics and several nuisance parameters. Our likelihood functions include the latest data from $\mathit{Planck}$, direct and indirect detection experiments, and the LHC. For DM masses below 100 GeV, we find that it is impossible to satisfy all constraints simultaneously while maintaining EFT validity at LHC energies. For new physics scales around 1 TeV, our results are influenced by several small excesses in the LHC data and depend on the prescription that we adopt to ensure EFT validity. Furthermore, we find large regions of viable parameter space where the EFT is valid and the relic density can be reproduced, implying that WIMPs can still account for the DM of the universe while being consistent with the latest data.


Introduction
Despite years of searching, the identity of dark matter (DM) remains a mystery. Nevertheless, the large number of past, present and future probes of its particle interactions makes it essential to regularly revisit the constraints on the most popular theoretical candidates, in order to guide future searches.
A favoured paradigm for the particle nature of dark matter is that of Weakly Interacting Massive Particles (WIMPs), due to the fact that it allows for a simple thermal mechanism to produce DM with the cosmologicallyobserved abundance [1]. Such models have also attracted attention due to the large number of possible signals they predict, none of which have been definitively observed so far. Although this has led some to make claims of the demise of WIMPs [2], others have argued that such predictions are premature [3].
A relatively agnostic approach to WIMP model building is to pursue a bottom-up, Effective Field Theory (EFT) approach, in which one enumerates all of the allowed higher dimensional operators which lead to interactions between DM and Standard Model (SM) particles. Any result described by an EFT can in general be explained by many high-energy theories. In this way, the EFT description is a model-independent one, as it does not depend on the Ultraviolet (UV) completion that describes an effective operator. This is, however, a double-edged sword: because an effective operator does not encode any information about the UV completion, it has no constraining power in distinguishing between the range of UV theories that can map to it -nor can all UV-complete theories be mapped to an EFT description for the energies we are interested in here.
A common approach to the analysis of EFTs for DM in the literature has been to consider a single operator at a time [35][36][37][38][39][40][41] and compare experimental bounds on the new physics scale Λ with the values implied by the observed DM relic density. This method, however, severely limits the scope of the analysis and potentially leads to overly-aggressive exclusions, not only because it neglects (potentially destructive) interferences between different operators [42], but also because the relic density constraint can be considerably relaxed when several operators contribute to the DM annihilation cross-section. The first global study of EFTs for scalar, fermionic and vector DM taking interference effects into account was performed in Ref. [43], but no collider constraints were included in the analysis and no couplings to gluons were considered. More recently, Ref. [44] applied Bayesian methods to perform a global analysis of scalar DM, for which only a small number of effective operators need to be considered and collider constraints can be neglected. Examples of global studies considering subspaces of a general DM EFT include Refs. [45][46][47].
In the present work, we exploit the computational power of the GAMBIT framework [48] to perform the first global analysis of a very general set of effective operators up to dimension 7 that describe the interactions between a Dirac fermion DM particle (or a DM sub-component) and quarks or gluons. Such a set-up arises for example in many extensions of the SM gauge group, such as gauged baryon number [49] or other anomaly-free gauge extensions that require additional stable fermions [50,51]. Our novel approach of considering many operators simultaneously enables us to study parameter regions where several types of DM interactions need to be combined in order to satisfy all constraints. Our analysis substantially improves upon the previous state-of-the-art in both the statistical rigour with which the DM EFT parameter space is interrogated, and in the new combinations of constraints that are simultaneously applied. We also increase the level of detail with which individual constraints are modelled, summarised as follows.
First, we include a much improved calculation of direct detection constraints using the GAMBIT module DarkBit [52]. We consider the renormalization group (RG) evolution of all effective operators from the electroweak to the hadronic scale and then match the relativistic operators onto the non-relativistic effective theory [53] relevant for DM-nucleon scattering. We then calculate event rates in direct detection experiments to leading order in the chiral expansion, including the contributions from operators that are naively suppressed in the non-relativistic limit, and determine the resulting constraints using detailed likelihood functions for a large number of recent experiments. In the process, we include a number of nuisance parameters to account for uncertainties in nuclear form factors and the astrophysical distribution of DM.
Second, we consider the most recent constraints on DM annihilations using gamma rays and the Cosmic Microwave Background (CMB). To include the latter, we employ the recently released GAMBIT module Cos-moBit [54], which uses detailed spectra to calculate effective functions for the efficiency of the injected energy deposition and obtain constraints on the DM annihilation cross-section while varying cosmological parameters. For the calculation of annihilation cross-sections we make use of the new GAMBIT Universal Model Machine (GUM) [55,56] to automatically generate the relevant code based on the EFT Lagrangian.
Third, we combine the above detailed astrophysical and cosmological constraints with a state-of-the-art implementation of LHC constraints on WIMP dark matter. A central concern for any study of EFTs is the range of validity of the EFT approach [57][58][59][60][61][62][63][64]. This is particularly true when considering constraints from the LHC, which may probe energies above the assumed scale of new physics. A naive application of the EFT in such a case may lead to unphysical predictions, such as unitarity violation. Whenever this is the case it becomes essential to adopt some form of truncation to ensure that only reliable predictions are used to calculate experimental constraints.
In the present work we address these challenges in two key ways. First, we separate the scale of new physics Λ from the individual Wilson coefficients C (rather than scanning over a combination such as C/Λ 2 ), such that the former can be directly interpreted as the scale where the EFT breaks down and the latter can be constrained by perturbativity. Second, we check the impact of a phenomenological nuisance parameter that describes the possible modification of LHC spectra at energies beyond the range of EFT validity. The nuisance parameter smoothly interpolates between an abrupt truncation and no truncation at all.
Our analysis reveals viable parameter regions for general WIMP models across a wide range of new physics scales, including very small values of Λ (Λ < 200 GeV), where there are no relevant LHC constraints and very large values of Λ (Λ > 1.5 TeV), where LHC constraints are largely robust. Of particular interest are the intermediate values of Λ (Λ ∼ 700-900 GeV), for which our DM EFT partly accommodates several small LHC data excesses that could be interesting to analyse in more detail in the context of specific UV completions or simplified models. However, our analysis also reveals that there cannot be a large hierarchy between Λ and the DM mass m χ . In particular, even with the most general set of operators we consider, it is impossible to simultaneously have a small DM mass (m χ 100 GeV) and a large new physics scale (Λ > 200 GeV). In other words, for light DM to be consistent with all constraints, it is necessary for the new physics scale to be so low that the EFT approach breaks down for the calculation of LHC constraints. For heavier DM, on the other hand, thermal production of DM in the early universe would exceed the observed abundance whenever Λ is more than one order of magnitude larger than m χ (up to the unitarity bound at a few hundred TeV [65], where the maximum possible value of Λ approaches m χ ).
This work is organised as follows. We introduce the DM EFT description in Sec. 2. In Sec. 3, we discuss the constraints used in this study, and our methods for computing likelihoods and observables. We present our results in Sec. 4. Finally, we present our conclusions in Sec. 5. The samples from our scans and the corresponding GAMBIT input files, and plotting scripts can be downloaded from Zenodo [66].

Dark Matter Effective Field Theory
In this study, we consider possible interactions of SM fields with a Dirac fermion DM field, χ, that is a singlet under the SM gauge group. For phenomenological reasons discussed in detail in Sec. 3, we focus on interactions between χ and the quarks or gluons of the SM. We assume that the mediators that generate these interactions are heavier than the scales probed by the experiments under consideration. Following the notation of Refs. [67,68], the interaction Lagrangian for the theory can be written as where Q such that the free parameters of the theory are the DM mass m χ , the scale of new physics Λ, and the set of dimensionless Wilson coefficients {C (d) a }. For sufficiently large Λ, the phenomenology at small energies is dominated by the operators of lowest dimension, and we therefore limit ourselves to d ≤ 7. However, even this leaves a relatively large set of operators. The DM EFT that is valid below the electroweak (EW) scale (with the Higgs, W , Z and the top quark integrated out) contains 2 dimension five, 4 dimension six, and 22 dimension seven operators (not counting flavour multiplicities), while the DM EFT above the EW scale for a singlet Dirac fermion DM has 4 dimension five, 12 dimension six, and 41 dimension seven operators (again, not counting flavour multiplicities) [68]. The large set of possible operators poses a challenge for a global statistical analysis where bounds on Λ and {C (d) a } are derived from experimental observations (see Sec. 3 for details). An added complexity is that we consider both processes where the typical energy transfer is above the EW scale (such as collider searches and indirect detection) as well as processes in which the energy release is small (direct detection). The consistent implementation of these bounds requires the combination of both DM EFTs, together with the appropriate matching conditions between the two.
To make the problem tractable we focus in our numerical analysis on a subset of DM EFT operators -the dimension six operators involving DM, χ, and SM quark fields, q, The difference between the DM EFT below the EW scale and the DM EFT above the EW scale is in this case very simple: above the EW scale the quark flavours run over all SM quarks, including the top quark, while below the EW scale the top quark is absent. While the above set of operators does not span the full dimension six bases of the two DM EFTs, it does collect the most relevant operators. The full dimension six operator basis contains operators where quarks are replaced by the SM leptons. These are irrelevant for the collider and direct detection constraints we consider, and are thus omitted for simplicity. The basis of dimension six operators for the DM EFT above the EW scale contains, in addition, operators that are products of DM and Higgs currents. These are expected to be tightly constrained by direct detection to have very small coefficients such that they are irrelevant in other observables, and are thus also dropped for simplicity.
To explore to what extent the numerical analyses would change, if the set of considered DM EFT operators were enlarged, we also perform global fits including, in addition to the dimension six operators (3)-(6), a set of dimension seven operators that comprise interactions with the gluon field either through the QCD field strength tensor G a µν or its dual G µν = 1 2 µνρσ G ρσ , as well as operators constructed from scalar, pseudoscalar and tensor bilinears: 10,q = m q (χiσ µν γ 5 χ)(qσ µν q) .
The definition of the operators describing interactions with the gluons, Q 1-4 , includes a loop factor since in most new physics models these operators are generated at one loop. Similarly, the couplings to scalar and tensor quark bilinears, Q 5-10,q , include a conventional factor of the quark mass m q , since they have the same flavour structure as the quark mass terms (coupling left-handed and right-handed quark fields). The m q suppression of these operators is thus naturally encountered in new physics models that satisfy low energy flavour constraints, such as minimal flavour violation and its extensions. Note that, unless explicitly stated otherwise, m q always refers to the running mass in the modified minimal subtraction (MS) scheme.
The complete dimension-seven basis below the EW scale contains eight additional operators with derivatives acting on the DM fields [68]. To simplify the discussion we do not include these operators in our analysis, partially because they do not lead to new chiral structures in the SM currents. Moreover, the direct detection constraints on these additional operators are expressible in terms of the operators that we do include in the global fits due to the non-relativistic nature of the scattering process.
Note that the operators Q (7) 5-10,q are not invariant under EW gauge transformations, and are thus replaced in the DM EFT above the EW scale by operators of the form (χχ)(q L q R )H, where H is the Higgs doublet. In all the processes we consider, H can be replaced by its vacuum expectation value -either because the emission of the Higgs boson is phase-space suppressed or suppressed by small Yukawa couplings, or both. This means that, up to renormalization group effects (to be discussed in Sec. 2.1), the operators Q 5-10,q can also be used in our fitting procedure above the EW scale.
In principle, analogous operators to Q 1-10,q exist for leptons instead of quarks [69,70] and weak gauge bosons instead of gluons [71][72][73]. 1 In general, these play a much smaller role in the phenomenology and will not be considered here. Similarly, throughout this work the Wilson coefficients of any dimension five operators are set to zero at the UV scale.
The Wilson coefficients of the operators defined above depend implicitly on the energy scale of the process under consideration. In our fits, all Wilson coefficients are specified at the new physics scale Λ. If this scale is larger than the top mass, Λ > m t , all six quarks are active degrees of freedom and the Wilson coefficients need to be specified for q = u, d, s, c, b, t. For Λ < m t , the top quarks are integrated out, and only the Wilson coefficients for q = u, d, s, c, b need to be specified. This is done automatically in our fitting procedures, such that effectively both EFTs are used in the fit, according to the numerical value of the scale Λ.
Although, a priori, the Wilson coefficients for each quark flavour are independent, we will restrict ourselves to the assumption of minimal flavour violation (which implies C i,t ), and the assumption of isospin invariance (which implies C i,u ). 2 Hence, each operator comes with only one free parameter in addition to the global parameters Λ and m χ . Under these assumptions, the two EFTs above and below the EW scale have the same number of free parameters.
1 Furthermore, there are two additional dimension-6 operators describing DM-photon interactions: the anapole moment and the charge radius [74]. For a recent discussion of LHC constraints on these operators we refer to Ref. [75]. 2 These constraints also ensure that the dimension-six operators do not explicitly break electroweak symmetry, which requires C 1,u − C (6) 3,u = C

Running and mixing
For many applications, the RG running of the Wilson coefficients (i.e. their dependence on the energy scale µ) can be neglected. In fact, the operators Q 2,q , Q (7) 5,q and Q (7) 6,q have vanishing anomalous dimension, while Q 3,q , Q 4,q , Q 7,q , Q 8,q as well as Q 1-4 exhibit no running at one-loop order in QCD [77]. Nevertheless, there are two cases when the effects of running can be important: 1. Mixing: Different operators can mix with each other under RG evolution, such that operators assumed negligible at one scale may give a relevant contribution at a different scale. This is particularly important in the context of direct detection, because for certain operators the DM-nucleon scattering crosssection is strongly suppressed in the non-relativistic limit. In such a case, the dominant contribution to direct detection may arise from operators induced only at the loop level [78,79]. In our case, the dominant effects arise from the top quark Yukawa and are discussed below. 2. Threshold corrections: Whenever the scale µ drops below the mass of one of the quarks, the number of active degrees of freedom is reduced and a finite correction to various operators arises. In our context, the only effect is the matching of the operators Q 5-8,q onto the operators Q 1-4 at the heavy quark thresholds, which is given by Mixing of the tensor operators Q 9,q , Q 10,q above the EW scale and subsequent matching gives rise to the dimension-five dipole operators where F µν is the electromagnetic field strength tensor and e is the electromagnetic charge. These operators give an important contribution to direct detection experiments and are thus kept. 3 3 Note that as per our assumptions the Wilson coefficients In the present work we include these effects as follows. To calculate the Wilson coefficients at the hadronic scale µ = 2 GeV (relevant for direct detection) we make use of the public code DirectDM v2.2.0 [67,68], which calculates the RG evolution of the operators defined above, including threshold corrections and mixing effects. The code furthermore performs a matching of the resulting operators at µ = 2 GeV onto the basis of non-relativistic effective operators relevant for DM direct detection (see Sec. 3.1).
DirectDM currently requires as input the Wilson coefficients in the five-flavour scheme given at the scale m Z = 91.1876 GeV. For Λ < m t (five-flavour EFT), we can therefore directly pass the Wilson coefficients defined above to DirectDM. For Λ > m t (six-flavour EFT), there are three additional effects that are considered. First, as pointed out in Ref. [80], the operators Q (7) 9,10,t give a contribution to the dipole operators Q (5) 1,2 at the one-loop level, which is given by 4 Second, as pointed out first in Ref. [81], the operator with an axial-vector top-quark current Q 3,t mixes into the operators Q after integrating out the Z boson at the weak scale. Here, s w ≡ sin θ w with θ w the weak mixing angle, and v = 246 GeV is the Higgs field vacuum expectation value. The flavour universal UV contributions C (6) 1,u/d (Λ) largely compensate the mixing effect in the fit; the remnant effect, due to the isospin-breaking Z couplings, is small.
Third, in order to match the EFT with six active quark flavours onto the five-flavour scheme, we need to integrate out the top quark and apply the top quark threshold corrections given in Eq. (17). We neglect any other effects of RG evolution between the scales Λ and m Z , i.e. all Wilson coefficients other than C are directly passed to DirectDM. 5 For the purpose of calculating the LHC constraints, we neglect the effects of running and do not consider loop-induced mixing between different operators, which 4 For historical reasons, in the numerical code log(m 2 t /Λ 2 ) instead of log(m 2 Z /Λ 2 ) was used. The effect on the numerical results is negligible. 5 Small remnant effects of the bottom and charm Yukawa coupling are taken into account below the EW scale via double weak insertions [79] that are included in the DirectDM code. Fig. 1: Illustration of our approach for studying DM EFTs compared to a more naive approach, in which one only uses the experiment that yields the strongest bound on C/Λ 2 . The resulting exclusion is indicated by the red shaded region. By independently varying Λ, we can include additional information from experiments that give weaker bounds on C/Λ 2 but for which the EFT has a larger range of validity. The additional exclusion obtained in this way is indicated by the blue shaded region. The region of parameter space that corresponds to the non-perturbative values of Wilson coefficient C is excluded in either approach (shaded brown).
is a good approximation for the operators Q  5-10,q mixing effects are known to be important in principle [82], but these operators are currently unconstrained by the LHC in the parameter region where the EFT is valid (see Sec. 2.2). Likewise we also calculate DM annihilation cross-sections at tree level. In particular, in these calculations we neglect the running of the strong coupling (α s ) and use the pole quark masses (m pole q ) instead of the running quark masses. Moreover, we neglect a small loop-level contribution from the operators Q

EFT validity
A central concern when employing an EFT to capture the effects of new physics is that the scale of new physics must be sufficiently large compared to the energy scales of interest for the EFT description to be valid. Unfortunately, the point at which the EFT breaks down is difficult to determine from the low-energy theory alone. Considerations of unitarity violation make it possible to determine the scale where the EFT becomes unphysical, but in many cases the EFT description already fails at lower energies, in particular if the UV completion is weakly coupled.
To address this issue in the present study, we simultaneously vary the overall scale Λ, which corresponds to the energy where new degrees of freedom become relevant and the EFT description breaks down, and the Wilson coefficients C (d) a for each operator. Doing so introduces a degeneracy, because cross sections are invariant under the rescaling Λ → αΛ and C (d) a . However, the advantage of this approach is that the parameter Λ can be used to determine which constraints can be trusted in the EFT limit. This is illustrated in Fig. 1, which compares our approach of varying Λ and C (d) a separately to the naive approach where only C We emphasize that this approach assumes the same new-physics scale for all effective operators, even though they may be generated through different mechanisms, and hence at different scales, in the UV. In practice, one should think of Λ as the minimum of all of these scales, i.e. the energy at which new degrees of freedom first become relevant. These new degrees of freedom may not contribute to all processes, such that some effective operators may provide an accurate description even at energies above Λ. Whether or not this is the case cannot be determined from the low-energy viewpoint, such that we conservatively limit the EFT validity to energies below Λ.
For the purpose of direct detection constraints, the only requirement on Λ is that it is larger than the hadronic scale, so that the effective operators can be written in terms of free quarks and gluons. This is the case for Λ 2 GeV, which will always be satisfied in the present study. However, in order to evaluate direct detection constraints, it is necessary to determine the relic abundance of DM particles, which depends on the cross-sections for the processes χχ → qq or χχ → gg, just as in the case of indirect detection constraints (see Sec. 3.3). For this calculation to be meaningful in the EFT framework, we require Λ > 2m χ . Parameter points with smaller values of Λ will thus be invalidated. A dedicated study of direct detection constraints for Λ < 2m χ will be left for future work.
In the context of LHC searches for DM, EFT validity requires that the invariant mass of the DM pair produced in a collision satisfies m χχ < Λ [83]. To obtain robust constraints, only events with smaller energy transfer should be included in the calculation of likelihoods. The problem with this prescription is that m χχ does not directly correspond to any observable quantity (such as the missing energy / E T of the event) and hence the impact of varying Λ on predicted LHC spectra is difficult to assess. One possible way to address this issue would be to generate new LHC events for each parameter point and include only those events with small enough m χχ in the likelihood calculation, but this is not computationally feasible in the context of a global scan.
In the present work, we adopt the following simpler approach: Rather than comparing Λ to the invariant mass of the DM pair, we compare it to the typical overall energy scale of the event, which can be estimated by the amount of missing energy produced. In other words, we do not modify the missing energy spectrum for / E T < Λ and only apply the EFT validity requirement for larger values of / E T . This approach is less conservative than the one advocated, for instance in Refs. [63,64], where the energy scale of the event is taken to be the partonic centre-of-mass energy √ŝ , but it has the crucial advantage that it can be applied after event generation, since the differential cross-section with respect to missing energy dσ/d / E T is exactly the quantity that is directly compared to data. 6 In the following, we will consider two different prescriptions for how to impose the EFT validity. The first one is to introduce a hard cut-off, i.e. to set dσ/d / The second, more realistic, prescription is to introduce a smooth cut-off that leads to a non-zero but steeply falling missing energy spectrum above Λ. For this we make the replacement for / E T > Λ. Here a is a free parameter that depends on the specific UV completion. The limits a → 0 and a → ∞ correspond to no truncation and an abrupt truncation above the cut-off, respectively. For the case that the EFT results from the exchange of an s-channel mediator with mass close to Λ, one finds a ≈ 2 [30]. Rather than taking inspiration from a specific UV completion, we will instead keep a as a free parameter in the interval [0, 4] and find the value that gives the best fit to data at each parameter point. This approach typically leads to conservative LHC bounds in the sense that much stronger exclusions may be obtained in specific UV completions, if the heavy particles that generate the effective DM interactions can be directly produced at the LHC. However, this truncation procedure can lead to unrealistic spectral shapes with sharp features that may be tuned to fit fluctuations in the data. As will be discussed in more detail in Sec. 4, any explanation of data excesses through this approach must be interpreted with care.
Without upper bounds on the Wilson coefficients, any requirement on EFT validity could be satisfied by 6 We emphasize that mχχ and / E T are not strongly correlated in the sense that there are events with both / E T mχχ (if the DM pair is emitted approximately in the longitudinal direction) and / E T mχχ (if the two DM particles are light and approximately collinear). Since our approach does not modify the spectrum for / E T < Λ, we risk overestimating the differential cross-section in this regime. However, the sensitivity of the LHC to DM EFTs typically stems from events with large / E T , where our prescription is more appropriate. making both Λ and the Wilson coefficients arbitrarily large. We therefore require C (d) a < 4π, which is necessary for a perturbative UV completion and ensures that there is no unitarity violation in the validity range of the EFT [62].
One drawback of this prescription is that the EFT validity requirement depends on the normalisation of the effective operators. For example, we have written Q (7) 1-2 with a prefactor α s /(12π) and Q 3-4 with a prefactor α s /(8π) to reflect the fact that in many UV completions, these operators would be generated at the one-loop level. If these operators are instead generated at tree level (e.g. from a strongly interacting theory), it would be more appropriate to write the prefactor as 4πα s . With the latter convention any constraint on the new physics scale Λ becomes stronger by a factor (48π 2 ) 1/3 ≈ 5.3 for Q (7) 1-2 and by a factor (32π 2 ) 1/3 ≈ 4.6 for Q 3-4 , meaning that much larger values of Λ are experimentally testable and the range of EFT validity is substantially increased. We have confirmed explicitly that the results presented in Sec. 4 do not depend on the specific definition of the Wilson coefficients for Q 1-4 . 7

Parameter ranges
In this study we focus on the following parameter regions. In order to be able to neglect QCD resonances in the process χχ → qq, we restrict ourselves to m χ > 5 GeV. In order to have a sufficiently large separation of scales between the new physics scale Λ and the hadronic scale, we also require Λ > 20 GeV. As discussed in Sec. 2.2, we furthermore impose the bound |C (d) a | < 4π on all Wilson coefficients and the bound Λ > 2m χ . The upper bounds on m χ and Λ depend on the details of the scans that we perform and will be discussed in Sec. 4.

Constraints
In this section we describe the constraints relevant for our model. A summary of all likelihoods included in our scans is provided in Table 1. For each likelihood that directly constrains the interactions of the DM particle we also quote the background-only log-likelihood ln L bg obtained when setting all Wilson coefficients to zero. For the remaining likelihoods we instead quote the maximum achievable value of the log-likelihood ln L max . The sum of all these contributions, ln L ideal = −105.3 will be used to calculate log-likelihood differences below. 7 We note that the explicit factor of mq in the definition of Q 5-10,q not only affects the EFT validity but also directly the resulting phenomenology. Hence our results cannot be easily translated to operators with non-trivial flavour structure. Nuisances (see Table 4) 5.141 Table 1: Likelihoods included in our scans and their respective values for the background-only hypothesis. For each likelihood that directly constrains the interactions of the DM particle we also quote the background-only log-likelihood ln L bg obtained when setting all Wilson coefficients to zero. For the remaining likelihoods we instead quote the maximum achievable value of the log-likelihood ln L max .

Direct detection
Direct detection experiments search for the scattering of DM particles from the Galactic halo off nuclei in an ultra-pure target by measuring the energy E R of recoiling nuclei. The differential event rate with respect to recoil energy is given by where ρ 0 is the local DM density, m T is the target nucleus mass, f (v) is the local DM velocity distribution and is the minimal DM velocity to cause a recoil carrying away a kinetic energy is the reduced mass of the DM-nucleus system. The local DM density and velocity distribution are not very well known and introduce sizeable uncertainties in the prediction of experimental signals (see the discussion of nuisance parameters in Sec. 3.6). Nevertheless, the greatest challenge in the present context is the calculation of the differential scattering cross-section dσ/dE R . For this purpose, one needs to map the effective interactions between relativistic DM particles and quarks or gluons defined above onto effective interactions between non-relativistic DM particles and nucleons N = p, n. Table 2: A full list of dimension-6 and 7 operators included in this study, and the types of interactions they induce. For the DM-nucleon scattering cross-section, we distinguish between spin-independent (SI) and spin-dependent (SD) interactions, with the former receiving a large coherent enhancement and the latter vanishing for nuclei with zero spin. We use "unsuppressed" ("suppressed") to denote tree-level contributions that do not vanish (that vanish) in the zero-velocity limit, while "loop-induced" implies that an unsuppressed interaction is induced at the one-loop level. For the annihilation cross-section we use "s-wave" ("p-wave") to denote annihilations that do not vanish (that vanish) in the zero-velocity limit. Note that if the s-wave contribution is helicity suppressed (i.e. proportional to m 2 q /m 2 χ ), the p-wave contribution may dominate in the relic density calculation.

Dimension-6 operators
The EFT of non-relativistic interactions can be written as where the operators O N i depend only on the DM spin S χ , the nucleon spin S N , the momentum transfer q and the DM-nucleon relative velocity v [4,53,99].
The non-relativistic operators can be divided into four categories according to whether or not they depend on the nucleon spin S N , such that scattering is suppressed for nuclei with vanishing spin, and whether or not they depend on q and/or v, such that scattering is suppressed in the non-relativistic limit. Specifically, O N 1 leads to spin-independent (SI) unsuppressed scattering, lead to SD momentum-suppressed scattering, which is typically unobservable. For the relativistic operators included in this study we give the dominant type of interaction they induce in the non-relativistic limit in Table 2.
The coefficients c N i (q 2 ) can be directly calculated from the Wilson coefficients of the relativistic operators at µ = 2 GeV. The explicit dependence on the momentum transfer q = √ 2m T E R is a result of two effects. First, under RG evolution some of the effective DMquark operators mix into the DM dipole operators Q (5) 1,2 (see Eq. (20)). These operators then induce long-range interactions, i.e. contributions to the c N i (q 2 ) that scale as q −2 . Since the momentum transfer can be very small in direct detection experiments, these contributions can be important in spite of their loop suppression. Second, the coefficients include nuclear form factors, obtained by evaluating expectation values of quark currents like N |qγ µ q|N . These form factors can be calculated in chiral perturbation theory and exhibit a pion pole for axial and pseudoscalar currents, i.e. a divergence for q → m π [100,101].
All of these effects are fully taken into account in DirectDM, which calculates the coefficients c N i (q 2 ) for given Wilson   and target element of interest. DDCalc also performs the velocity integrals needed for the calculation of the differential event rate, and the convolution with energy resolution and detector acceptance needed to predict signals in specific experiments: where M is the detector mass, T exp is the exposure time and φ(E R ) is the acceptance function. By combining DirectDM and DDCalc, we can obtain likelihoods for a wide range of direct detection experiments. In the present analysis, we include constraints from the most recent XENON1T analysis [93], LUX 2016 [88], PandaX 2016 [91] and 2017 [92] analyses, CDMSlite [84], CRESST-II [85] and CRESST-III [86], PICO-60 2017 [89] and 2019 [90], and DarkSide-50 [87].
The hadronic inputs to DirectDM v2.2.0 [67] were updated with the most recent N f = 2 + 1 lattice QCD results, following the FLAG quality requirements [107], see Table 3. All the inputs are evaluated at µ = 2 GeV. The hadronic matrix elements for protons and neutrons are related using isospin conservation.
For operators with vector quark currents, the least well known are the hadronic matrix elements involving the strange quark, while the matrix elements for operators with u, d quark vector currents have negligible errors to the precision we are working with. Since the strange quark vector current vanishes at q 2 = 0, the first non-vanishing contribution is obtained only at next-to-leading order in the chiral expansion, and depends on the strange quark charge radius, r 2 s = −0.0045(14) fm 2 [104,105]. For the nuclear magnetic moment induced by the strange quark, µ s = −0.036 (21) [104,105], we inflate the errors according to the Particle Data Group prescription.
The scalar form factors at zero recoil are obtained from expressions in Ref. [103], namely where the upper (lower) sign is for the proton (neutron). We use a rather conservative estimate σ πN = (50 ± 15) MeV [101] that covers the spread between the lattice QCD [111][112][113][114][115][116][117][118]  The matrix elements of tensor currents are described by three sets of form factors, but only two, g q T and B q/N T,10 (0), enter the chirally leading expressions. For g q T , the only N f = 2 + 1 result from Ref. [117] does not satisfy the FLAG quality requirements, so we use the N f = 2 + 1 + 1 results from Ref. [106] instead; the difference between the N f = 2 + 1 and N f = 2 + 1 + 1 results is expected to be small. For B q/N T,10 (0), we use the results from the constituent quark model in Ref. [109].

Relic abundance of DM
The Early Universe time evolution of the number density of the χ particles, n χ , is governed by the Boltzmann equation [124] dn χ dt where n χ,eq is the number density in equilibrium, H(t) is the Hubble rate and σv rel is the thermally averaged cross-section times the relative (Møller) velocity, given where K 1,2 are the modified Bessel functions and v lab is the velocity of one of the annihilating (anti-)DM particles in the rest frame of the other (for a discussion, see also Ref. [125]). We stress that there is no additional factor of 1/2 in the above equations. However, the fact that DM consists of Dirac particles implies that the total contribution to the observed DM density is given by n χ + nχ = 2n χ (disregarding the possibility of an initial asymmetry [126]). We compute tree-level annihilation cross-sections using CalcHEP v3.6.27 [127,128], where the implementation of the four-fermion interactions is generated by GUM [55, 56] from UFO files via the tool ufo_to_mdl (described in App. B). To ensure the EFT picture is valid, we invalidate points where Λ ≤ 2m χ . We obtain the relic density of χ by numerically solving Eq. (28) at each parameter point, assuming the standard cosmological history 8 and using the routines implemented in DarkSUSY v6.2.2 [130,131] via DarkBit. We then compare the prediction to the relic density constraint from Planck 2018: Ω DM h 2 = 0.120 ± 0.001 [98]. We include a 1% theoretical error on the computed values of the relic density, which we combine in quadrature with the observed error on the Planck measured value. More details on this prescription can be found in Refs. [48,52].
We note that our uncertainty estimate does not include uncertainties in the calculation of the annihilation cross-section very close to quark thresholds, which may be considerably larger. Moreover, our approach does not capture the potential effect of additional degrees of freedom on σv rel during freeze-out. The resulting effects, such as resonances or coannihilations could both increase and decrease the resulting value of Ω χ (see e.g. Ref. [132,133]), so the relic density constraint should be interpreted with care for Λ ∼ 2m χ , i.e. close to the EFT validity boundary (see Sec. 2.2).
The very nature of the EFT construction implies additional degrees of freedom above the energy scale Λ. Given the potential for a rich dark sector containing χ, and in particular, the possibility of additional DM candidates not captured by the EFT, we will by default not demand that the particle χ constitutes all of the observed DM, i.e. we allow for the possibility of other DM species to contribute to the observed relic density. In practice, this means that we modify the relic density constraint in such a way that the likelihood is flat if the predicted value is smaller than the observed one. In this case, we rescale all predicted direct and indirect detection signals by and f 2 χ , respectively. In doing so, we assume that the fraction f χ is the same in all astrophysical systems and that any additional DM population does not contribute to signals in these experiments. In a second set of scans we then impose a stricter requirement, namely that the DM particle under consideration saturates the DM relic abundance (f χ ≈ 1) rather than imposing the relic density as an upper bound (f χ ≤ 1). 9

Indirect detection with gamma rays
If DM is held in thermal equilibrium in the early universe via collisions with SM particles, then it can still annihilate today, especially in regions of high DM density. As with the relic abundance calculation, in order for the effective picture to hold for DM annihilation, we must impose Λ > 2m χ .
Gamma rays from dwarf spheroidal galaxies (dSphs) are a particularly robust way of constraining annihilation signals from DM [134]. In general, for a given energy bin i, the DM-induced γ-ray flux from target k can be written in the factorised form Φ i · J k , where details of the particle physics processes are encoded in Φ i , and details of the astrophysics are encoded in J k . See the DarkBit manual [52] for more details.
In general, only operators that lead to s-wave annihilation (Q (6) 1,3,4,q , Q 2,4,6,8,9,10,q ) give rise to observable gamma-ray signals; see for instance, Table 2. For the operators Q (6) 2,q and Q 1,3,5,7,q , the leading contribution to the annihilation cross-section is p-wave suppressed, i.e. proportional to v 2 rel . As DM in dSphs is extremely cold, with v 2 1/2 ∼ 10 −4 , this factor is very small, and the resulting limits are exceedingly weak. We therefore neglect p-wave contributions to all annihilation processes here.
For s-wave annihilation, one obtains where f χ is the DM fraction defined in eq. (30), (σv) 0,j denotes the zero-velocity limit of the cross-section for χχ → j and N γ,j is the number of photons, per annihilation, resulting from the final state channel j. The prefactor 1/4 accounts for the Dirac nature of the DM particles (under the assumption that n χ = nχ). Again, we use CalcHEP to compute annihilation cross-sections, with the CalcHEP model files generated by ufo_to_mdl via GUM (see App. B). The photon yields dN γ,j /dE used in DarkBit are based on tabulated Pythia runs, as provided by DarkSUSY.
The J-factor for each dSph k is simply the line-ofsight integral over the DM distribution assuming an NFW density profile and the solid angle Ω, where D k is the distance to the dSph. In our analysis we use the Pass-8 combined analysis of 15 dSphs after 6 years of Fermi-LAT data [96]. We use the gamLike v1.0.1 interface within DarkBit [52] to compute the likelihood for the gamma-ray observations, ln L exp , constructed from the product Φ i · J k and summed over all targets and energy bins, We also include a contribution from profiling over the J-factors of each dSph, ln L J = k ln L(J k ) [52,96], such that the full likelihood reads Gamma rays from the Galactic centre region provide a promising complementary way of constraining a signal from annihilating DM. While the J-factor is expected to be significantly higher than for dSphs, however, this conclusion is largely based on the result of numerical simulations of gravitational clustering rather than on the direct analysis of kinematical data. The reason for this is that the gravitational potential within the solar circle is dominated by baryons, not by DM, which adds additional uncertainty due to a dominant component of astrophysical gamma rays from this target region. As a result, Galactic centre observations with Fermi-LAT are somewhat less competitive than the dSph limits discussed above [135]. The upcoming Cherenkov Telescope Array (CTA), on the other hand, has a good chance of probing thermally produced DM up to particle masses of several TeV [136]. We will not include the projected CTA likelihoods in our scans, but indicate the reach of CTA when discussing our results.

Solar capture
The presence of non-zero elastic scattering cross-sections with nuclei combined with self-annihilation to heavy SM states leads to an additional, unique signature of DM in the form of high-energy neutrinos from the Sun. If Milky Way DM scatters with solar nuclei and loses enough kinetic energy to fall below the local escape velocity, it will become gravitationally bound. As long as it is above the evaporation mass threshold 4 GeV, captured DM will thermalize in a small region near the solar centre, and annihilate to SM products which then produce neutrinos via regular decay processes. These are distinct from the neutrinos from Solar fusion, as they are expected to have much higher energies than the ∼ MeV scales of fusion processes. Leading constraints have been obtained by Super-Kamiokande down to a few GeV [137], and by the IceCube South Pole Neutrino Observatory, between 20 and 10 4 GeV [138]. For typical annihilation cross-sections, the captured DM population reaches an equilibrium that is determined by the capture rate. For each likelihood evaluation, we obtain the nonrelativistic effective operators (Eq. (25)) as described in Sec. 3.1, using DirectDM to obtain the non-relativistic Wilson coefficients. These are passed to the public code Capt'n General [139], which computes the DM capture rate via the integral over the solar radius r and DM halo velocity u: where w(r) = u 2 + v 2 esc, (r) is the DM velocity at position r, and is the probability of scattering from velocity w to a velocity less than the local Solar escape velocity v esc, (r), dσ i /dE R is the DM-nucleus scattering cross-section, n i (r) is the number density of species i with atomic mass m N,i , and µ i = m χ /m N,i . Version 2.1 of Capt'n General uses the method described in detail in Ref. [140], separating the DM-nucleus cross-section into factors proportional to non-relativistic Wilson coefficients, powers of w and exchanged momentum q, and operator-dependent nuclear response functions computed in Ref. [140] for the 16 most abundant elements in the Sun. Solar parameters are based on the Barcelona Group's AGSS09ph Standard Solar Model [141,142]. Annihilation cross-sections are computed as described in Sec. 3.2, via CalcHEP. Once the equilibrium population of DM in the Sun has been obtained, crosssections and annihilation rates are passed to DarkSUSY, which computes the neutrino yields as a function of energy. These are finally passed to nulike v1.0.9 [143,144], which computes event-level likelihoods based on a reanalysis of the 79-string IceCube search for DM annihilation in the Sun [97].

Cosmic Microwave Background
Additional constraints on the DM annihilation crosssection arise from the early universe, more specifically from observations of the Cosmic Microwave Background (CMB). Annihilating DM particles inject energy into the primordial plasma, which affects the reionisation history and alters the optical depth τ . The magnitude of this effect depends on the specific annihilation channel and how efficiently the injected energy is deposited. These details can be encoded in an effective efficiency coefficient f eff , which depends on the injected yields of photons, electrons and positrons, and thus on the DM mass and its branching ratios into different final states [145]. The CMB is then sensitive to the following parameter combination: where σv rel ≈ (σv) 0 to a very good approximation during recombination; we thus also neglect p-wave contributions to all annihilation processes here. In order to calculate p ann for a given parameter point, one first needs to calculate the injected spectrum of photons, electrons and positrons and then convolve the result with suitable transfer functions that link the energy injection rate to the energy deposition rate [146]. The first part of this calculation has been automated within DarkSUSY and is accessible via DarkBit. The second part relies on DarkAges [147] (which is part of the ExoCLASS branch of CLASS) and is accessible via CosmoBit [54], see App. C for further details. 10 10 It is noteworthy that DarkAges calculates f eff (z) as a redshiftdependent function instead of a single redshift-independent coefficient f eff , as it is implicitly assumed in Eq. (37). In order to compress the function f eff (z) into this coefficient, it is convolved with a weighting function W (z) that encodes the CMB sensitivity to energy injection through s-wave annihilation as a function of redshift [145].
As the Planck collaboration only quotes the 95% credible interval for p ann [98], the remaining challenge is to obtain a likelihood for p ann from cosmological data. Although this likelihood can, in principle, be calculated for each parameter point individually using the Cos-moBit interface to CLASS and the Planck likelihoods, carrying out such a large number of calculations would be prohibitively slow, in particular if the cosmological parameters of the ΛCDM model are to be varied simultaneously. In the present work, we therefore adopt a simpler approach, where we first calculate the likelihood when varying p ann (while profiling over the ΛCDM and cosmological nuisance parameters). This approach yields where p 28 ann ≡ p ann / 10 −28 cm 3 s −1 GeV −1 . In arriving at this result, we have included the Planck TT,TE,EE+lowE+lensing likelihoods (using the 'lite' likelihood for multipoles ≥ 30, which only require one additional nuisance parameter [148]), as well as the BAO data of 6dF [149], SDSS DR7 MGS [150], and the SDSS BOSS DR12 galaxy sample [151]. This profile likelihood, which reproduces the 95% credible interval obtained by the Planck collaboration [98], can then be used in all subsequent scans, so that only p ann needs to be calculated for each parameter point and it is no longer necessary to call CLASS or plc.

Charged cosmic rays
Finally, DM particles annihilating in the Galactic halo also produce positrons, antiprotons and, to a lesser degree, heavier anti-nuclei that could in principle be observed in the spectrum of charged cosmic rays. Positrons quickly lose their energy through synchrotron radiation, and are thus a robust probe of exotic contributions from the local Galactic environment; the resulting bounds on DM annihilating to quarks or gluons are, however, much weaker than the other indirect detection constraints discussed here [152]. 11 Anti-nuclei, on the other hand, probe a significant fraction of the entire Galactic halo because energy losses are much less efficient in this case. For antiprotons, this generally leads to competitive constraints on DM annihilation signals [155][156][157], but it also means that such bounds necessarily strongly depend on uncertainties relating to modelling the production and propagation of cosmic rays in the Galactic 11 We note that this conclusion would be radically different for unsuppressed direct annihilation to leptons [153,154], which would result from leptonic operators analogous to Q halo. In addition to the dozen (or more) free parameters in the diffusion-reacceleration equations, there exist significant uncertainties on the energy dependence of the nuclear cross-sections responsible for the conventional antiproton flux [158] and possible correlated systematics [159]. A full statistical analysis, which would require a treatment of the large number of (effective) propagation parameters as nuisance parameters in our scans, is prohibitive in terms of computational costs [160] and hence beyond the scope of this work.

Collider physics
The effective operators defined in Sec. 2 allow for the pair production of WIMPs in the proton-proton collisions at the LHC. If one of the incoming partons radiates a jet through initial state radiation (ISR), one can observe the process pp → χχj as a single jet associated with missing transverse energy ( / E T ). In this study, we include the CMS [95] and ATLAS [94] monojet analyses based on 36 fb −1 and 139 fb −1 of data from Run II, respectively. ATLAS and CMS have performed a number of further searches for other types of ISR, leading for example to mono-photon signatures, but these are known to give weaker bounds on DM EFTs than monojet searches [24,161,162].
The expected number of events in a given bin of the / E T distribution is where L = 36 fb −1 or 139 fb −1 is the total integrated luminosity, σ the total production cross-section and the factor ( A) is the efficiency times acceptance for passing the kinematic selection requirements for the analysis. Both σ and ( A) can be obtained via Monte Carlo simulation, but given the dimensionality of the DM EFT parameter space it is computationally too expensive to perform these simulations on the fly during the parameter scan, as would be the standard approach to collider simulations within ColliderBit in GAMBIT.
Starting from UFO files generated using FeynRules v2.0 [163], we have therefore produced separate interpolations of σ and A based on the output of Monte Carlo simulations with MadGraph_aMC@NLO v2.6.6 [164] (v2.9.2) for the CMS (ATLAS) analysis, interfaced to Pythia v8.1 [165] for parton showering and hadronization. The matching between MadGraph and Pythia is performed according to the CKKW prescription, and the detector response is simulated using Delphes v3.4.2 [166]. The ColliderBit code extension that enables σ and ( A) interpolations to be used as an alternative to direct Monte Carlo simulation will be generalised and documented in the next major version of ColliderBit.
We only include the dimension-6 and 7 EFT operators (C (6) i and C (7) i=1,...,4 ) which are relevant for collider searches. Other operators give a negligible contribution due to either being suppressed by the parton distribution functions (in the case of heavy quarks), or by a factor of the fermion mass (small in the case of light quarks).
To reduce the computation time for our study, we generate events in discrete grids of the Wilson coefficients and DM mass. Separate grids are defined for each set of operators that do not interfere, such that the total number of events will simply be the sum of the contributions calculated from each grid. At dimension-6, there is interference between operators Q 1,q /Q (6) 4,q and Q 2,q /Q 3,q . For these Wilson coefficients, we parametrize the tabulated grids in terms of a mixing angle θ, defined via C 1,2 = sin θ and C 3,4 = cos θ. The CMS and ATLAS analyses have 22 and 13 exclusive signal regions, respectively, corresponding to the individual bins in the missing transverse energy distributions. As discussed below, the publicly available information makes it possible to combine all signal regions for the CMS analysis, while for the ATLAS analysis, only a single signal region can be used at once. To maximize the sensitivity of the ATLAS analysis, we combine the three highest missing energy bins, for which systematic uncertainties in the background estimation (and hence their correlations) are negligible, such that the highest bin in our analysis corresponds to all events with / E T > 1000 GeV. 12 Once the predicted yields for all bins have been evaluated, taking into account the EFT validity constraint as described in Sec. 2.2, we compute a likelihood for each analysis as follows.
For the CMS analysis, we follow the "simplified likelihood" method [167], since the required covariance matrix was published by CMS. In this approach, the full experimental likelihood function is approximated by a standard convolved Poisson-Gaussian form, with the systematic uncertainties on the background predictions treated as correlated Gaussian distributions: For each signal region i, the observed yield, expected signal yield and expected background yield are given by n i , s i and b i , respectively. The deviation from the nominal expected yield due to systematic uncertainties 12 We note that this combination also reduces the impact of a local ∼ 2.5σ excess in the third-highest bin, which would otherwise strongly bias our analysis.
is given by γ i . The correlations between the different γ i are encoded in the covariance matrix Σ provided by CMS, where we also add the signal yield uncertainties in quadrature along the diagonal. We follow the procedure in Ref. [167] in treating the γ i nuisance parameters as linear corrections to the expected yields. For every point in our scans of the DM EFT parameter space, we profile Eq. (39) over the 22 nuisance parameters in γ to obtain a likelihood solely in terms of the set of DM EFT signal estimates s: In the case of the ATLAS analysis, for which such a covariance matrix is not available, the conservative course of action is to calculate a likelihood using only the signal region with the best expected sensitivity. The ATLAS likelihood is therefore given by where L ATLAS (s i ,γ i ) is the single-bin equivalent of Eq. (39), and i refers to the signal region with the best expected sensitivity, i.e. the signal region that would give the lowest likelihood in the case n i = b i . The total LHC log-likelihood is then given by ln L LHC = ln L CMS + ln L ATLAS . However, due to the per-point signal region selection required in the evaluation of ln L ATLAS , the variation in typical yields between the different signal regions would manifest as a large variation in the effective likelihood normalization between different parameter points. To avoid this we follow the standard approach in ColliderBit of using the log-likelihood difference ∆ ln L LHC = ln L LHC (s) − ln L LHC (s = 0) (42) as the LHC log-likelihood contribution in the parameter scan [168]. When presenting the results of a global fit we identify the maximum-likelihood point Θ best-fit in the DM EFT parameter space and map out the 1σ and 2σ confidence regions defined using the likelihood ratio L(Θ)/L(Θ best-fit ). Thus, in cases where some region of the DM EFT parameter space can accommodate a modest excess in the collider data, other DM EFT parameter regions that might still perform better than the SM, or that are experimentally indistinguishable from SM, can appear as excluded. While this is perfectly reasonable, given that the comparison is to the best-fit DM EFT point and not to the SM expectation, it is also interesting to study the global fit results under the assumption that mild excesses in the collider data indeed do not originate from a true new physics signal. A simple and pragmatic approach is then to replace ∆ ln L LHC with a capped version, This will assign the same log-likelihood value, ∆ ln L cap LHC = 0, for all DM EFT parameter points whose prediction fit the collider data as well as, or better than, the SM prediction (s = 0) does. Thus, analogous to how exclusion limits from LHC searches are constructed to only exclude new physics scenarios that predict too many signal events, the capped likelihood only penalizes parameter points for performing worse than the background-only scenario. The result obtained from using ∆ ln L cap LHC in a fit is therefore close to the result one would obtain by constructing a joint exclusion limit for the LHC searches, and applying this limit as a hard cut on the parameter space favoured by the other observables. The main difference is that the capped LHC likelihood incorporates a continuous likelihood penalty. 13 A more detailed introduction to the capped likelihood construction can be found in Ref. [169].
Below we will present some results using this capped LHC likelihood, and some using the full LHC likelihood in Eq. (42). In light of the discussion above, the two sets of results should be interpreted as answering slightly different questions: The fit results with the full LHC likelihood show what DM EFT scenario is in best agreement with the complete set of current data, and how much worse other DM EFT scenarios perform in comparison. The results with the capped LHC likelihood map out the DM EFT parameter space that is preferred by the non-collider observables and not excluded by a combination of the LHC searches.

Nuisance parameter likelihoods
In our scans we also vary a set of relevant nuisance parameters related to the DM observables and SM measurements. Most of these nuisance parameters are directly constrained by dedicated measurements, which we include through appropriate likelihood functions. In some cases, however, several conflicting measurements exist, indicating additional systematic uncertainties in the methodology. In these cases we constrain the nuisance parameters through effective likelihoods intended to give a conservative constraint on the allowed ranges. The nuisance parameters and 3σ ranges used in this  study are summarised in Table 4. We briefly cover each nuisance likelihood in turn below. We follow the default prescription in DarkBit for the local DM density ρ 0 , where the likelihood is given by a log-normal distribution with central value ρ 0 = 0.40 GeV cm −3 and error σ ρ0 = 0.15 GeV cm −3 . We scan over an asymmetric range in ρ 0 to reflect the log-normal distribution -see Ref. [52] for more details.
We follow the same treatment of the Milky Way halo as in the GAMBIT Higgs portal study [110]. We utilise Gaussian likelihoods for parameters describing the Maxwell-Boltzmann velocity distribution, specifically the peak of the distribution v peak = 240 ± 8 km s −1 [170], and the galactic escape velocity v esc = 528 ± 25 km s −1 , based on the Gaia data [171].
We employ a Gaussian likelihood for the running top quark mass in the MS scheme with a central value m t (m t ) = 162.9 GeV and an error 2.0 GeV [172]. 14 The top pole mass (m pole t ) is then computed using the following formula: We use only the one-loop QCD corrections in this shift in order to be consistent with the procedure carried out in Ref.
[172]. We have checked that the above expression gives the expected result for the top pole mass and matches well with Ref.
[172]. For direct detection, we employ nuisance parameter likelihoods for a number of hadronic input parameters that are used to evaluate form factors at the nuclear scale. Specifically, we use a product of four Gaussian likelihoods to include the constraints on σ πN , ∆s, g s T and r 2 s quoted in put parameters are fixed to the central values given in Table 3.

Results
We now present the results obtained from comprehensive scans of the parameter space introduced above. These scans were carried out with the differential evolution sampler Diver v1.4.0 [173] using a population of 5 × 10 4 and a convergence threshold of either 10 −5 or 3 × 10 −5 . As we will analyse our scan results using profile likelihood maps, the sole aim of the scans is to map out the likelihood function in sufficient detail across the highlikelihood regions of parameter space. In particular, no statistical interpretation is associated with the density of parameter samples, and we can therefore combine samples from scans that use different metrics on the parameter space. To ensure that all parameter regions are properly explored, we perform two different types of scans: -Full: We explore DM masses up to the unitarity bound (5 GeV < m χ < 150 TeV and 20 GeV < Λ < 300 TeV). 15 In these scans, m χ and Λ are scanned on a logarithmic scale, while the Wilson coefficients are scanned on both a linear and a logarithmic scale (i.e. we combine the samples from both scanning strategies to achieve a thorough exploration of the whole parameter space). -Restricted: We consider the parameter region where experimental constraints are most relevant 15 We note that for the largest values of mχ and Λ that we consider in these scans our approach of specifying all operators in the broken phase of electroweak symmetry and ignoring the effects of running between µ = Λ and µ = m Z becomes questionable. The constraints that we obtain above the TeV scale are therefore only approximate and should be interpreted with care.
(m χ < 500 GeV and Λ < 2 TeV). In these scans the DM mass is scanned on a linear scale, the scale Λ on a logarithmic scale and the Wilson coefficients on a scale that is logarithmic on [−4π, −10 −6 ], linear on [−10 −6 , 10 −6 ] and logarithmic on [10 −6 , 4π]. This approach was found to achieve the optimum resolution of the LHC constraints while simultaneously ensuring that enough viable samples are also found for small Λ when some or all of the Wilson coefficients are tightly constrained.
All nuisance parameters are scanned on a linear scale. In the first set of scans, we fix the Wilson coefficients for all dimension-7 operators to zero, so that there are 6 model parameters and 8 nuisance parameters. The second set of scans then includes all 14 Wilson coefficients, bringing the total number of parameters up to 24. We furthermore consider a number of variations in the constraints that we include in our scans: -We perform scans where the DM particle is allowed to be a sub-component (f χ ≤ 1) and scans where we require that the DM relic density be saturated (f χ ≈ 1), see Sec. 3.2; -We perform scans with both the capped LHC likelihood and the full LHC likelihood (see Sec. 3.5); -When considering the full LHC likelihood, we furthermore apply two different prescriptions for imposing the EFT validity: a hard cut-off and a smooth cut-off (see Sec. 2.2).
Unless explicitly stated otherwise, our default choices for the discussion below are to allow a DM sub-component and consider the capped LHC likelihood with a hard cut-off.

Capped LHC likelihood
Let us begin with the case that the LHC likelihood is capped, i.e. it cannot exceed the likelihood of the background-only hypothesis. We first consider only dimension-6 operators with different requirements for the DM relic density, and then also include dimension-7 operators.

Dimension-6 operators only (relic density upper bound)
Our main results for this case are shown in Fig. 2 in terms of the DM mass and the new physics scale Λ. The left panel corresponds to the full parameter range, whereas the right panel provides a closer look at the most interesting parameter region. We find a large viable parameter space but also a number of notable features. For large values of m χ and Λ, the allowed parameter space is determined by the EFT validity requirement Λ > 2m χ and the relic density requirement which, combined with the perturbativity bound on the Wilson coefficients, implies an upper bound on Λ for given m χ . These different constraints are compatible only for m χ < 150 TeV, implying an upper bound on the scale of new physics of Λ < 300 TeV. This limit corresponds to the well-known unitarity bound for thermal freeze-out [65]. The zoomed-in version in the right panel reveals a number of additional features. In the top-left corner (small m χ , large Λ), there are strong constraints from the LHC, which make it impossible to satisfy the relic density requirement. These constraints become weaker as Λ decreases and the EFT can only be trusted for smaller values of / E T . The various sharp features correspond to the points where Λ crosses the boundary of a specific / E T bin, leading to a jump in the likelihood. In our conservative approach, LHC constraints are completely absent for Λ < 200 GeV. Finally, we find that there is a slight upward fluctuation in Fermi-LAT data, which can be fitted for m χ = 5.0 GeV and f 2 χ σv 0 = 1.1 × 10 −27 cm 3 s −1 . 16 We emphasize that a great advantage of our approach is that we treat the new-physics scale Λ as an independent parameter, which is kept explicit in Fig. 2 (rather than being profiled out like the individual Wilson coefficients). This makes it possible in a straight-forward way to distinguish those parameter regions where the EFT predictions can be considered robust and those parameter regions where additional constraints may apply. As discussed in Sec. 2.2, the EFT is expected to be valid if Λ is sufficiently greater than the largest p T bin considered in the LHC analyses, i.e. Λ > 1.3 TeV. Conversely, for Λ < 200 GeV we conservatively suppress constraints from the LHC, such that the viable parameter regions found in this range must be interpreted with great care. For intermediate values of Λ, LHC constraints are being applied but may depend on the specific UV completion. Which of these parameter regions is considered most interesting depends on the specific context and is left to the reader.
A complementary perspective is provided in Fig. 3, which shows the allowed parameter regions in terms of the DM mass, the relic density and the rescaled annihilation cross-section. A number of additional features become apparent in these plots. First, for m χ 100 GeV it is impossible to saturate the observed DM relic density, Ω DM h 2 = 0.12, due to the combined constraints from direct and indirect detection experiments. However, these The shaded region (corresponding to Λ ≤ 2mχ) is excluded by the EFT validity requirement. In the right panel, the parameter ranges have been restricted to the most interesting region. Note that the position of the best-fit points in the two panels is somewhat arbitrary, as there is a degeneracy between Λ and C constraints are suppressed for DM sub-components, such that it is possible to have very small relic densities in this mass region. For m χ > 100 GeV (corresponding to Λ > 200 GeV), on the other hand, constraints from the LHC become relevant, which are not suppressed for DM sub-components. These constraints are then again relaxed for m χ 1 TeV as the LHC energy becomes insufficient to produce a pair of DM particles.
For m χ 1 TeV, we find that there is a direct correspondence between Ω χ h 2 and the rescaled annihilation cross-section f 2 χ σv 0 . This is because the operators that induce p-wave annihilations (in particular C (6) 2 ) are strongly constrained by the LHC and direct detection experiments, and the annihilation cross-section is therefore always dominated by the s-wave contribution. For larger DM masses, it becomes possible for the p-wave contribution to dominate the relic density calculation, such that the total annihilation cross-section is velocitydependent and becomes tiny in the present universe. While indirect detection experiments presently cannot probe the relevant parameter space for TeV-scale DM, it is worth stressing that CTA will be able to do so for operators that induce s-wave annihilation. We illustrate this in Fig. 3 by indicating the sensitivity of CTA to a DM signal from the galactic center [136] (for simplicity based on the assumption of bb final states, noting that any hadronic DM annihilation channel results in very similar gamma-ray spectra at these energies). We note that the CTA sensitivity indicated in Fig. 3 is based on assuming a standard Einasto profile as expected for WIMP DM; if the DM density in the galactic center is instead roughly constant, the sensitivity can worsen by up to about one order of magnitude [136].
Let us finally consider the allowed parameter space in terms of the Wilson coefficients. The coefficient C gives rise to spin-independent scattering, which is very strongly constrained by direct detection experiments. Thus, this coefficient is required to be so small that it cannot give a sizeable contribution to any other process. The coefficient C  2 , which gives rise to a momentum-suppressed spin-independent scattering (see Table 2).
In the right panel of Fig. 4, we show the allowed parameter region in terms of C (6) 3 , which induces scattering that is simultaneously momentum-suppressed and spin-dependent, such that direct detection constraints are very weak. Correspondingly, we find that this coeffi-  Fig. 2, we consider only dimension-6 operators and cap the LHC likelihood at the value of the background-only hypothesis. The solid red line in the middle panel denotes the "initial construction" projection sensitivity of Cherenkov Telescope Array (CTA) towards the Galactic Centre (GC) for the bb final state [136]. cient is largely unconstrained for Λ < 200 GeV. We also identify this coefficient as giving the main contribution for fitting the Fermi-LAT excess. For larger values of Λ, on the other hand, the constraints are very similar to the ones for C (6) 2,4 as the LHC only has limited sensitivity to distinguish the spin structure of the operators.

Dimension-6 operators only (relic density saturated)
Next we consider the case where the relic density constraint is imposed not only as an upper limit but as an actual measurement, i.e., the DM particle under consideration is required to account for all of the DM in the universe via the effective interactions that we consider. We show in Fig. 5 the allowed parameter space in the restricted m χ -Λ plane when considering a capped LHC likelihood, i.e. the same likelihoods as in Fig. 2 apart from the modified relic density requirement. As expected from the top row of Fig. 3, it is not possible to saturate the observed relic density for m χ 100 GeV. The reason is that for such small DM masses the relic density requirement is incompatible with Fermi-LAT and CMB bounds on the annihilation cross-section for operators that predict dominantly s-wave annihilation (Q (6) 1,q and Q (6) 3,q ), and incompatible with direct detection and LHC constraints for Q (6) 2,q and Q (6) 4,q . Constraints from direct and indirect detection experiments are also responsible for the preference for larger DM masses visible in Fig. 5. In particular, the Fermi-LAT likelihood pushes the best-fit point towards the boundary m χ = 500 GeV. We find the likelihood of the best-fit point to be slightly worse than for the background-only hypothesis: 2∆ ln L ≡ 2(ln L best−fit − ln L ideal ) = −0.5. Extending the range of the scan to larger DM masses would allow the model to fully evade the Fermi-LAT constraint. This would shift the bestfit point and the allowed parameter regions to slightly larger DM masses without changing the remaining conclusions (see also Fig. 3).
For a complementary view of the parameter space, we show in Fig. 6 the predicted number of signal events in the next-generation direct detection experiment LZ [174] as a function of the DM mass. Due to the various different operators contributing to the DM-nucleus scattering, the predicted number of signal events is a more useful quantity to consider than the DM-nucleon scattering cross-section at zero momentum transfer. The predicted number of events corresponds to nuclear recoil energies in the search window [6 keV, 30 keV] and assumes an exposure of 5.6 × 10 6 kg days and 50% acceptance for nuclear recoils (see Ref. [110] for details on our implementation of LZ). plane (right) for the restricted parameter ranges. As in Fig. 2, we only consider dimension-6 operators and cap the LHC likelihood at the value of the background-only hypothesis. The contour lines show the 1σ and 2σ confidence regions. Note that the position of the best-fit points (white stars) in the two panels is somewhat arbitrary, as there is a degeneracy between Λ and C  We find that of the order of 10 events are predicted around the best-fit point, which requires a non-zero contribution from the operator Q leading to spinindependent (but momentum-suppressed) scattering. However, the predicted number of events varies significantly within the allowed region of parameter space and can be as small as 0.1 at 68% confidence level. In this case the main contribution arises from the mixing of the operator Q as given in Eq. (21). 17 While such an event number is too low to be detected with next-generation experiments, it is still well above the neutrino background and should be observable with more ambitious future detectors such as DARWIN [175] or DarkSide-20k [176].
Another interesting approach would be to not perform a relic density calculation at all and simply assume that f χ = 1 is achieved through some modification of early universe cosmology. In this case it would also be possible to consider Λ < 2m χ since the calculation of the annihilation cross-section is unnecessary. However, since none of the other likelihoods that we consider give a strong preference for a DM signal, there would then be no lower bound on the interaction strength of the DM particle, i.e. it would be possible for all Wilson coefficients to vanish simultaneously. Hence we expect all combinations of m χ and Λ to be viable in this approach, and we do not explore this direction further in the present work.

Operators up to dimension 7 (relic density upper bound)
We now turn to the case where we simultaneously consider all dimension-6 operators as well as the dimension-7 operators involving DM particles and quarks or gluons introduced in Sec. 2. We remind the reader that we neglect additional dimension-7 operators involving Higgs bosons that would arise in theories respecting unbroken electroweak symmetry (which are phenomenologically irrelevant) as well as operators with derivative interactions (which largely give redundant information). Even with these restrictions our analysis requires 24-dimensional (16 model + 8 nuisance) parameter scans.
In Fig. 7, we show the allowed regions in the m χ -Λ plane (left) and in the m χ -Ω χ h 2 plane (right) when using the capped LHC likelihood. As before, we find that the parameter region at small m χ and Λ can fit the slight Fermi-LAT excess with best-fit values: m χ = 5.5 GeV and f 2 χ σv 0 = 1.9 × 10 −27 cm 3 s −1 . As the inclusion of additional parameters can only increase the profile likelihood, we expect the allowed regions of parameter space to be larger than the ones found above. Interestingly, the differences between the left panel of Fig. 7 and the right panel of Fig. 2 are rather minimal. In other words, the inclusion of the 10 additional dimension-7 operators does not open up new parameter space in terms of m χ and Λ. This is of course expected for the parameter region with large m χ and small Λ (bottom-right), which is excluded by the EFT validity constraint but surprising for the region with small m χ and large Λ (top-left), which is excluded by the combination of the LHC constraints and the relic density requirement.
The reason why this parameter space remains inaccessible is that the gluon operators Q 1-4 are strongly constrained by the LHC for Λ > 200 GeV and can therefore not contribute significantly to the annihilation crosssection. The quark operators Q 5-10,q , on the other hand, are unconstrained by the LHC, but for m χ < m t , the resulting annihilation cross-section is suppressed by a factor m 2 b m 2 χ /Λ 6 , and therefore too small to produce a relic abundance that evades the upper bound from the relic density requirement given the perturbativity bound on the Wilson coefficients. Comparing the right panel of Fig. 7 to the allowed parameter regions from Fig. 3 (indicated by the grey dashed lines) does however reveal a number of differences. First of all, it is now possible to saturate the relic density bound for small m χ (and small Λ), thanks to the contribution of Q 3,q and Q 7,q , which both give suppressed signals in direct and indirect detection experiments and are therefore largely unconstrained. Moreover, for m χ > m t , we find that the predicted relic abundance can be substantially smaller than for the case with only dimension-6 operators, thanks to the contribution from the dimension-7 DM-quark operators Q 5-10,q . The additional freedom in the annihilation cross-section also implies that the impact of imposing a strict relic density requirement is reduced compared to the case of dimension-6 operators only and will therefore not be discussed in further detail here.
We emphasize that global fits with 24 free parameters are computationally quite challenging, in particular when the best-fit region is not strongly constrained by data. As a result the contours in Fig. 7 are less smooth than for the case of dimension-6 operators only. This is particularly obvious in the right panel for DM masses around 150 GeV. In this region many operators are strongly constrained by LHC data while annihilations into top quarks are kinematically forbidden. This makes it challenging to find parameter points that satisfy the relic density constraint, leading to comparably poor sampling. We have confirmed explicitly that this is not a physical effect, i.e. the allowed parameter region should be smooth and extend to Ω χ h 2 = 0.12 everywhere.

Dimension-6 operators only (relic density upper bound)
We now move onto the case where the full (rather than capped) LHC likelihood is included in the scans. Fig. 8 shows the allowed parameter regions in terms of m χ and Λ for the case where we introduce a hard cut-off in the missing energy spectrum for / E T > Λ (left panel), and the case where we introduce a smooth cut-off (right panel), as discussed in Sec. 2.2. We see that in both cases, the results differ from Fig. 2, i.e. there is a preference for higher Λ values. This preference arises due to data excesses in a few high-/ E T bins in the ATLAS and CMS monojet searches.
The difference in the above two results can be understood as follows. For / E T < Λ, the missing energy spectrum arising from DM is harder than the background,  while for / E T > Λ, we either set it to zero or assume that it drops rapidly. Thus, the ratio of signal-to-background is largest for / E T ≈ Λ, enabling our model to (partially) fit local excesses in the data. This is illustrated in Fig. 9, which shows the missing energy spectra for background and signal in CMS when applying different EFT validity prescriptions. As seen in the distribution of pulls in the bottom panel, the CMS search observes a couple of 1σ-2σ data excesses in bins around / E ≈ 700 GeV (purple bars). By including a DM signal prediction on top of the SM background, these excesses can be reduced, thus reducing the pulls and improving the overall fit to the data (green bars). However, unless the signal spectrum dies off sufficiently fast above / E ≈ 700 GeV, the model will be penalized for causing larger pulls in the highest-/ E T bins, as seen for instance for the unmodified signal spectrum (lightest green bars, corresponding to a = 0).
For the case where we impose a hard cut-off (left panel in Fig. 8), we find (at the 1σ level) separate parameter regions preferred by the CMS analysis (Λ ≈ 700 GeV) and the ATLAS analysis (Λ 1 TeV), with the overall best-fit point corresponding to the latter and being preferred relative to the background-only hypothesis by 2∆ ln L = 2.2. When allowing for a smooth cut-off, on the other hand, the best-fit solution produces a partially improved fit to both excesses simultaneously, by suppressing the signal distribution approximately proportional to ( / E T /Λ) −1 . In this case, the best-fit point has 2∆ ln L = 2.6. 18 We refrain from translating these numbers into p-values, which would require extensive Monte Carlo simulations. For both choices of cut-off, the best-fit point predicts an annihilation cross-section that is slightly larger than the thermal cross-section, such that the DM particles in this case would only constitute a DM sub-component.
We emphasise that the preference for a non-zero signal contribution is to some degree an artefact of the way in which we have implemented the EFT validity requirement. Realistic UV completions typically do not introduce sharp features in the missing energy spectrum, making it harder to fit excesses observed in individual bins. Nevertheless, our findings emphasise the need to analyse missing energy searches at the LHC in terms of specific models in order to assess whether the signal preference found in the EFT approach can be recovered (at least partially) in a more complete setting.

Dimension-6 operators only (relic density saturated)
We have also run scans with the full LHC likelihood and requiring the DM relic density to be saturated (see Fig. 10). We find the expected changes with respect to Fig. 8, namely that small DM masses are disfavoured. For the case of a hard cut-off, the position of the bestfit point is unaffected, while for a smooth cut-off, it is pushed to slightly larger values of m χ and Λ. The respective preferences are reduced slightly to 2∆ ln L = 1.9 and 2.0. We also find that the best-fit point requires several Wilson coefficients to be non-zero. While the LHC signal can be fitted by either C (6) 2 or C (6) 4 , the relic density can only be reproduced with a contribution from C (6) 3 . This is because Q (6) 2,4 lead to suppressed annihilation rates in the early universe, compared to Q 3 , while Q 1 is strongly constrained by direct detection (see also Table 2). A summary of the various best-fit points from our scans with dimension-6 operators only is given in Table 5. We note that essentially all of our scans require a non-zero contribution from C (6) 3 at the best-fit point in order to satisfy the relic density requirement. This is an interesting finding given that this operator is present only for Dirac fermion DM but not for Majorana fermion DM. In other words, we expect our results to change considerably for the case of Majorana fermion DM. Satisfying the relic density constraint with dimension-6 operators only while evading experimental constraints will be very challenging in this case.

Operators up to dimension 7
In Fig. 11, we finally show the case where the full LHC likelihood is included when simultaneously considering all dimension-6 and dimension-7 operators, using either a hard cut-off (left) or profiling over possible smooth cut-offs (right). In the former case we find that the result looks very similar to the case of dimension-6 operators only (left panel of Fig. 8) and also the likelihood at the best-fit point is very similar. In the latter case we find that it is now possible to simultaneously accommodate the upward fluctuations in the Fermi-LAT data (as in Fig. 2) and in the LHC data (as in Fig. 8). Doing so requires a small new-physics scale Λ ∼ 80 GeV together with a rather soft cut-off a ≈ 1.7 of the / E T spectrum above Λ. The resulting best-fit point has 2∆ ln L = 2.9, which is the highest likelihood found in any of our scans.
A closer analysis reveals that the contribution of the dimension-6 operators is in fact not necessary to accommodate the small LHC excesses, because sufficiently large contributions can also be obtained from the gluon operators. For example, the operator Q 4 is essentially unconstrained by direct detection and can induce sizeable LHC signals if C 4 takes values close to the perturbativity bound. While it is challenging to satisfy the relic density requirement using only gluon operators, the allowed parameter space expands substantially when including a contribution from the dimension-7 DM-quark operators Q 5-8,q . As a result, the allowed regions in m χ -Λ parameter space look very similar to the ones shown in Fig. 11 even when the Wilson coefficients of all dimension-6 operators are set to zero. For the same reason we expect no significant difference between Dirac and Majorana DM particles in this case. This complex interplay between different operators only becomes apparent in a global analysis and would be missed when studying individual operators separately.

Conclusions and outlook
In this work we have presented the first global analysis of the full set of effective operators up to dimension 7 involving a Dirac fermion DM particle and quarks or gluons. Key to enabling such an analysis were a number of technical developments: -We have fully automated the calculation of direct detection constraints, including mixing under RG evolution and matching onto non-relativistic effective operators at the hadronic scale, and indirect detection constraints, including cosmological constraints on energy injection; -We have adopted a novel approach to address the issue of EFT validity at the LHC. Rather than performing a simple truncation procedure, we introduce 2 ) 2 + (C 2 ) 2 + (C 4 ) 2 /Λ 2 = 0.23 Table 5: Best-fit points from our various scans involving dimension-6 operators with restricted parameter ranges (5 GeV ≤ mχ ≤ 500 GeV and 20 GeV ≤ Λ ≤ 2 TeV). For most scans, there are degeneracies between different parameters around the best-fit point.
In these cases, we only quote the combination that is well-constrained rather than each parameter individually. Parameters not stated explicitly are compatible with zero. a smooth cut-off for / E T > Λ and treat this parameter as a nuisance parameter to ensure that no artificially strong exclusions arise from the tails of the predicted distributions; -We employ highly efficient likelihood calculations and sampling algorithms that make it possible to scan over up to 24 parameters (the DM mass m χ , the new physics scale Λ, 14 Wilson coefficients and 8 nuisance parameters).
In combination, these developments enable us, for the first time, to include interference effects between different operators in all parts of the analysis.
Our main result is that it is typically possible to suppress the scattering and annihilation cross-sections in the non-relativistic limit, and thereby evade direct and indirect detection constraints while satisfying the relic density requirement. Doing so does not require finely tuned cancellations or interference effects but is a direct consequence of the spin structure of the operators that we consider. The LHC, however, plays a special role, because the production of relativistic DM particles is less sensitive to the specific spin structure of the operator. As a result, we find generally strong constraints on small DM masses and large Λ, both for the case of dimension-6 operators only and also when including dimension-7 operators. Moreover, when allowing excesses in individual LHC bins to be fitted (rather than artificially capping the LHC likelihood), we find a slight preference for a DM signal with a relatively low new physics scale. Given that the magnitude of this excess is sensitive to the precise EFT validity prescription that we adopt, we have not attempted to quantify its significance within the EFT.
We find that it is typically not necessary to have simultaneous contributions from many different operators in order to find viable regions of parameter space. Indeed, large viable regions of parameter space are found both for the case when we consider only dimension-6 operators and only dimension-7 operators. These sets of operators can easily be generated by integrating out a heavy mediator with spin 1 or spin 0, respectively. However, we typically do require sizeable contributions from operators that violate parity and/or CP, reflecting the pressure on the simplest WIMP models from the non-observation of a DM signal in direct and indirect detection experiments (see Ref. [110] for a similar discussion in the context of Higgs portal models).
A particularly interesting observation is that it is generally not possible to have a large hierarchy between the DM mass and the new physics scale without violating the relic density requirement. In particular, for m χ 100 GeV, constraints from the LHC require Λ 200 GeV, meaning that the EFT is no longer valid at LHC energies and additional new degrees of freedom should be kinematically accessible. Moreover, the wellknown unitarity bound on the DM mass implies a robust upper bound on the scale of new physics of the order of 300 TeV. We also note that for masses in the TeV range CTA will have a unique chance of probing part of the currently inaccessible parameter space that is spanned between the EFT validity and the relic density constraints.
We emphasise that it is generally possible for the DM particle under consideration to constitute only a DM sub-component (in which case, constraints from direct and indirect detection experiments are correspondingly suppressed), but large regions of viable parameter space also remain when requiring the relic density to be saturated. In future studies, it will be interesting to modify the way in which the relic density calculation is included. For example, one could consider an initial particle-antiparticle asymmetry in the dark sector, which would make it possible to saturate the relic density in parameter regions that would normally predict an underabundance, while at the same time suppressing constraints from indirect detection experiments. A more radical approach would be to not perform a relic density calculation at all and simply assume that the observed relic abundance (with f χ = 1) is achieved through some unspecified modification of standard cosmology. A detailed analysis of direct detection constraints on such a scenario is in preparation.
An exciting direction for future investigation is to embed the EFTs considered here into a more complete approach based on UV-complete (or simplified) models. Almost all of the machinery developed for the present work will also be directly applicable in this case. The main difference arises in the interpretation of the LHC signals. If the mediator of the DM interactions is kinematically accessible at LHC energies, it will be essential to not only consider the resulting changes in the missing energy spectra, but also additional signatures arising from visible decays of the mediator [177,178] (see Ref. [179] for a recent discussion of how to connect DM EFTs and UV-complete models). Furthermore, close to the EFT validity boundary the presence of the mediator will also modify the results of the relic density calculation, thus affecting the target couplings for these signals. It will also be interesting to see to what extent the slight LHC excesses can be accommodated in such a set-up.
Another important extension of the present work will be to also consider operators coupling DM to leptons as well as electroweak gauge and Higgs bosons in order to embed our approach into a framework that respects the unbroken electroweak symmetry. Given that the relevant RG evolution is known (and already implemented in DirectDM) and that the relevant annihilation cross-sections and injection spectra can be calculated automatically, such an extension does not pose any conceptual difficulties regarding direct or indirect detection constraints and relic density calculations. Again, the most challenging part will be to include all relevant collider constraints (which in this case stem also from LEP). Given that these constraints are typically weaker than the corresponding ones for quarks, it will be interesting to see whether some of the conclusions found in the present work can be relaxed and additional viable parameter space opens up.
Finally, it will be very interesting to consider DM EFTs with non-trivial flavour structure, for example with couplings predominantly to the third generation. In such a set-up, one generally expects sizeable flavourchanging neutral currents and hence it will be essential to connect the EFTs used to study DM to the ones employed in flavour physics. Such a study would be particularly exciting given the recently observed anomalies in various flavour observables (see e.g. Refs. [180][181][182]). Moreover, the effects of electroweak operator mixing on the direct detection bounds are expected to be much more pronounced in such scenarios.
Of course, the most important outstanding task is to collect more data that may shed light on the nature of DM. Upcoming LHC analyses will improve the sensitivity to missing energy signatures of DM, the next generation of direct detection experiments [174,183,184] will be able to probe substantially smaller scattering cross-sections, and ongoing [185][186][187] and planned [136] indirect detection experiments will probe the freeze-out paradigm with unprecedented precision. Our present work has shown that this effort is highly worthwhile given the wide regions of parameter space that cannot currently be excluded in a model-independent way. Reducing the vast number of viable possibilities to explain DM therefore remains a key challenge for years to come.
to interface a new model to DirectDM. For more background on the technical aspects of the GAMBIT framework, please refer to the original GAMBIT manual [48], and the GUM paper [55,56].
DirectDM matches Wilson coefficients of a relativistic EFT onto a non-relativistic EFT valid at the nuclear scale. The GAMBIT implementation interfaces with the Python version of this package.
Relativistic Wilson coefficients can be defined at the 3-, 4-or 5-quark scale, with the capability DD_rel_WCs_flavscheme. For a given model, a new module function providing this capability should be written, returning the type map_str_dbl (std::map<std::string,double>). Once this capability has been fulfilled, GAMBIT uses the module function DD_nonrel_WCs_flavscheme to call the DirectDM backend via the convenience function get_NR_WCs_flav. This provides the capability DD_nonrel_WCs which can be connected to the DDCalc backend.
This module function providing the capability DD_nonrel_WCs depends on the capability WIMP_properties, of native GAMBIT type WIMPprops. WIMP_properties supplies the particle information about the WIMP candidate, such as its spin, mass, and whether or not it is self-conjugate, extracted from the particle database and either the spectrum or model parameters.
As an example, consider a simplified model where a vector mediator governs the interaction between dtype quarks and a fermionic DM candidate χ, with the following interaction Lagrangian, The model implementation within GAMBIT will contain four free parameters: the couplings g χ and g b , the DM mass m χ , and the mediator mass m V . The model definition for the above simplified model looks like: The operator in DirectDM corresponding to this interaction is Q 1,q = (χγ µ χ)(qγ µ q). We identify the relevant coefficient to pass to DirectDM as g χ g b /m 2 V . This is simply implemented in DarkBit by the following source code: For a full definition of the operator basis used in Di-rectDM, we refer the reader to Refs. [67,68].
When DirectDM is used, the user must also scan over the model nuclear_params_ChPT_sigmapiN, which contains (nuisance) parameters used in the matching and running routines in DirectDM. These are defined in Table 1 of Ref. [67]. We provide a YAML file containing the default values used in DirectDM (see the nuclear_params_ChPT_sigmapiN_fixed.yaml file in the $GAMBIT/yaml_files/include/ directory).

Appendix B UFO to CalcHEP
ufo_to_mdl is a simple Python tool distributed with GAMBIT v2.1 and above, and is integrated in the GUM framework. ufo_to_mdl is located at $GAMBIT/gum/src/ufo_to_mdl.py. It can also be run as a standalone tool, using either Python2 or Python3. Below we briefly describe the motivation for ufo_to_mdl and how to use it.
The purpose of ufo_to_mdl is to generate CalcHEP input (.mdl files) from UFO files. The motivation for this tool's creation is that FeynRules does not generate four-fermion CalcHEP output, but it can create such output for MadGraph. In fact, at the time of writing, LanHEP [189] is the only package that supports automatic generation of four-fermion contact interactions for CalcHEP files. ufo_to_mdl allows the user to study fourfermion interactions using CalcHEP (and correspondingly, micrOMEGAs), effectively creating a pathway from FeynRules to CalcHEP for effective theories of this kind. In the context of GAMBIT and the GUM pipeline, ufo_ to_mdl allows the user to study EFTs of DM using the routines provided by micrOMEGAs and CalcHEP inside of the GAMBIT framework, such as relic density calculations, direct detection rates, and indirect detection via the Process Catalogue (see the DarkBit manual [52] for details).
Usage of ufo_to_mdl is straightforward. There are two modes ufo_to_mdl can be operated in: comparison mode and conversion mode. The mode integrated into the GUM pipeline is the comparison mode, which compares two directories containing .ufo and .mdl files generated by FeynRules: python ufo_to_mdl.py <UFODir> <MDLDir> This ensures that all vertices in the MadGraph files are present in the CalcHEP files. ufo_to_mdl does not explicitly check that the vertex functions and Lorentz indices are in agreement; it solely checks the particle content of the vertices. If there are vertices missing from the CalcHEP files, 19 ufo_to_mdl generates these vertices and writes a set of corrected CalcHEP files to a new directory <MDLDir>_ufo2mdl.
In the case of four-fermion operators, ufo_to_mdl adds an additional auxiliary field to the particle content, and creates two 3-field interactions by way of this new auxiliary mediator particle, following the prescription described in Chapter 8 of the CalcHEP manual [128]. An auxiliary field has no momentum dependence and serves only to split the vertex into a form in which CalcHEP can use. The order of fields generated by ufo_to_mdl will be identical to those in the MadGraph files, i.e. a vertex i (χΓ χ χ) ψΓ ψ ψ (47) would be broken up into two vertices, i (χΓ χ χ) φ and i ψΓ ψ ψ φ , where Γ χ is a generic Dirac structure contracted with the field χ, and φ is the auxiliary field, with Lorentz indices corresponding to Γ (either scalar, vector or tensor). As a result, operators in FeynRules files should be written pairwise.
As noted above, ufo_to_mdl can also be used as a standalone tool independent of the GUM pipeline. Running ufo_to_mdl in conversion mode, with only a directory containing MadGraph files as input, python ufo_to_mdl.py <UFODir> will generate .mdl files from scratch and save them in a new directory with name <UFODir>_ufo2mdl. The version of ufo_to_mdl released with GAMBIT v2.1 does not support non-trivial colour structures and will throw an error if it is asked to generate a vertex with implicit colour structure.

Appendix C CMB energy injection
In order to provide CMB constraints from energy injection through decays and annihilation of DM, the yields dN/dE of photons, positrons and electrons produced in these processes need to be known. With GAMBIT v2.1, the existing capabilities for the calculation of photon yields (GA_Yield 20 ) were generalised and capabilities that calculate the yields of positrons (positron_Yield) and electrons (electron_Yield) were introduced. To support future analyses of charged cosmic rays, we also introduced the capabilities antiproton_Yield and antideuteron_Yield that calculate the yields of antiprotons and anti-deuterons, respectively. These capabilities are, however, not used for the CMB energy injection calculations.
Once the yields are known, they need to be passed to DarkAges via the capability energy_injection_spectrum to derive the effective efficiency function f eff (z). For maximal flexibility, we have implemented the function energy_injection_spectrum_ProcessCatalog that automatically provides the inputs for DarkAges based on the model-dependent TH_ProcessCatalog, and the yields for photons, electrons and positrons. Once these capabilities have been provided, no further input from the user is needed.
To enable CMB energy injection constraints, the user also needs to declare that the model in question can be mapped to one of the energy injection "flag" models (AnnihilatingDM_general or DecayingDM_general) and their parameters. This can be done via a friend relationship to the appropriate "flag" model. Assuming that the model under consideration contains annihilating DM particles, the user has to define a relation to AnnihilatingDM_general, and its two parameters sigmav and mass. It is important to note that the model AnnihilatingDM_general implicitly assumes that the DM particle constitutes all of DM (f χ = 1) and that it is self-conjugate. In case that the particle is not self-conjugate, the parameter sigmav needs to be rescaled by κ = 1/2. Likewise, if the DM candidate does not constitute all of DM, sigmav needs to be rescaled by f 2 χ .
To define the translation function, the user has to make sure that the definition of the model AnnihilatingDM_general is known, i.e. the following header is included: #include →"gambit/Models/models/CosmoEnergyInjection.hpp" Furthermore, the translation function and its dependen-