labelnameListing language=C++, basicstyle=, basewidth=0.53em,0.44em,
numbers=none, tabsize=2, breaklines=true, escapeinside=@@,
showstringspaces=false, numberstyle=, keywordstyle=, stringstyle=,
identifierstyle=, commentstyle=, directivestyle=, emphstyle=, frame=single,
rulecolor=, rulesepcolor=, literate=  $\sim$ 1, moredelim=*[directive] #,
moredelim=*[directive] #  language=C++, basicstyle=, basewidth=0.53em,0.44em,
numbers=none, tabsize=2, breaklines=true, escapeinside=*@@*,
showstringspaces=false, numberstyle=, keywordstyle=, stringstyle=,
identifierstyle=, commentstyle=, directivestyle=, emphstyle=, frame=single,
rulecolor=, rulesepcolor=, literate=  $\sim$ 1,
moredelim=**[is][#define]BeginLongMacroEndLongMacro  language=C++,
basicstyle=, basewidth=0.53em,0.44em, numbers=none, tabsize=2,
breaklines=true, escapeinside=@@, numberstyle=, showstringspaces=false,
keywordstyle=, stringstyle=, identifierstyle=, commentstyle=, directivestyle=,
emphstyle=, frame=single, rulecolor=, rulesepcolor=, literate=  $\sim$ 1,
moredelim=*[directive] #, moredelim=*[directive] #  language=Python,
basicstyle=, basewidth=0.53em,0.44em, numbers=none, tabsize=2,
breaklines=true, escapeinside=@@, showstringspaces=false, numberstyle=,
keywordstyle=, stringstyle=, identifierstyle=, commentstyle=, emphstyle=,
frame=single, rulecolor=, rulesepcolor=, literate =  $\sim$ 1 as  as 3
.set.set4  language=Fortran, basicstyle=, basewidth=0.53em,0.44em,
numbers=none, tabsize=2, breaklines=true, escapeinside=@@,
showstringspaces=false, numberstyle=, keywordstyle=, stringstyle=,
identifierstyle=, commentstyle=, emphstyle=, morekeywords=and, or, true,
false, frame=single, rulecolor=, rulesepcolor=, literate=  $\sim$ 1
language=bash, basicstyle=, numbers=none, tabsize=2, breaklines=true,
escapeinside=@@, frame=single, showstringspaces=false, numberstyle=,
keywordstyle=, stringstyle=, identifierstyle=, commentstyle=, emphstyle=,
frame=single, rulecolor=, rulesepcolor=, morekeywords=gambit, cmake, make,
mkdir, gum, python, wget, tar, cp, pippi, mpirun, deletekeywords=test,
literate = /gambit/gambit6 gambit/gambit/6 gum/gum/4 /include/include8
cmake/cmake/6 .cmake.cmake6 .gum.gum6 .tar.tar4 source/source/7 typetype5
$\sim$ 1 mathmath4  language=bash, basicstyle=, numbers=none, tabsize=2,
breaklines=true, escapeinside=*@@*, frame=single, showstringspaces=false,
numberstyle=, keywordstyle=, stringstyle=, identifierstyle=, commentstyle=,
emphstyle=, frame=single, rulecolor=, rulesepcolor=, morekeywords=gambit,
cmake, make, mkdir, deletekeywords=test, literate = gambit gambit7
/gambit/gambit6 gambit/gambit/6 /include/include8 cmake/cmake/6 .cmake.cmake6
$\sim$ 1  language=, basicstyle=, identifierstyle=, numbers=none, tabsize=2,
breaklines=true, escapeinside=*@@*, showstringspaces=false, frame=single,
rulecolor=, rulesepcolor=, literate=  $\sim$ 1  language=bash,
escapeinside=@@, keywords=true,false,null, otherkeywords=, keywordstyle=,
basicstyle=, identifierstyle=, sensitive=false, commentstyle=,
morecomment=[l]#, morecomment=[s]/**/, stringstyle=, moredelim=**[s][],:,
moredelim=**[l][]:, morestring=[b]’, morestring=[b]", literate = ---\---3 >>1
gtr>1 grt>1 ||1 \-  \- 3 } }1 { {1 [ [1 ] ]1  $\sim$ 1, breakindent=0pt,
breakatwhitespace, columns=fullflexible  language=bash, escapeinside=@@,
keywords=true,false,null,all, otherkeywords=, keywordstyle=, basicstyle=,
identifierstyle=, sensitive=false, commentstyle=, morecomment=[l]#,
morecomment=[s]/**/, stringstyle=, moredelim=**[l][]:, morestring=[b]’,
morestring=[b]", literate = ---\---3 grt>1 gtr>1 />>1 /<<1 lss<1 pls+1 mns-1
||1 \- - 3 } }1 { {1 [ [1 ] ]1  $\sim$ 1, breakindent=0pt, breakatwhitespace,
columns=fullflexible  language=Mathematica, basicstyle=,
basewidth=0.53em,0.44em, numbers=none, tabsize=2, breaklines=true, postbreak=,
escapeinside=@@, numberstyle=, showstringspaces=false, numberstyle=,
keywordstyle=, stringstyle=, identifierstyle=, commentstyle=, directivestyle=,
emphstyle=, frame=single, rulecolor=, rulesepcolor=, literate=  $\sim$ 1,
moredelim=*[directive] #, moredelim=*[directive] #, mathescape=false  ∎

11institutetext:  School of Physics and Astronomy, Monash University,
Melbourne, VIC 3800, Australia22institutetext:  Department of Physics and
Institute of Theoretical Physics, Nanjing Normal University, Nanjing, Jiangsu
210023, China33institutetext:  Arthur B. McDonald Canadian Astroparticle
Physics Research Institute, Kingston ON K7L 3N6, Canada44institutetext:
Department of Physics, Engineering Physics and Astronomy, Queen’s University,
Kingston ON K7L 3N6, Canada55institutetext:  Center for Cosmology, Particle
Physics and Phenomenology, Université catholique de Louvain, B-1348 Louvain-
la-Neuve, Belgium66institutetext:  Department of Physics, Imperial College
London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ,
UK77institutetext:  School of Mathematics and Physics, The University of
Queensland, St. Lucia, Brisbane, QLD 4072, Australia88institutetext:
Department of Physics, University of Oslo, N-0316 Oslo, Norway99institutetext:
Department of Physics, University of Cincinnati, Cincinnati, Ohio 45221,
USA1010institutetext:  Department of Physics, Weber State University, 1415
Edvalson St., Dept. 2508, Ogden, UT 84408, USA1111institutetext:  Bureau of
Meteorology, Melbourne, VIC 3001, Australia1212institutetext:  Institute for
Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University,
D-52056 Aachen, Germany1313institutetext:  Kavli Institute for Cosmology,
University of Cambridge, Madingley Road, Cambridge, CB3 0HA,
UK1414institutetext:  Cavendish Laboratory, University of Cambridge, JJ
Thomson Avenue, Cambridge, CB3 0HE, UK1515institutetext:  Univ Lyon, Univ Lyon
1, CNRS/IN2P3, Institut de Physique des 2 Infinis de Lyon, UMR 5822, 69622
Villeurbanne, France1616institutetext:  Theoretical Physics Department, CERN,
CH-1211 Geneva 23, Switzerland1717institutetext:  Physikalisches Institut der
Rheinischen Friedrich-Wilhelms-Universität Bonn, 53115 Bonn,
Germany1818institutetext:  Oskar Klein Centre for Cosmoparticle Physics,
AlbaNova University Centre, SE-10691 Stockholm, Sweden1919institutetext:  ARC
Centre of Excellence for Dark Matter Particle Physics & CSSM, Department of
Physics, University of Adelaide, Adelaide, SA 50052020institutetext:  Istituto
Nazionale di Fisica Nucleare, Sezione di Torino, via P. Giuria 1, I–10125
Torino, Italy2121institutetext:  Perimeter Institute for Theoretical Physics,
Waterloo ON N2L 2Y5, Canada2222institutetext:  DESY, Notkestraße 85, D-22607
Hamburg, Germany

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

The GAMBIT Collaboration: Peter Athronmonash,nanjing  Neal Avis
Kozarmcdonald,queens
Csaba Balázsmonash  Ankit Beniwallouvain,a  Sanjay Bloorimperial,uq,b  Torsten
Bringmannoslo
Joachim Brodcincy  Christopher Changuq  Jonathan M. Cornellwsu
Ben Farmerbom  Andrew Fowlienanjing  Tomás E. Gonzalomonash,aachen,c  Will
Handleykicc,cambridge  Felix Kahlhoeferaachen,d  Anders Kvellestadoslo  Farvah
Mahmoudilyon,cern
Markus T. Primbonn  Are Raklevoslo  Janina J. Renkimperial,okc  Andre
Scaffidiadelaide,infn  Pat Scottimperial,uq  Patrick Stöckeraachen  Aaron C.
Vincentmcdonald,queens,perimeter  Martin Whiteadelaide  Sebastian Wilddesy
Jure Zupancincy

(Received: date / Accepted: date)

###### Abstract

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 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 _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.

††journal: Eur. Phys. J. C

###### Contents

  1. 1 Introduction
  2. 2 Dark Matter Effective Field Theory
    1. 2.1 Running and mixing
    2. 2.2 EFT validity
    3. 2.3 Parameter ranges
  3. 3 Constraints
    1. 3.1 Direct detection
    2. 3.2 Relic abundance of DM
    3. 3.3 Indirect detection with gamma rays
    4. 3.4 Other indirect detection constraints
    5. 3.5 Collider physics
    6. 3.6 Nuisance parameter likelihoods
  4. 4 Results
    1. 4.1 Capped LHC likelihood
    2. 4.2 Full LHC likelihood
  5. 5 Conclusions and outlook
  6. A DirectDM interface
  7. B UFO to CalcHEP
  8. C CMB energy injection

##  1 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 cosmologically-observed
abundance Lee:1977ua . 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 Arcadi:2017kky , others have argued that such predictions are
premature Leane:2018kjk .

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.

In spite of these limitations, the bottom-up approach is well-advised given
the lack of direct evidence pointing to the properties of DM. The EFT approach
in particular is highly suitable for low-velocity environments such as direct
detection Fan:2010gt ; Agrawal:2010fh ; Fitzpatrick:2010br ; Crivellin:2014gpa
; DEramo:2016gos ; Hoferichter:2016nvd ; Kahlhoefer:2016eds  and indirect
detection Goodman:2010qn ; Beltran:2008xg ; Cheung:2011nt ; Harnik:2008uu ;
DeSimone:2013gj ; Karwin:2016tsw ; Carpenter:2016thc . At higher energy
scales, the EFT approach starts breaking down, such that simplified models
have become the theories of choice for the interpretation of LHC searches
Abdallah:2015ter ; Kahlhoefer:2017dnp  (see also Refs. Alanne:2017oqj ;
Alanne:2020xcb  for a hybrid approach called “Extended Dark Matter EFT”).
Nevertheless, there is an extensive literature on EFTs at colliders Bai:2010hh
; Dreiner:2013vla ; Zhou:2013fla ; Fox:2012ee ; Rajaraman:2011wf ;
Goodman:2010ku ; Fox:2011pm ; Beltran:2010ww ; Buchmueller:2013dya ;
Belyaev:2016pxe ; Pobbe:2017wrj  including studies by ATLAS ATLAS:2012ky  and
CMS Chatrchyan:2012me , which may help to shed light on the nature of DM when
interpreted with care.

A common approach to the analysis of EFTs for DM in the literature has been to
consider a single operator at a time Buckley:2011kk ; Cheung:2012gi ;
MarchRussell:2012hi ; Zheng:2010js ; Belyaev:2018pqr ; Bertuzzo:2017lwt ;
DelNobile:2013sia  and compare experimental bounds on the new physics scale
${\rm\Lambda}$ 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
Kumar:2013iva , 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.
Balazs:2014rsa , but no collider constraints were included in the analysis and
no couplings to gluons were considered. More recently, Ref. Liem:2016xpm
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. Matsumoto:2014rxa ; Blennow:2015yca ;
Matsumoto:2016hbs .

In the present work, we exploit the computational power of the GAMBIT
framework gambit  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 Duerr:2014wra  or other anomaly-free gauge
extensions that require additional stable fermions Dudas:2013sia ;
Bauer:2018egk . 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 DarkBit . 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 Fitzpatrick:2012ix  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 CosmoBit CosmoBit , 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) Gonzalo:2021cnq ; GUM  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 Shoemaker:2011vi ; Busoni:2013lha ; Busoni:2014sya ;
Busoni:2014haa ; Endo:2014mja ; Bell:2016obu ; Racco:2015dxa ;
Bruggisser:2016nzw . 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 ${\rm\Lambda}$ from the individual Wilson
coefficients $\mathcal{C}$ (rather than scanning over a combination such as
$\mathcal{C}/{\rm\Lambda}^{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
${\rm\Lambda}$ (${\rm\Lambda}<200\,\mathrm{GeV}$), where there are no relevant
LHC constraints and very large values of ${\rm\Lambda}$
(${\rm\Lambda}>1.5\,\mathrm{TeV}$), where LHC constraints are largely robust.
Of particular interest are the intermediate values of ${\rm\Lambda}$
(${\rm\Lambda}\sim 700\text{--}900\,\mathrm{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 ${\rm\Lambda}$ and the DM mass $m_{\chi}$. In
particular, even with the most general set of operators we consider, it is
impossible to simultaneously have a small DM mass ($m_{\chi}\lesssim
100\,\mathrm{GeV}$) and a large new physics scale
(${\rm\Lambda}>200\,\mathrm{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 ${\rm\Lambda}$
is more than one order of magnitude larger than $m_{\chi}$ (up to the
unitarity bound at a few hundred TeV Griest:1989wd , where the maximum
possible value of ${\rm\Lambda}$ approaches $m_{\chi}$).

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 Zenodo_DMEFT .

##  2 Dark Matter Effective Field Theory

In this study, we consider possible interactions of SM fields with a Dirac
fermion DM field, $\chi$, that is a singlet under the SM gauge group. For
phenomenological reasons discussed in detail in Sec. 3, we focus on
interactions between $\chi$ 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.
Bishara:2017nnn ; Brod:2017bsw , the interaction Lagrangian for the theory can
be written as

$\mathcal{L}_{\rm{int}}=\sum_{a,d}\dfrac{\mathcal{C}_{a}^{(d)}}{{\rm\Lambda}^{d-4}}\mathcal{Q}_{a}^{(d)}\,,$
(1)

where $\mathcal{Q}_{a}^{(d)}$ is the DM-SM operator, $d\geq 5$ is the mass
dimension of the operator, $\mathcal{C}_{a}^{(d)}$ is the dimensionless Wilson
coefficient associated to $\mathcal{Q}_{a}^{(d)}$, and ${\rm\Lambda}$ is the
scale of new physics (which can be identified with the mediator mass). The
full Lagrangian for the theory is then

$\mathcal{L}=\mathcal{L}_{\rm{SM}}+\mathcal{L}_{\rm{int}}+\overline{\chi}\left(i\not{\partial}-m_{\chi}\right)\chi\,,$
(2)

such that the free parameters of the theory are the DM mass $m_{\chi}$, the
scale of new physics ${\rm\Lambda}$, and the set of dimensionless Wilson
coefficients $\\{\mathcal{C}_{a}^{(d)}\\}$.

For sufficiently large ${\rm\Lambda}$, the phenomenology at small energies is
dominated by the operators of lowest dimension, and we therefore limit
ourselves to $d\leq 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) Brod:2017bsw . The
large set of possible operators poses a challenge for a global statistical
analysis where bounds on ${\rm\Lambda}$ and $\\{\mathcal{C}_{a}^{(d)}\\}$ 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, $\chi$, and SM
quark fields, $q$,

$\displaystyle\mathcal{Q}_{1,q}^{(6)}$
$\displaystyle=(\overline{\chi}\gamma_{\mu}\chi)(\overline{q}\gamma^{\mu}q)\,,$
(3) $\displaystyle\mathcal{Q}_{2,q}^{(6)}$
$\displaystyle=(\overline{\chi}\gamma_{\mu}\gamma_{5}\chi)(\overline{q}\gamma^{\mu}q)\,,$
(4) $\displaystyle\mathcal{Q}_{3,q}^{(6)}$
$\displaystyle=(\overline{\chi}\gamma_{\mu}\chi)(\overline{q}\gamma^{\mu}\gamma_{5}q)\,,$
(5) $\displaystyle\mathcal{Q}_{4,q}^{(6)}$
$\displaystyle=(\overline{\chi}\gamma_{\mu}\gamma_{5}\chi)(\overline{q}\gamma^{\mu}\gamma_{5}q)\,.$
(6)

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}_{\mu\nu}$ or its dual
$\widetilde{G}_{\mu\nu}=\frac{1}{2}\epsilon_{\mu\nu\rho\sigma}G^{\rho\sigma}$,
as well as operators constructed from scalar, pseudoscalar and tensor
bilinears:

$\displaystyle\mathcal{Q}_{1}^{(7)}$
$\displaystyle=\frac{\alpha_{s}}{12\pi}(\overline{\chi}\chi)G^{a\mu\nu}G^{a}_{\mu\nu}\,,$
(7) $\displaystyle\mathcal{Q}_{2}^{(7)}$
$\displaystyle=\frac{\alpha_{s}}{12\pi}(\overline{\chi}i\gamma_{5}\chi)G^{a\mu\nu}G^{a}_{\mu\nu}\,,$
(8) $\displaystyle\mathcal{Q}_{3}^{(7)}$
$\displaystyle=\frac{\alpha_{s}}{8\pi}(\overline{\chi}\chi)G^{a\mu\nu}\widetilde{G}^{a}_{\mu\nu}\,,$
(9) $\displaystyle\mathcal{Q}_{4}^{(7)}$
$\displaystyle=\frac{\alpha_{s}}{8\pi}(\overline{\chi}i\gamma_{5}\chi)G^{a\mu\nu}\widetilde{G}^{a}_{\mu\nu}\,,$
(10) $\displaystyle\mathcal{Q}_{5,q}^{(7)}$
$\displaystyle=m_{q}(\overline{\chi}\chi)(\overline{q}q)\,,$ (11)
$\displaystyle\mathcal{Q}_{6,q}^{(7)}$
$\displaystyle=m_{q}(\overline{\chi}i\gamma_{5}\chi)(\overline{q}q)\,,$ (12)
$\displaystyle\mathcal{Q}_{7,q}^{(7)}$
$\displaystyle=m_{q}(\overline{\chi}\chi)(\overline{q}i\gamma_{5}q)\,,$ (13)
$\displaystyle\mathcal{Q}_{8,q}^{(7)}$
$\displaystyle=m_{q}(\overline{\chi}i\gamma_{5}\chi)(\overline{q}i\gamma_{5}q)\,,$
(14) $\displaystyle\mathcal{Q}_{9,q}^{(7)}$
$\displaystyle=m_{q}(\overline{\chi}\sigma^{\mu\nu}\chi)(\overline{q}\sigma_{\mu\nu}q)\,,$
(15) $\displaystyle\mathcal{Q}_{10,q}^{(7)}$
$\displaystyle=m_{q}(\overline{\chi}i\sigma^{\mu\nu}\gamma_{5}\chi)(\overline{q}\sigma_{\mu\nu}q)\,.$
(16)

The definition of the operators describing interactions with the gluons,
$\mathcal{Q}_{1\text{--}4}^{(7)}$, 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,
$\mathcal{Q}_{5\text{--}10,q}^{(7)}$, 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
($\overline{\rm{MS}}$) scheme.

The complete dimension-seven basis below the EW scale contains eight
additional operators with derivatives acting on the DM fields Brod:2017bsw .
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 $\mathcal{Q}_{5\text{--}10,q}^{(7)}$ are not invariant
under EW gauge transformations, and are thus replaced in the DM EFT above the
EW scale by operators of the form $(\bar{\chi}\chi)(\bar{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 $\mathcal{Q}_{5\text{--}10,q}^{(7)}$ can also be
used in our fitting procedure above the EW scale.

In principle, analogous operators to $\mathcal{Q}_{1\text{--}10,q}^{(7)}$
exist for leptons instead of quarks Kopp:2009et ; Fox:2011fx  and weak gauge
bosons instead of gluons Weiner:2012gm ; Frandsen:2012db ; Paz:2020pbc
.111Furthermore, there are two additional dimension-6 operators describing DM-
photon interactions: the anapole moment and the charge radius Kavanagh:2018xeh
. For a recent discussion of LHC constraints on these operators we refer to
Ref. Arina:2020mxo . 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 ${\rm\Lambda}$. If this
scale is larger than the top mass, ${\rm\Lambda}>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 ${\rm\Lambda}<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
${\rm\Lambda}$.

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
$\mathcal{C}_{i,d}^{(d)}=\mathcal{C}_{i,s}^{(d)}=\mathcal{C}_{i,b}^{(d)}$ and
$\mathcal{C}_{i,u}^{(d)}=\mathcal{C}_{i,c}^{(d)}=\mathcal{C}_{i,t}^{(d)}$),
and the assumption of isospin invariance (which implies
$\mathcal{C}_{i,d}^{(d)}=\mathcal{C}_{i,u}^{(d)}$).222These constraints also
ensure that the dimension-six operators do not explicitly break electroweak
symmetry, which requires
$\mathcal{C}_{1,u}^{(6)}-\mathcal{C}_{3,u}^{(6)}=\mathcal{C}_{1,d}^{(6)}-\mathcal{C}_{3,d}^{(6)}$
Haisch:2016usn .  Hence, each operator comes with only one free parameter in
addition to the global parameters ${\rm\Lambda}$ and $m_{\chi}$. Under these
assumptions, the two EFTs above and below the EW scale have the same number of
free parameters.

###  2.1 Running and mixing

For many applications, the RG running of the Wilson coefficients (i.e. their
dependence on the energy scale $\mu$) can be neglected. In fact, the operators
$\mathcal{Q}_{1,q}^{(6)}$, $\mathcal{Q}_{2,q}^{(6)}$,
$\mathcal{Q}_{5,q}^{(7)}$ and $\mathcal{Q}_{6,q}^{(7)}$ have vanishing
anomalous dimension, while $\mathcal{Q}_{3,q}^{(6)}$,
$\mathcal{Q}_{4,q}^{(6)}$, $\mathcal{Q}_{7,q}^{(7)}$,
$\mathcal{Q}_{8,q}^{(7)}$ as well as $\mathcal{Q}_{1\text{--}4}^{(7)}$ exhibit
no running at one-loop order in QCD Hill:2014yxa . Nevertheless, there are two
cases when the effects of running can be important:

  1. 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 cross-section 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 Bishara:2018vix ; Brod:2018ust . In
our case, the dominant effects arise from the top quark Yukawa and are
discussed below.

  2. 2.

Threshold corrections: Whenever the scale $\mu$ 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 $\mathcal{Q}_{5\text{--}8,q}^{(7)}$ onto the
operators $\mathcal{Q}_{1\text{--}4}^{(7)}$ at the heavy quark thresholds,
which is given by

$\displaystyle\mathcal{C}_{1}^{(7)}$
$\displaystyle=\mathcal{C}_{1}^{(7)}-\mathcal{C}_{5,q}^{(7)}\,,$
$\displaystyle\mathcal{C}_{2}^{(7)}$
$\displaystyle=\mathcal{C}_{2}^{(7)}-\mathcal{C}_{6,q}^{(7)}\,,$
$\displaystyle\mathcal{C}_{3}^{(7)}$
$\displaystyle=\mathcal{C}_{3}^{(7)}+\mathcal{C}_{7,q}^{(7)}\,,$
$\displaystyle\mathcal{C}_{4}^{(7)}$
$\displaystyle=\mathcal{C}_{4}^{(7)}+\mathcal{C}_{8,q}^{(7)}\,.$ (17)

Mixing of the tensor operators
$\mathcal{Q}_{9,q}^{(7)},\mathcal{Q}_{10,q}^{(7)}$ above the EW scale and
subsequent matching gives rise to the dimension-five dipole operators

$\displaystyle\mathcal{Q}_{1}^{(5)}$
$\displaystyle=\frac{e}{8\pi^{2}}(\overline{\chi}\sigma_{\mu\nu}\chi)F^{\mu\nu}\,,$
(18) $\displaystyle\mathcal{Q}_{2}^{(5)}$
$\displaystyle=\frac{e}{8\pi^{2}}(\overline{\chi}i\sigma_{\mu\nu}\gamma_{5}\chi)F^{\mu\nu}\,,$
(19)

where $F_{\mu\nu}$ 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.333Note that as per our
assumptions the Wilson coefficients
$\mathcal{C}_{1}^{(5)},\mathcal{C}_{2}^{(5)}$ are taken to be zero at scale
${\rm\Lambda}$ and are only generated by the RG effects.

In the present work we include these effects as follows. To calculate the
Wilson coefficients at the hadronic scale $\mu=2\,\mathrm{GeV}$ (relevant for
direct detection) we make use of the public code DirectDM v2.2.0
Bishara:2017nnn ; Brod:2017bsw , 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
$\mu=2\,\mathrm{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\,\mathrm{GeV}$. For
${\rm\Lambda}<m_{t}$ (five-flavour EFT), we can therefore directly pass the
Wilson coefficients defined above to DirectDM. For ${\rm\Lambda}>m_{t}$ (six-
flavour EFT), there are three additional effects that are considered. First,
as pointed out in Ref. Haisch:2013uaa , the operators
$\mathcal{Q}_{9,10,t}^{(7)}$ give a contribution to the dipole operators
$\mathcal{Q}_{1,2}^{(5)}$ at the one-loop level, which is given by444For
historical reasons, in the numerical code $\log(m_{t}^{2}/{\rm\Lambda}^{2})$
instead of $\log(m_{Z}^{2}/{\rm\Lambda}^{2})$ was used. The effect on the
numerical results is negligible.

$\mathcal{C}_{1,2}^{(5)}(m_{Z})=\frac{4\,m_{t}^{2}}{{\rm\Lambda}^{2}}\log\left(\frac{m_{Z}^{2}}{{\rm\Lambda}^{2}}\right)\,\mathcal{C}_{9,10;t}^{(7)}({\rm\Lambda})\,.$
(20)

Second, as pointed out first in Ref. Crivellin:2014qxa , the operator with an
axial-vector top-quark current $\mathcal{Q}_{3,t}^{(6)}$ mixes into the
operators $\mathcal{Q}_{1,q}^{(6)}$ with light quark vector currents. The
relevant effects are given by Bishara:2018vix

$\begin{split}\mathcal{C}_{1,u/d}^{(6)}&(m_{Z})=\mathcal{C}_{1,u/d}^{(6)}({\rm\Lambda})\\\\[5.0pt]
&+\frac{2s_{w}^{2}\mp(3-6s_{w}^{2})}{8\pi^{2}}\frac{m_{t}^{2}}{v^{2}}\log\left(\frac{m_{Z}^{2}}{{\rm\Lambda}^{2}}\right)\,\mathcal{C}_{3,t}^{(6)}({\rm\Lambda})\,,\end{split}$
(21)

after integrating out the $Z$ boson at the weak scale. Here,
$s_{w}\equiv\sin\theta_{w}$ with $\theta_{w}$ the weak mixing angle, and
$v=246\,$GeV is the Higgs field vacuum expectation value. The flavour
universal UV contributions $\mathcal{C}_{1,u/d}^{(6)}({\rm\Lambda})$ 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 ${\rm\Lambda}$ and $m_{Z}$, i.e. all Wilson
coefficients other than $\mathcal{C}_{1,2}^{(5)}$ and
$\mathcal{C}_{1\text{--}4}^{(7)}$ are directly passed to DirectDM.555Small
remnant effects of the bottom and charm Yukawa coupling are taken into account
below the EW scale via double weak insertions Brod:2018ust  that are included
in the DirectDM code.

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 is a good approximation for the operators
$\mathcal{Q}_{1\text{--}4,q}^{(6)}$ and $\mathcal{Q}_{1\text{--}4}^{(7)}$. For
the operators $\mathcal{Q}_{5\text{--}10,q}^{(7)}$ mixing effects are known to
be important in principle Haisch:2012kf , 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 $(\alpha_{s})$ and use the pole quark masses
$(m_{q}^{\text{pole}})$ instead of the running quark masses. Moreover, we
neglect a small loop-level contribution from the operators
$\mathcal{Q}_{5\text{--}8,q}^{(7)}$ to the operators
$\mathcal{Q}_{1\text{--}4}^{(7)}$.

###  2.2 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.

Figure 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/{\rm\Lambda}^{2}$. The resulting exclusion is indicated
by the red shaded region. By independently varying ${\rm\Lambda}$, we can
include additional information from experiments that give weaker bounds on
$C/{\rm\Lambda}^{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).

To address this issue in the present study, we simultaneously vary the overall
scale ${\rm\Lambda}$, which corresponds to the energy where new degrees of
freedom become relevant and the EFT description breaks down, and the Wilson
coefficients $\mathcal{C}_{a}^{(d)}$ for each operator. Doing so introduces a
degeneracy, because cross sections are invariant under the rescaling
${\rm\Lambda}\to\alpha{\rm\Lambda}$ and
$\mathcal{C}_{a}^{(d)}\to\alpha^{d-4}\mathcal{C}_{a}^{(d)}$. However, the
advantage of this approach is that the parameter ${\rm\Lambda}$ 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 ${\rm\Lambda}$
and $\mathcal{C}_{a}^{(d)}$ separately to the naive approach where only
$\mathcal{C}_{a}^{(d)}/{\rm\Lambda}^{d-4}$ is constrained.

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 ${\rm\Lambda}$ 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
${\rm\Lambda}$. 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 ${\rm\Lambda}$.

For the purpose of direct detection constraints, the only requirement on
${\rm\Lambda}$ 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 ${\rm\Lambda}\gtrsim 2\,\mathrm{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 $\chi\chi\to qq$ or
$\chi\chi\to 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 ${\rm\Lambda}>2m_{\chi}$. Parameter points with smaller values of
${\rm\Lambda}$ will thus be invalidated. A dedicated study of direct detection
constraints for ${\rm\Lambda}<2m_{\chi}$ 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_{\chi\chi}<{\rm\Lambda}$ Berlin:2014cfa . 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_{\chi\chi}$ does
not directly correspond to any observable quantity (such as the missing energy
$\not{E}_{T}$ of the event) and hence the impact of varying ${\rm\Lambda}$ 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_{\chi\chi}$ 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 ${\rm\Lambda}$ 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 $\not{E}_{T}<{\rm\Lambda}$ and only apply the EFT
validity requirement for larger values of $\not{E}_{T}$. This approach is less
conservative than the one advocated, for instance in Refs. Racco:2015dxa ;
Bruggisser:2016nzw , where the energy scale of the event is taken to be the
partonic centre-of-mass energy $\sqrt{\hat{s}}$, but it has the crucial
advantage that it can be applied _after_ event generation, since the
differential cross-section with respect to missing energy
$\mathrm{d}\sigma/\mathrm{d}\not{E}_{T}$ is exactly the quantity that is
directly compared to data.666We emphasize that $m_{\chi\chi}$ and
$\not{E}_{T}$ are not strongly correlated in the sense that there are events
with both $\not{E}_{T}\ll m_{\chi\chi}$ (if the DM pair is emitted
approximately in the longitudinal direction) and $\not{E}_{T}\gg m_{\chi\chi}$
(if the two DM particles are light and approximately collinear). Since our
approach does not modify the spectrum for $\not{E}_{T}<{\rm\Lambda}$, 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
$\not{E}_{T}$, where our prescription is more appropriate.

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 $\mathrm{d}\sigma/\mathrm{d}\not{E}_{T}=0$ for $\not{E}_{T}>{\rm\Lambda}$.
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
${\rm\Lambda}$. For this we make the replacement

$\frac{\mathrm{d}\sigma}{\mathrm{d}\not{E}_{T}}\to\frac{\mathrm{d}\sigma}{\mathrm{d}\not{E}_{T}}\left(\frac{\not{E}_{T}}{{\rm\Lambda}}\right)^{-a},$
(22)

for $\not{E}_{T}>{\rm\Lambda}$. Here $a$ is a free parameter that depends on
the specific UV completion. The limits $a\to 0$ and $a\to\infty$ 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 ${\rm\Lambda}$, one finds $a\approx 2$ Buchmueller:2013dya
. 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 making both ${\rm\Lambda}$ and the Wilson
coefficients arbitrarily large. We therefore require
$\mathcal{C}_{a}^{(d)}<4\pi$, which is necessary for a perturbative UV
completion and ensures that there is no unitarity violation in the validity
range of the EFT Bell:2016obu .

One drawback of this prescription is that the EFT validity requirement depends
on the normalisation of the effective operators. For example, we have written
$\mathcal{Q}_{1\text{--}2}^{(7)}$ with a prefactor $\alpha_{s}/(12\pi)$ and
$\mathcal{Q}_{3\text{--}4}^{(7)}$ with a prefactor $\alpha_{s}/(8\pi)$ 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\pi\alpha_{s}$. With the latter
convention any constraint on the new physics scale ${\rm\Lambda}$ becomes
stronger by a factor $(48\pi^{2})^{1/3}\approx 5.3$ for
$\mathcal{Q}_{1\text{--}2}^{(7)}$ and by a factor $(32\pi^{2})^{1/3}\approx
4.6$ for $\mathcal{Q}_{3\text{--}4}^{(7)}$, meaning that much larger values of
${\rm\Lambda}$ 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 $\mathcal{Q}_{1\text{--}4}^{(7)}$.777We note that the
explicit factor of $m_{q}$ in the definition of
$\mathcal{Q}_{5\text{--}10,q}^{(7)}$ 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.

###  2.3 Parameter ranges

In this study we focus on the following parameter regions. In order to be able
to neglect QCD resonances in the process $\chi\bar{\chi}\to q\bar{q}$, we
restrict ourselves to $m_{\chi}>5\,\mathrm{GeV}$. In order to have a
sufficiently large separation of scales between the new physics scale
${\rm\Lambda}$ and the hadronic scale, we also require
${\rm\Lambda}>20\,\mathrm{GeV}$. As discussed in Sec. 2.2, we furthermore
impose the bound $|\mathcal{C}_{a}^{(d)}|<4\pi$ on all Wilson coefficients and
the bound ${\rm\Lambda}>2m_{\chi}$. The upper bounds on $m_{\chi}$ and
${\rm\Lambda}$ depend on the details of the scans that we perform and will be
discussed in Sec. 4.

##  3 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\mathcal{L}^{\text{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\mathcal{L}^{\text{max}}$. The sum of all these contributions,
$\ln\mathcal{L}^{\text{ideal}}=-105.3$ will be used to calculate log-
likelihood differences below.

Experiment | $\ln\mathcal{L}^{\text{bg}}$ | $\ln\mathcal{L}^{\text{max}}$
---|---|---
CDMSlite Agnese:2015nto  | $-16.68$ |
CRESST-II Angloher:2015ewa  | $-27.59$ |
CRESST-III Abdelhameed:2019hmk  | $-27.22$ |
DarkSide 50 Agnes:2018fwg  | $-0.09$ |
LUX 2016 LUXrun2  | $-1.47$ |
PICO-60 Amole:2017dex ; Amole:2019fdf  | $-1.496$ |
PandaX Tan:2016zwf ; Cui:2017nnn  | $-3.436$ |
XENON1T Aprile:2018dbl  | $-3.651$ |
ATLAS monojet Aad:2021egl  | $0$ |
CMS monojet Sirunyan:2017jix  | $0$ |
_Fermi_ -LAT LATdwarfP8  | $-33.245$ |
IceCube 79-string IC79_SUSY  | $0$ |
_Planck_ 2018: $p_{\text{ann}}$ Aghanim:2018eyx  | $-1.507$ |
_Planck_ 2018: $\Omega h^{2}$ Aghanim:2018eyx  |  | $5.989$
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\mathcal{L}^{\text{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\mathcal{L}^{\text{max}}$.

###  3.1 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_{\mathrm{R}}$ of recoiling nuclei. The differential event rate with respect
to recoil energy is given by

$\frac{\text{d}R}{\text{d}E_{\text{R}}}=\frac{\rho_{0}}{m_{T}\,m_{\chi}}\int_{v_{\text{min}}}^{\infty}vf(v)\frac{\text{d}\sigma}{\mathrm{d}E_{\text{R}}}\text{d}^{3}v\;,$
(23)

where $\rho_{0}$ is the local DM density, $m_{T}$ is the target nucleus mass,
$f(v)$ is the local DM velocity distribution and

$v_{\text{min}}(E_{\text{R}})=\sqrt{\frac{m_{T}E_{\text{R}}}{2\,\mu^{2}}}$
(24)

is the minimal DM velocity to cause a recoil carrying away a kinetic energy
$E_{\text{R}}$, where $\mu=m_{T}\,m_{\chi}/(m_{T}+m_{\chi})$ is the reduced
mass of the DM-nucleus system.

| SI scattering | SD scattering | Annihilations
---|---|---|---
Dimension-6 operators |  |  |
$\mathcal{Q}_{1,q}^{(6)}=(\overline{\chi}\gamma_{\mu}\chi)(\overline{q}\gamma^{\mu}q)$ | unsuppressed | — |  $s$-wave
$\mathcal{Q}_{2,q}^{(6)}=(\overline{\chi}\gamma_{\mu}\gamma_{5}\chi)(\overline{q}\gamma^{\mu}q)$ | suppressed | — |  $p$-wave
$\mathcal{Q}_{3,q}^{(6)}=(\overline{\chi}\gamma_{\mu}\chi)(\overline{q}\gamma^{\mu}\gamma_{5}q)$ | — | suppressed |  $s$-wave
$\mathcal{Q}_{4,q}^{(6)}=(\overline{\chi}\gamma_{\mu}\gamma_{5}\chi)(\overline{q}\gamma^{\mu}\gamma_{5}q)$ | — | unsuppressed |  $s$-wave $\propto m_{q}^{2}/m_{\chi}^{2}$
Dimension-7 operators |  |  |
$\mathcal{Q}_{1}^{(7)}=\dfrac{\alpha_{s}}{12\pi}(\overline{\chi}\chi)G^{a\mu\nu}G^{a}_{\mu\nu}$ | unsuppressed | — |  $p$-wave
$\mathcal{Q}_{2}^{(7)}=\dfrac{\alpha_{s}}{12\pi}(\overline{\chi}i\gamma_{5}\chi)G^{a\mu\nu}G^{a}_{\mu\nu}$ | suppressed | — |  $s$-wave
$\mathcal{Q}_{3}^{(7)}=\dfrac{\alpha_{s}}{8\pi}(\overline{\chi}\chi)G^{a\mu\nu}\widetilde{G}^{a}_{\mu\nu}$ | — | suppressed |  $p$-wave
$\mathcal{Q}_{4}^{(7)}=\dfrac{\alpha_{s}}{8\pi}(\overline{\chi}i\gamma_{5}\chi)G^{a\mu\nu}\widetilde{G}^{a}_{\mu\nu}$ | — | suppressed |  $s$-wave
$\mathcal{Q}_{5,q}^{(7)}=m_{q}(\overline{\chi}\chi)(\overline{q}q)$ | unsuppressed | — |  $p$-wave $\propto m_{q}^{2}/m_{\chi}^{2}$
$\mathcal{Q}_{6,q}^{(7)}=m_{q}(\overline{\chi}i\gamma_{5}\chi)(\overline{q}q)$ | suppressed | — |  $s$-wave $\propto m_{q}^{2}/m_{\chi}^{2}$
$\mathcal{Q}_{7,q}^{(7)}=m_{q}(\overline{\chi}\chi)(\overline{q}i\gamma_{5}q)$ | — | suppressed |  $p$-wave $\propto m_{q}^{2}/m_{\chi}^{2}$
$\mathcal{Q}_{8,q}^{(7)}=m_{q}(\overline{\chi}i\gamma_{5}\chi)(\overline{q}i\gamma_{5}q)$ | — | suppressed |  $s$-wave $\propto m_{q}^{2}/m_{\chi}^{2}$
$\mathcal{Q}_{9,q}^{(7)}=m_{q}(\overline{\chi}\sigma^{\mu\nu}\chi)(\overline{q}\sigma_{\mu\nu}q)$ | loop-induced | unsuppressed |  $s$-wave $\propto m_{q}^{2}/m_{\chi}^{2}$
$\mathcal{Q}_{10,q}^{(7)}=m_{q}(\overline{\chi}i\sigma^{\mu\nu}\gamma_{5}\chi)(\overline{q}\sigma_{\mu\nu}q)$ | loop-induced | suppressed |  $s$-wave $\propto m_{q}^{2}/m_{\chi}^{2}$
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_{q}^{2}/m_{\chi}^{2}$), the $p$-wave contribution may
dominate in the relic density calculation.

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
$\mathrm{d}\sigma/\mathrm{d}E_{\mathrm{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$. The EFT of non-relativistic interactions can
be written as

$\mathcal{L}_{\text{NR}}=\sum_{i,N}c_{i}^{N}(q^{2})\,\mathcal{O}^{N}_{i}\;,$
(25)

where the operators $\mathcal{O}^{N}_{i}$ depend only on the DM spin
$\vec{S}_{\chi}$, the nucleon spin $\vec{S}_{N}$, the momentum transfer
$\vec{q}$ and the DM-nucleon relative velocity $\vec{v}$ Fitzpatrick:2012ix ;
Anand:2013yka ; Fan:2010gt .

The non-relativistic operators can be divided into four categories according
to whether or not they depend on the nucleon spin $\vec{S}_{N}$, such that
scattering is suppressed for nuclei with vanishing spin, and whether or not
they depend on $\vec{q}$ and/or $\vec{v}$, such that scattering is suppressed
in the non-relativistic limit. Specifically, $\mathcal{O}^{N}_{1}$ leads to
spin-independent (SI) unsuppressed scattering, $\mathcal{O}^{N}_{4}$ leads to
spin-dependent (SD) unsuppressed scattering, $\mathcal{O}^{N}_{5}$,
$\mathcal{O}^{N}_{8}$, $\mathcal{O}^{N}_{11}$ lead to SI momentum-suppressed
scattering and $\mathcal{O}^{N}_{6}$, $\mathcal{O}^{N}_{7}$,
$\mathcal{O}^{N}_{9}$, $\mathcal{O}^{N}_{10}$, $\mathcal{O}^{N}_{12}$ 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_{i}^{N}(q^{2})$ can be directly calculated from the Wilson
coefficients of the relativistic operators at $\mu=2\,\mathrm{GeV}$. The
explicit dependence on the momentum transfer $q=\sqrt{2m_{T}E_{\text{R}}}$ is
a result of two effects. First, under RG evolution some of the effective DM-
quark operators mix into the DM dipole operators $\mathcal{Q}_{1,2}^{(5)}$
(see Eq. (20)). These operators then induce long-range interactions, i.e.
contributions to the $c_{i}^{N}(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 $\langle
N^{\prime}|\overline{q}\gamma^{\mu}q|N\rangle$. 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\to m_{\pi}$ Bishara:2016hek ;
Bishara:2017pfq .

Parameter | Value |  |  | Parameter | Value
---|---|---|---|---|---
$\sigma_{\pi N}$ |  $50(15)$ MeV Bishara:2017pfq  |  |  | $\mu_{p}$ |  $2.793$ PDG20
$Bc_{5}(m_{d}-m_{u})$ |  $-0.51(8)$ MeV Crivellin:2013ipa  |  |  | $\mu_{n}$ |  $-1.913$ PDG20
$g_{A}$ |  $1.2756(13)$ PDG20  |  |  | $\mu_{s}$ |  $-0.036(21)$ Djukanovic:2019jtp ; Sufian:2016pex
$m_{G}$ |  $836(17){\rm~{}MeV}$ Bishara:2017nnn  |  |  | $g_{T}^{u}$ |  $0.784(30)$ Gupta:2018lvp
$\sigma_{s}$ |  $52.9(7.0)\rm{~{}MeV}$ Aoki:2019cca  |  |  | $g_{T}^{d}$ |  $-0.204(15)$ Gupta:2018lvp
$\Delta u+\Delta d$ |  $0.440(44)\,$ Liang:2018pis  |  |  | $g_{T}^{s}$ |  $-27(16)\cdot 10^{-3}$ Gupta:2018lvp
$\Delta s$ |  $-0.035(9)$ Liang:2018pis  |  |  | $B_{T,10}^{u/p}$ |  $3.0(1.5)$ Pasquini:2005dk
$B_{0}m_{u}$ |  $0.0058(5)\rm{~{}GeV}^{2}$ Bishara:2017nnn  |  |  | $B_{T,10}^{d/p}$ |  $0.24(12)$ Pasquini:2005dk
$B_{0}m_{d}$ |  $0.0124(5)\rm{~{}GeV}^{2}$ Bishara:2017nnn  |  |  | $B_{T,10}^{s/p}$ |  $0.0(2)$ Pasquini:2005dk
$B_{0}m_{s}$ |  $0.249(9)\rm{~{}GeV}^{2}$ Bishara:2017nnn  |  |  | $r_{s}^{2}$ |  $-0.115(35){\rm~{}GeV}^{-2}$ Sufian:2016pex ; Djukanovic:2019jtp
Table 3: The hadronic input parameters as used in DirectDM v2.2.0.

All of these effects are fully taken into account in DirectDM, which
calculates the coefficients $c_{i}^{N}(q^{2})$ for given Wilson coefficients
$\mathcal{C}_{a}^{(d)}$ at a higher scale (see App. A). These coefficients are
then passed onto DDCalc v2.2.0 DarkBit ; HP , which calculates the
differential cross-section for each operator $\mathcal{O}^{N}_{i}$ (including
interference) 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:

$N_{p}=M\,T_{\text{exp}}\int\phi(E_{\mathrm{R}})\,\frac{\mathrm{d}R}{\mathrm{d}E_{\mathrm{R}}}\,\mathrm{d}E_{\mathrm{R}}\,,$
(26)

where $M$ is the detector mass, $T_{\text{exp}}$ is the exposure time and
$\phi(E_{\mathrm{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 Aprile:2018dbl , LUX 2016
LUXrun2 , PandaX 2016 Tan:2016zwf  and 2017 Cui:2017nnn  analyses, CDMSlite
Agnese:2015nto , CRESST-II Angloher:2015ewa  and CRESST-III
Abdelhameed:2019hmk , PICO-60 2017 Amole:2017dex  and 2019 Amole:2019fdf , and
DarkSide-50 Agnes:2018fwg .

The hadronic inputs to DirectDM v2.2.0 Bishara:2017nnn  were updated with the
most recent $N_{f}=2+1$ lattice QCD results, following the FLAG quality
requirements Aoki:2019cca , see Table 3. All the inputs are evaluated at
$\mu=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_{s}^{2}=-0.0045(14)\,$fm2 Sufian:2016pex ;
Djukanovic:2019jtp . For the nuclear magnetic moment induced by the strange
quark, $\mu_{s}=-0.036(21)$ Djukanovic:2019jtp ; Sufian:2016pex , we inflate
the errors according to the Particle Data Group prescription.

The scalar form factors at zero recoil are obtained from expressions in Ref.
Crivellin:2013ipa , namely

$\sigma_{u(d)}^{N}=\frac{\sigma_{\pi N}}{2}(1\mp\xi)\pm
Bc_{5}\,(m_{d}-m_{u})\left(1\mp\frac{1}{\xi}\right),$ (27)

where the upper (lower) sign is for the proton (neutron). We use a rather
conservative estimate $\sigma_{\pi N}=(50\pm 15)$ MeV Bishara:2017pfq  that
covers the spread between the lattice QCD Horsley:2011wr ; Durr:2015dna ;
Yang:2015uis ; Abdel-Rehim:2016won ; Bali:2016lvx ; Alexandrou:2019brg ;
Yamanaka:2018uud ; Borsanyi:2020bpd  and pionic atom determinations
Alarcon:2011zs ; Hoferichter:2015dsa ; Dmitrasinovic:2016hup ;
RuizdeElvira:2017stg ; Friedman:2019zhc ; Yang:2015uis ; Abdel-Rehim:2016won ;
Yamanaka:2018uud . The other two parameters are
$\xi\equiv(m_{d}-m_{u})/(m_{d}+m_{u})=0.36\,\pm 0.04$ and
$Bc_{5}\,(m_{d}-m_{u})=(-0.51\pm 0.08)$ MeV Crivellin:2013ipa .

The matrix elements of tensor currents are described by three sets of form
factors, but only two, $g_{T}^{q}$ and $B_{T,10}^{q/N}(0)$, enter the chirally
leading expressions. For $g_{T}^{q}$, the only $N_{f}=2+1$ result from Ref.
Yamanaka:2018uud  does not satisfy the FLAG quality requirements, so we use
the $N_{f}=2+1+1$ results from Ref. Gupta:2018lvp  instead; the difference
between the $N_{f}=2+1$ and $N_{f}=2+1+1$ results is expected to be small. For
$B_{T,10}^{q/N}(0)$, we use the results from the constituent quark model in
Ref. Pasquini:2005dk .

###  3.2 Relic abundance of DM

The Early Universe time evolution of the number density of the $\chi$
particles, $n_{\chi}$, is governed by the Boltzmann equation Gondolo:1990dk

$\displaystyle\frac{dn_{\chi}}{dt}+3Hn_{\chi}=-\langle\sigma
v_{\textrm{rel}}\rangle\left(n_{\chi}n_{\bar{\chi}}-n_{\chi,\textrm{eq}}n_{\bar{\chi},\textrm{eq}}\right)\,,$
(28)

where $n_{\chi,\textrm{eq}}$ is the number density in equilibrium, $H(t)$ is
the Hubble rate and $\langle\sigma v_{\textrm{{rel}}}\rangle$ is the thermally
averaged cross-section times the relative (Møller) velocity, given by

$\displaystyle\langle\sigma
v_{\textrm{rel}}\rangle=\int^{\infty}_{4m_{\chi}^{2}}\\!ds\,\frac{\sqrt{s-4m_{\chi}^{2}}(s-2m_{\chi}^{2})\,K_{1}\left(\sqrt{s}/T\right)}{8m_{\chi}^{4}\,TK_{2}^{2}\left(m_{\chi}/T\right)}\,\sigma
v_{\rm lab}\,,$ (29)

where $K_{1,2}$ are the modified Bessel functions and $v_{\rm 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. Binder:2017rgn ). 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_{\chi}+n_{\bar{\chi}}=2n_{\chi}$
(disregarding the possibility of an initial asymmetry Kaplan:2009ag ).

We compute tree-level annihilation cross-sections using CalcHEP v3.6.27
Pukhov:2004ca ; Belyaev:2012qa , where the implementation of the four-fermion
interactions is generated by GUM Gonzalo:2021cnq ; GUM  from UFO files via the
tool ufo_to_mdl (described in App. B). To ensure the EFT picture is valid, we
invalidate points where ${\rm\Lambda}\leq 2m_{\chi}$. We obtain the relic
density of $\chi$ by numerically solving Eq. (28) at each parameter point,
assuming the standard cosmological history888For a recent review on the
effects of non-standard cosmological scenarios, see Ref. Arbey:2021gdg . and
using the routines implemented in DarkSUSY v6.2.2 darksusy ; darksusy4  via
DarkBit. We then compare the prediction to the relic density constraint from
_Planck_ 2018: $\Omega_{\textrm{DM}}\,h^{2}=0.120\pm 0.001$ Aghanim:2018eyx .
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. gambit ; DarkBit .

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 $\langle\sigma
v_{\textrm{rel}}\rangle$ during freeze-out. The resulting effects, such as
resonances or coannihilations could both increase and decrease the resulting
value of $\Omega_{\chi}$ (see e.g. Ref. Bell:2013wua ; Baker:2015qna ), so the
relic density constraint should be interpreted with care for ${\rm\Lambda}\sim
2m_{\chi}$, 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 ${\rm\Lambda}$. Given the potential for a rich dark
sector containing $\chi$, and in particular, the possibility of additional DM
candidates not captured by the EFT, we will by default _not_ demand that the
particle $\chi$ 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

$f_{\chi}\equiv(\Omega_{\chi}+\Omega_{\bar{\chi}})/\Omega_{\mathrm{DM}}=2\Omega_{\chi}/\Omega_{\mathrm{DM}}$
(30)

and $f_{\chi}^{2}$, respectively. In doing so, we assume that the fraction
$f_{\chi}$ 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_{\chi}\approx 1$) rather than imposing the relic density as an upper bound
($f_{\chi}\leq 1$).999Note that, since we include uncertainties in both the
relic density calculation and the _Planck_ measurement, $\Omega_{\chi}h^{2}$
can deviate slightly from 0.120 even when we require that the DM relic
abundance is saturated. In this case we set
$f_{\chi}=\text{min}(\Omega_{\chi}h^{2}/0.120,1)$, which can therefore
slightly deviate from (but never exceed) unity.

###  3.3 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
${\rm\Lambda}>2m_{\chi}$.

Gamma rays from dwarf spheroidal galaxies (dSphs) are a particularly robust
way of constraining annihilation signals from DM Bringmann:2012ez . In
general, for a given energy bin $i$, the DM-induced $\gamma$-ray flux from
target $k$ can be written in the factorised form $\Phi_{i}\cdot J_{k}$, where
details of the particle physics processes are encoded in $\Phi_{i}$, and
details of the astrophysics are encoded in $J_{k}$. See the DarkBit manual
DarkBit  for more details.

In general, only operators that lead to $s$-wave annihilation
($\mathcal{Q}_{1,3,4,q}^{(6)},\mathcal{Q}_{2,4,6,8,9,10,q}^{(7)}$) give rise
to observable gamma-ray signals; see for instance, Table 2. For the operators
$\mathcal{Q}_{2,q}^{(6)}$ and $\mathcal{Q}_{1,3,5,7,q}^{(7)}$, the leading
contribution to the annihilation cross-section is $p$-wave suppressed, i.e.
proportional to $v_{\mathrm{rel}}^{2}$. As DM in dSphs is extremely cold, with
$\langle v^{2}\rangle^{1/2}\sim 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

$\displaystyle\Phi_{i}$
$\displaystyle=\frac{f_{\chi}^{2}}{4}\sum_{j}\frac{(\sigma v)_{0,j}}{4\pi
m_{\chi}^{2}}\int_{\Delta E_{i}}dE\,\frac{dN_{\gamma,j}}{dE}\,,$ (31)

where $f_{\chi}$ is the DM fraction defined in eq. (30), $(\sigma v)_{0,j}$
denotes the zero-velocity limit of the cross-section for $\chi\bar{\chi}\to j$
and $N_{\gamma,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_{\chi}=n_{\bar{\chi}}$).
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_{\gamma,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-of-sight integral over the
DM distribution assuming an NFW density profile and the solid angle $\Omega$,

$\displaystyle J_{k}$
$\displaystyle=\int_{\Delta\Omega_{k}}d\Omega\int_{\mathrm{l.o.s.}}ds\,\rho_{\rm
DM}^{2}\simeq D_{k}^{-2}\int d^{3}x\,\rho_{\rm DM}^{2}\,,$ (32)

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 LATdwarfP8 .
We use the gamLike v1.0.1 interface within DarkBit DarkBit  to compute the
likelihood for the gamma-ray observations, $\ln\mathcal{L}_{\rm{exp}}$,
constructed from the product $\Phi_{i}\cdot J_{k}$ and summed over all targets
and energy bins,

$\displaystyle\ln\mathcal{L}_{\rm{exp}}=\sum^{\rm{N_{dSphs}}}_{k=1}\sum^{\rm{N_{eBins}}}_{i=1}\ln\mathcal{L}_{ki}\left(\Phi_{i}\cdot
J_{k}\right)\,.$ (33)

We also include a contribution from profiling over the $J$-factors of each
dSph, $\ln\mathcal{L}_{J}=\sum_{k}\ln\mathcal{L}(J_{k})$ DarkBit ; LATdwarfP8
, such that the full likelihood reads

$\displaystyle\ln\mathcal{L}_{\rm{dSphs}}^{\rm{prof.}}=\underset{\\{J_{k}\\}}{\textrm{max}}\left(\ln\mathcal{L}_{\rm{exp}}+\ln\mathcal{L}_{J}\right)\,.$
(34)

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 TheFermi-LAT:2017vmf . 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
Acharyya:2020sbj . We will not include the projected CTA likelihoods in our
scans, but indicate the reach of CTA when discussing our results.

###  3.4 Other indirect detection constraints

#### 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 $\simeq 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 $\sim$ MeV scales of fusion processes. Leading
constraints have been obtained by Super-Kamiokande down to a few GeV
Choi:2015ara , and by the IceCube South Pole Neutrino Observatory, between 20
and 104 GeV Aartsen:2016zhm . 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 non-relativistic
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 Kozar:2021iur , which computes the DM capture rate
via the integral over the solar radius $r$ and DM halo velocity $u$:

$C_{\odot}(t)=4\pi\int_{0}^{R_{\odot}}r^{2}\int_{0}^{\infty}\frac{f(u)}{u}\,w\Omega(w,r)\,du\,dr\,,$
(35)

where $w(r)=\sqrt{u^{2}+v^{2}_{\text{esc},\odot}(r)}$ is the DM velocity at
position $r$, and

$\displaystyle\Omega(w)$
$\displaystyle=w\sum_{i}n_{i}(r,t)\frac{\mu_{i,+}^{2}}{\mu_{i}}\Theta\left(\frac{\mu_{i}v^{2}}{\mu_{i,-}^{2}}-u^{2}\right)$
$\displaystyle\hskip
14.22636pt\times\int_{m_{\chi}u^{2}/2}^{m_{\chi}w^{2}\mu_{i}/2\mu_{i,+}^{2}}\frac{d\sigma_{i}}{dE_{\text{R}}}\,dE_{\text{R}}\,,$
(36)

is the probability of scattering from velocity $w$ to a velocity less than the
local Solar escape velocity $v_{\mathrm{esc},\odot}(r)$,
${d\sigma_{i}}/{dE_{\text{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 $\mu_{i}=m_{\chi}/m_{N,i}$. Version 2.1 of Capt’n General uses the method
described in detail in Ref. Catena:2015uha , 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. Catena:2015uha  for the 16 most abundant
elements in the Sun. Solar parameters are based on the Barcelona Group’s
AGSS09ph Standard Solar Model Vinyoles:2016djt ; AGSS .

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,
cross-sections 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 IC79 ; IC22Methods , which computes event-level likelihoods
based on a re-analysis of the 79-string IceCube search for DM annihilation in
the Sun IC79_SUSY .

#### Cosmic Microwave Background

Additional constraints on the DM annihilation cross-section 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
$\tau$. 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_{\text{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 Slatyer15a
. The CMB is then sensitive to the following parameter combination:

$p_{\text{ann}}\equiv f_{\chi}^{2}f_{\text{eff}}\frac{\langle\sigma
v_{\mathrm{rel}}\rangle}{m_{\chi}}\;,$ (37)

where $\langle\sigma v_{\mathrm{rel}}\rangle\approx(\sigma 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_{\text{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 Slatyer:2015kla . The
first part of this calculation has been automated within DarkSUSY and is
accessible via DarkBit. The second part relies on DarkAges Stocker:2018avm
(which is part of the ExoCLASS branch of CLASS) and is accessible via CosmoBit
CosmoBit , see App. C for further details.101010It is noteworthy that DarkAges
calculates $f_{\text{eff}}(z)$ as a redshift-dependent function instead of a
single redshift-independent coefficient $f_{\rm eff}$, as it is implicitly
assumed in Eq. (37). In order to compress the function $f_{\text{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 Slatyer15a .

As the _Planck_ collaboration only quotes the 95% credible interval for
$p_{\text{ann}}$ Aghanim:2018eyx , the remaining challenge is to obtain a
likelihood for $p_{\text{ann}}$ from cosmological data. Although this
likelihood can, in principle, be calculated for each parameter point
individually using the CosmoBit 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
${\rm\Lambda}$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_{\text{ann}}$ (while profiling over the ${\rm\Lambda}$CDM and
cosmological nuisance parameters). This approach yields

$\mathcal{L}(p_{\text{ann}})=\mathcal{L}_{0}\exp\left[-\left(\frac{p_{\text{ann}}^{28}+0.48}{2.45}\right)^{2}\right]\;,$
(38)

where $p_{\text{ann}}^{28}\equiv
p_{\text{ann}}/\left(10^{-28}\,\mathrm{cm^{3}\,s^{-1}\,GeV^{-1}}\right)$. In
arriving at this result, we have included the _Planck_ TT,TE,EE+lowE+lensing
likelihoods (using the ‘lite’ likelihood for multipoles $\ell\geq 30$, which
only require one additional nuisance parameter Aghanim:2019ame ), as well as
the BAO data of 6dF 2011MNRAS.416.3017B , SDSS DR7 MGS 2015MNRAS.449..835R ,
and the SDSS BOSS DR12 galaxy sample Alam:2016hwk . This profile likelihood,
which reproduces the 95% credible interval obtained by the _Planck_
collaboration Aghanim:2018eyx , can then be used in all subsequent scans, so
that only $p_{\text{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
Kopp:2013eka .111111We note that this conclusion would be radically different
for unsuppressed direct annihilation to leptons Bergstrom:2013jra ;
Ibarra:2013zia , which would result from leptonic operators analogous to
$\mathcal{Q}_{1,3}^{(6)}$.  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 Bergstrom:1999jc ;
Bringmann:2006im ; Cuoco:2016eej , 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 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 Heisig:2020nse
and possible correlated systematics Boudaud:2019efq . 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 Johannesson:2016rlh  and hence beyond the scope
of this work.

###  3.5 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\rightarrow\chi\chi j$ as a single jet associated with missing
transverse energy ($\not{E}_{T}$). In this study, we include the CMS
Sirunyan:2017jix  and ATLAS Aad:2021egl  monojet analyses based on
$36\,\rm{fb}^{-1}$ and $139\,\rm{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 Bauer:2017fsw ;
Zhou:2013fla ; Brennan:2016xjh .

The expected number of events in a given bin of the $\not{E}_{T}$ distribution
is

$N=L\times\sigma\times(\epsilon A)\;,$

where $L=36\,\text{fb}^{-1}$ or $139\,\text{fb}^{-1}$ is the total integrated
luminosity, $\sigma$ the total production cross-section and the factor
$(\epsilon A)$ is the efficiency times acceptance for passing the kinematic
selection requirements for the analysis. Both $\sigma$ and $(\epsilon 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 Alloul:2013bka , we
have therefore produced separate interpolations of $\sigma$ and $\epsilon A$
based on the output of Monte Carlo simulations with MadGraph_aMC@NLO v2.6.6
Alwall:2011uj  (v2.9.2) for the CMS (ATLAS) analysis, interfaced to Pythia
v8.1 Sjostrand:2007gs  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 DELPHES3 . The
ColliderBit code extension that enables $\sigma$ and $(\epsilon 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 ($\mathcal{C}_{i}^{(6)}$
and $\mathcal{C}_{i=1,...,4}^{(7)}$) 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
$\mathcal{Q}_{1,q}^{(6)}$/$\mathcal{Q}_{4,q}^{(6)}$ and
$\mathcal{Q}_{2,q}^{(6)}$/$\mathcal{Q}_{3,q}^{(6)}$. For these Wilson
coefficients, we parametrize the tabulated grids in terms of a mixing angle
$\theta$, defined via $\mathcal{C}_{1,2}^{(6)}=\sin\theta$ and
$\mathcal{C}_{3,4}^{(6)}=\cos\theta$.

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
$\not{E}_{T}>1000\,\mathrm{GeV}$.121212We note that this combination also
reduces the impact of a local $\sim 2.5\sigma$ excess in the third-highest
bin, which would otherwise strongly bias our analysis. 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
Collaboration:2242860 , 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:

$\begin{split}\mathcal{L}_{\text{CMS}}(\bm{s},\bm{\gamma})=&\prod_{i=1}^{22}\left[\frac{(s_{i}+b_{i}+\gamma_{i})^{n_{i}}\,e^{-(s_{i}+b_{i}+\gamma_{i})}}{n_{i}!}\right]\\\
&\hphantom{\int}\times\frac{1}{\sqrt{\det
2\pi\Sigma}}e^{-\frac{1}{2}\bm{\gamma}^{T}\bm{\Sigma^{-1}}\bm{\gamma}}\,.\end{split}$
(39)

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 is given by $\gamma_{i}$. The correlations between the different
$\gamma_{i}$ are encoded in the covariance matrix $\bm{\Sigma}$ provided by
CMS, where we also add the signal yield uncertainties in quadrature along the
diagonal. We follow the procedure in Ref. Collaboration:2242860  in treating
the $\gamma_{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 $\bm{\gamma}$ to obtain a
likelihood solely in terms of the set of DM EFT signal estimates $\bm{s}$:

$\mathcal{L}_{\text{CMS}}(\bm{s})\equiv\mathcal{L}_{\text{CMS}}(\bm{s},\hat{\hat{\bm{\gamma}}}).$
(40)

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

$\mathcal{L}_{\text{ATLAS}}(s_{i})\equiv\mathcal{L}_{\text{ATLAS}}(s_{i},\hat{\hat{\gamma_{i}}})\,,$
(41)

where $\mathcal{L}_{\text{ATLAS}}(s_{i},\hat{\hat{\gamma_{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\mathcal{L}_{\text{LHC}}=\ln\mathcal{L}_{\text{CMS}}+\ln\mathcal{L}_{\text{ATLAS}}$.
However, due to the per-point signal region selection required in the
evaluation of $\ln\mathcal{L}_{\text{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

$\Delta\ln\mathcal{L}_{\text{LHC}}=\ln\mathcal{L}_{\text{LHC}}(\bf{s})-\ln\mathcal{L}_{\text{LHC}}(\bf{s}=\bf{0})$
(42)

as the LHC log-likelihood contribution in the parameter scan ColliderBit .

Nuisance parameter |  |  Value ($\pm$ $3\sigma$ range)
---|---|---
Local DM density | $\rho_{0}$ |  $0.2$–$0.8$ GeV cm-3
Most probable speed | $v_{\rm{peak}}$ |  $240\,(24)$ km s-1
Galactic escape speed | $v_{\textrm{esc}}$ |  $528\,(75)$ km s-1
Running top mass ($\overline{\rm{MS}}$ scheme) | $m_{t}(m_{t})$ |  $162.9\,(6.0)$ GeV
Pion-nucleon sigma term | $\sigma_{\pi N}$ |  $50\,(45)$ MeV
Strange quark contrib. to nucleon spin | $\Delta s$ | $-0.035\,(0.027)$
Strange quark nuclear tensor charge | $g_{T}^{s}$ | $-0.027\,(0.048)$
Strange quark charge radius of the proton | $r_{s}^{2}$ | $-0.115\,(0.105){\rm~{}GeV}^{-2}$
Table 4: A list of nuisance parameters that are varied simultaneously with the
DM EFT model parameters in our scans (the hadronic parameters are given at
$\mu=2$ GeV). All parameters are scanned over their $3\sigma$ range using flat
parametrisation. For more details on the respective nuisance likelihoods, see
Sec. 3.6.

When presenting the results of a global fit we identify the maximum-likelihood
point $\bf{\Theta}_{\text{best-fit}}$ in the DM EFT parameter space and map
out the $1\sigma$ and $2\sigma$ confidence regions defined using the
likelihood ratio
$\mathcal{L}(\bf{\Theta})/\mathcal{L}(\bf{\Theta}_{\text{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 $\Delta\ln\mathcal{L}_{\text{LHC}}$ with a capped version,

$\displaystyle\Delta\ln\mathcal{L}_{\text{LHC}}^{\text{cap}}(\bf{s})=\min[\mathit{\Delta}\ln\mathcal{L}_{\text{LHC}}(\bf{s}),\mathit{\Delta}\ln\mathcal{L}_{\text{LHC}}(\bf{s}=\bf{0})].$
(43)

This will assign the same log-likelihood value,
$\Delta\ln\mathcal{L}_{\text{LHC}}^{\text{cap}}=0$, for all DM EFT parameter
points whose prediction fit the collider data as well as, or better than, the
SM prediction ($\bf{s}=\bf{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 $\Delta\ln\mathcal{L}_{\text{LHC}}^{\text{cap}}$ 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.131313A practical benefit of having a continuous likelihood penalty
rather than a hard cut is that it helps guide the parameter sampler towards
the viable regions in the high-dimensional DM EFT parameter space. A more
detailed introduction to the capped likelihood construction can be found in
Ref. EWMSSM .

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.

###  3.6 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\sigma$ 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
$\rho_{0}$, where the likelihood is given by a log-normal distribution with
central value $\rho_{0}=0.40$ GeV cm-3 and error $\sigma_{\rho_{0}}=0.15$ GeV
cm-3. We scan over an asymmetric range in $\rho_{0}$ to reflect the log-normal
distribution – see Ref. DarkBit  for more details.

We follow the same treatment of the Milky Way halo as in the GAMBIT Higgs
portal study HP . We utilise Gaussian likelihoods for parameters describing
the Maxwell-Boltzmann velocity distribution, specifically the peak of the
distribution $v_{\rm{peak}}=240\,\pm\,8$ km s-1 Reid:2014boa , and the
galactic escape velocity $v_{\rm{esc}}=528\pm 25$ km s-1, based on the _Gaia_
data Deason:2019kgj .

We employ a Gaussian likelihood for the running top quark mass in the
$\overline{\textrm{MS}}$ scheme with a central value $m_{t}(m_{t})=162.9$ GeV
and an error $2.0$ GeV Aad:2019mkw .141414This is based on taking an average
of the asymmetric uncertainty $m_{t}(m_{t})=162.9^{+2.3}_{-1.6}$ GeV; see
table 2 in Ref. Aad:2019mkw . The top pole mass $(m_{t}^{\textrm{pole}})$ is
then computed using the following formula:

$m_{t}^{\textrm{pole}}=m_{t}(m_{t})\left[1+\dfrac{4}{3}\dfrac{\alpha_{s}(m_{Z})}{\pi}\right]\,.$
(44)

We use only the one-loop QCD corrections in this shift in order to be
consistent with the procedure carried out in Ref. Aad:2019mkw . We have
checked that the above expression gives the expected result for the top pole
mass and matches well with Ref. Aad:2019mkw .

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 $\sigma_{\pi N}$, $\Delta s$, $g_{T}^{s}$ and
$r_{s}^{2}$ quoted in Table 3. The remaining hadronic input parameters are
fixed to the central values given in Table 3.

##  4 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 ScannerBit  using a population of $5\times
10^{4}$ and a convergence threshold of either $10^{-5}$ or $3\times 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 high-likelihood 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\,\mathrm{GeV}<m_{\chi}<150\,\mathrm{TeV}$ and
$20\,\mathrm{GeV}<{\rm\Lambda}<300\,\mathrm{TeV}$).151515We note that for the
largest values of $m_{\chi}$ and ${\rm\Lambda}$ 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
$\mu={\rm\Lambda}$ and $\mu=m_{Z}$ becomes questionable. The constraints that
we obtain above the TeV scale are therefore only approximate and should be
interpreted with care. In these scans, $m_{\chi}$ and ${\rm\Lambda}$ 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 ($m_{\chi}<500\,\mathrm{GeV}$ and
${\rm\Lambda}<2\,\mathrm{TeV}$). In these scans the DM mass is scanned on a
linear scale, the scale ${\rm\Lambda}$ on a logarithmic scale and the Wilson
coefficients on a scale that is logarithmic on $[-4\pi,-10^{-6}]$, linear on
$[-10^{-6},10^{-6}]$ and logarithmic on $[10^{-6},4\pi]$. 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
${\rm\Lambda}$ 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.

Figure 2: Profile likelihood in the $m_{\chi}$–${\rm\Lambda}$ plane when
considering only dimension-6 operators and capping the LHC likelihood at the
value of the background-only hypothesis. The white contours indicate the
$1\sigma$ and $2\sigma$ confidence regions and the best-fit point is indicated
by the white star. The shaded region (corresponding to ${\rm\Lambda}\leq
2m_{\chi}$) 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 ${\rm\Lambda}$ and
$\mathcal{C}_{3,4}^{(6)}$ and hence the likelihood is essentially constant
across the entire yellow region (see also Fig. 4).

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_{\chi}\leq 1$) and scans where we require that the DM relic density be
saturated ($f_{\chi}\approx 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.

###  4.1 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 ${\rm\Lambda}$. 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_{\chi}$ and
${\rm\Lambda}$, the allowed parameter space is determined by the EFT validity
requirement ${\rm\Lambda}>2m_{\chi}$ and the relic density requirement which,
combined with the perturbativity bound on the Wilson coefficients, implies an
upper bound on ${\rm\Lambda}$ for given $m_{\chi}$. These different
constraints are compatible only for $m_{\chi}<150\,\mathrm{TeV}$, implying an
upper bound on the scale of new physics of ${\rm\Lambda}<300\,\mathrm{TeV}$.
This limit corresponds to the well-known unitarity bound for thermal freeze-
out Griest:1989wd .

The zoomed-in version in the right panel reveals a number of additional
features. In the top-left corner (small $m_{\chi}$, large ${\rm\Lambda}$),
there are strong constraints from the LHC, which make it impossible to satisfy
the relic density requirement. These constraints become weaker as
${\rm\Lambda}$ decreases and the EFT can only be trusted for smaller values of
$\not{E}_{T}$. The various sharp features correspond to the points where
${\rm\Lambda}$ crosses the boundary of a specific $\not{E}_{T}$ bin, leading
to a jump in the likelihood. In our conservative approach, LHC constraints are
completely absent for ${\rm\Lambda}<200\,\mathrm{GeV}$. Finally, we find that
there is a slight upward fluctuation in _Fermi_ -LAT data, which can be fitted
for $m_{\chi}=5.0\,\mathrm{GeV}$ and $f_{\chi}^{2}\langle\sigma
v\rangle_{0}=1.1\times 10^{-27}\,\mathrm{cm^{3}\,s^{-1}}$.161616We emphasize
that, although the best-fit point lies close to the boundary of the parameter
space, there is no preference for even smaller values of the DM mass and hence
our findings would not change when extending the scan range.

We emphasize that a great advantage of our approach is that we treat the new-
physics scale ${\rm\Lambda}$ 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 ${\rm\Lambda}$ is sufficiently
greater than the largest $p_{T}$ bin considered in the LHC analyses, i.e.
${\rm\Lambda}>1.3\,\mathrm{TeV}$. Conversely, for
${\rm\Lambda}<200\,\mathrm{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 ${\rm\Lambda}$, 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.

Figure 3: Profile likelihood in terms of the DM mass, the relic density and
the rescaled annihilation cross-section. As in 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 $b\bar{b}$ final state
Acharyya:2020sbj .

Figure 4: Profile likelihood in the ${\rm\Lambda}$–$\mathcal{C}^{(6)}_{4}$
plane (left) and the ${\rm\Lambda}$–$\mathcal{C}^{(6)}_{3}$ 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\sigma$ and $2\sigma$ 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
${\rm\Lambda}$ and $\mathcal{C}_{3,4}^{(6)}$ and hence the likelihood is
essentially constant across the entire yellow region.

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_{\chi}\lesssim 100\,\mathrm{GeV}$ it is impossible
to saturate the observed DM relic density, $\Omega_{\text{DM}}h^{2}=0.12$, due
to the combined constraints from direct and indirect detection experiments.
However, these constraints are suppressed for DM sub-components, such that it
is possible to have very small relic densities in this mass region. For
$m_{\chi}>100\,\mathrm{GeV}$ (corresponding to
${\rm\Lambda}>200\,\mathrm{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_{\chi}\gtrsim 1\,\mathrm{TeV}$ as
the LHC energy becomes insufficient to produce a pair of DM particles.

For $m_{\chi}\lesssim 1\,\mathrm{TeV}$, we find that there is a direct
correspondence between $\Omega_{\chi}h^{2}$ and the rescaled annihilation
cross-section $f_{\chi}^{2}\langle\sigma v\rangle_{0}$. This is because the
operators that induce $p$-wave annihilations (in particular
$\mathcal{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 velocity-
dependent 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 Acharyya:2020sbj
(for simplicity based on the assumption of $b\bar{b}$ 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 Acharyya:2020sbj
.

Let us finally consider the allowed parameter space in terms of the Wilson
coefficients. The coefficient $\mathcal{C}^{(6)}_{1}$ 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
$\mathcal{C}^{(6)}_{4}$, on the other hand, gives rise to spin-dependent
interactions, for which constraints are significantly weaker. We show the
allowed parameter regions for this coefficient in the left panel of Fig. 4.
The observed mirror symmetry results from the fact that all experimental
predictions (and hence the likelihoods) are invariant under a global sign
change of all Wilson coefficients. For ${\rm\Lambda}<200\,\mathrm{GeV}$, all
constraints are furthermore invariant under the rescaling
$\mathcal{C}\to\alpha^{2}\mathcal{C}$, ${\rm\Lambda}\to\alpha{\rm\Lambda}$,
which explains why the allowed parameter region grows with increasing
${\rm\Lambda}$. For ${\rm\Lambda}>200\,\mathrm{GeV}$, LHC constraints become
relevant and strongly constrain the magnitude of the coefficient. Very similar
results are obtained for the coefficient $\mathcal{C}^{(6)}_{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
$\mathcal{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 coefficient is largely
unconstrained for ${\rm\Lambda}<200\,\mathrm{GeV}$. We also identify this
coefficient as giving the main contribution for fitting the _Fermi_ -LAT
excess. For larger values of ${\rm\Lambda}$, on the other hand, the
constraints are very similar to the ones for $\mathcal{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_{\chi}$–${\rm\Lambda}$ 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_{\chi}\lesssim 100\,\mathrm{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 ($\mathcal{Q}_{1,q}^{(6)}$ and
$\mathcal{Q}_{3,q}^{(6)}$), and incompatible with direct detection and LHC
constraints for $\mathcal{Q}_{2,q}^{(6)}$ and $\mathcal{Q}_{4,q}^{(6)}$.

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_{\chi}=500\,\mathrm{GeV}$. We find the likelihood of the best-fit
point to be slightly worse than for the background-only hypothesis:
$2\Delta\ln\mathcal{L}\equiv 2(\ln\mathcal{L}^{\mathrm{best-
fit}}-\ln\mathcal{L}^{\text{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 best-fit 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 Akerib:2018lyp  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\,\mathrm{keV},30\,\mathrm{keV}]$ and assumes an exposure of $5.6\times
10^{6}\,\mathrm{kg\,days}$ and 50% acceptance for nuclear recoils (see Ref. HP
for details on our implementation of LZ).

Figure 5: Same as the right panel of Fig. 2 but requiring the DM relic density
to be saturated (rather than imposing an upper bound only).

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
$\mathcal{Q}_{2}^{(6)}$ leading to spin-independent (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 $\mathcal{Q}_{3}^{(6)}$ into the operator $\mathcal{Q}_{1}^{(6)}$
as given in Eq. (21).171717We note that this mixing effect could in principle
be cancelled by contributions from additional effective operators not included
in our analysis, such that even smaller event rates may be achievable. 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 Aalbers:2016jon
or DarkSide-20k Aalseth:2017fik .

Another interesting approach would be to not perform a relic density
calculation at all and simply assume that $f_{\chi}=1$ is achieved through
some modification of early universe cosmology. In this case it would also be
possible to consider ${\rm\Lambda}<2m_{\chi}$ 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_{\chi}$ and
${\rm\Lambda}$ to be viable in this approach, and we do not explore this
direction further in the present work.

Figure 6: Profile likelihood in terms of $m_{\chi}$ and the predicted number
of signal events in the LZ experiment when $\chi$ accounts for all of the
observed DM abundance (as in Fig. 5).

Figure 7: Profile likelihood in the $m_{\chi}$–${\rm\Lambda}$ plane (left) and
in terms of $m_{\chi}$ and the predicted relic density (right) when
considering all dimension-6 and dimension-7 operators, and capping the LHC
likelihood at the value of the background-only hypothesis. The white contours
show the $1\sigma$ and $2\sigma$ confidence regions and the white star marks
the best-fit point. For comparison, we also show the $1\sigma$ and $2\sigma$
confidence region contours (dashed grey lines) and best-fit point (grey star)
for the case of dimension-6 operators only in the right panel (see also Fig.
3).

#### 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_{\chi}$-${\rm\Lambda}$ plane
(left) and in the $m_{\chi}$–$\Omega_{\chi}h^{2}$ plane (right) when using the
capped LHC likelihood. As before, we find that the parameter region at small
$m_{\chi}$ and ${\rm\Lambda}$ can fit the slight _Fermi_ -LAT excess with
best-fit values: $m_{\chi}=5.5$ GeV and $f_{\chi}^{2}\langle\sigma
v\rangle_{0}=1.9\times 10^{-27}$ cm3 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_{\chi}$ and ${\rm\Lambda}$. This is of course
expected for the parameter region with large $m_{\chi}$ and small
${\rm\Lambda}$ (bottom-right), which is excluded by the EFT validity
constraint but surprising for the region with small $m_{\chi}$ and large
${\rm\Lambda}$ (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 $\mathcal{Q}_{1\text{--}4}^{(7)}$ are strongly constrained by the
LHC for ${\rm\Lambda}>200\,\mathrm{GeV}$ and can therefore not contribute
significantly to the annihilation cross-section. The quark operators
$\mathcal{Q}_{5\text{--}10,q}^{(7)}$, on the other hand, are unconstrained by
the LHC, but for $m_{\chi}<m_{t}$, the resulting annihilation cross-section is
suppressed by a factor $m_{b}^{2}m_{\chi}^{2}/{\rm\Lambda}^{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.

Figure 8: Profile likelihood in the $m_{\chi}$–${\rm\Lambda}$ parameter plane
when considering only dimension-6 operators and including the full LHC
likelihood. In the left (right) panel, we impose a hard (smooth) cut-off in
the predicted missing energy spectrum for $\not{E}_{T}>{\rm\Lambda}$ (see text
for details). Figure 9:  Top panel: Examples of missing energy spectra for the
CMS monojet search Sirunyan:2017jix , illustrating different choices for
imposing an EFT validity requirement on the signal prediction. For
$\not{E}_{T}>{\rm\Lambda}$, we scale the $\not{E}_{T}$ signal spectrum with
the factor $(\not{E}_{T}/{\rm\Lambda})^{-a}$ as described in Sec. 2.2. The
green distributions show the resulting signal predictions for four different
choices of $a$, from $a=0$ (lightest green), corresponding to no modification
of the spectrum, to $a\rightarrow\infty$ (darkest green), which removes any
signal contribution in $\not{E}_{T}$ bins above $\not{E}_{T}={\rm\Lambda}$.
The SM background prediction (purple) and the observed event counts (black
points) are taken from Ref. Sirunyan:2017jix . The last bin, starting at
$\not{E}_{T}=1250\,\mathrm{GeV}$, contains any overflow and is thus not
normalised to a given bin width. Bottom panel: A small bar chart per bin
showing the pulls, defined as
$(\textrm{data}-\textrm{prediction})/\textrm{uncertainty}$, resulting from
adding the indicated signal prediction on top of the SM background prediction.
The uncertainty includes the background uncertainty, signal uncertainty and
statistical data uncertainty, added in quadrature. The purple bars show the
pulls when only including the SM background prediction. The $\chi^{2}$ values
in the top panel legend correspond to the sum of the squared pulls in each
case. These values are intended for illustration only, i.e. they do not
correspond directly to $-2\ln\mathcal{L}_{\text{CMS}}$ where
$\mathcal{L}_{\text{CMS}}$ is defined in Eq. (39).

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_{\chi}$ (and small ${\rm\Lambda}$), thanks to the
contribution of $\mathcal{Q}_{3,q}^{(7)}$ and $\mathcal{Q}_{7,q}^{(7)}$, which
both give suppressed signals in direct and indirect detection experiments and
are therefore largely unconstrained. Moreover, for $m_{\chi}>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 $\mathcal{Q}_{5\text{--}10,q}^{(7)}$. 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\,\mathrm{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 $\Omega_{\chi}h^{2}=0.12$ everywhere.

###  4.2 Full LHC likelihood

#### 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_{\chi}$ and ${\rm\Lambda}$ for the case where we introduce a hard cut-
off in the missing energy spectrum for $\not{E}_{T}>{\rm\Lambda}$ (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 ${\rm\Lambda}$ values. This
preference arises due to data excesses in a few high-$\not{E}_{T}$ bins in the
ATLAS and CMS monojet searches.

Figure 10: Same as Fig. 8 but requiring the DM relic abundance to match the
total observed DM abundance.

The difference in the above two results can be understood as follows. For
$\not{E}_{T}<{\rm\Lambda}$, the missing energy spectrum arising from DM is
harder than the background, while for $\not{E}_{T}>{\rm\Lambda}$, we either
set it to zero or assume that it drops rapidly. Thus, the ratio of signal-to-
background is largest for $\not{E}_{T}\approx{\rm\Lambda}$, 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\sigma$–$2\sigma$ data excesses in bins around $\not{E}\approx
700\,\mathrm{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 $\not{E}\approx
700\,\mathrm{GeV}$, the model will be penalized for causing larger pulls in
the highest-$\not{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\sigma$ level) separate parameter regions preferred by the CMS
analysis (${\rm\Lambda}\approx 700\,\mathrm{GeV}$) and the ATLAS analysis
(${\rm\Lambda}\gtrsim 1\,\mathrm{TeV}$), with the overall best-fit point
corresponding to the latter and being preferred relative to the background-
only hypothesis by $2\Delta\ln\mathcal{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 $(\not{E}_{T}/{\rm\Lambda})^{-1}$.
In this case, the best-fit point has $2\Delta\ln\mathcal{L}=2.6$.181818We note
that in both cases, the likelihood is very flat around the maximum and hence
the precise location of the best-fit point is somewhat arbitrary. 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 best-fit point is unaffected, while for
a smooth cut-off, it is pushed to slightly larger values of $m_{\chi}$ and
${\rm\Lambda}$. The respective preferences are reduced slightly to
$2\Delta\ln\mathcal{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 $\mathcal{C}_{2}^{(6)}$ or $\mathcal{C}_{4}^{(6)}$, the
relic density can only be reproduced with a contribution from
$\mathcal{C}_{3}^{(6)}$. This is because $\mathcal{Q}_{2,4}^{(6)}$ lead to
suppressed annihilation rates in the early universe, compared to
$\mathcal{Q}_{3}^{(6)}$, while $\mathcal{Q}_{1}^{(6)}$ 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 $\mathcal{C}_{3}^{(6)}$ 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.

LHC likelihood | Relic density | $\mathbf{2\Delta\ln\mathcal{L}}$ | Best-fit $m_{\chi}$ | Best-fit ${\rm\Lambda}$ | Best-fit constrained coupling
---|---|---|---|---|---
| constraint |  | (GeV) | (GeV) | combination(s) $\mathbf{(TeV^{-2})}$
Capped | Upper bound | $0.3$ | $5.0$ | $<200$ | $|\mathcal{C}_{3}^{(6)}|/{\rm\Lambda}^{2}=67$
Capped | Saturated | $-0.5$ | $500$ | $>1000$ |  | $|\mathcal{C}_{2}^{(6)}|/{\rm\Lambda}^{2}=0.22$
---
$|\mathcal{C}_{3}^{(6)}|/{\rm\Lambda}^{2}=0.041$
Full (hard cut-off) | Upper bound | $2.2$ | $500$ | $>1250$ | $|\mathcal{C}_{3}^{(6)}|/{\rm\Lambda}^{2}=0.14$
Full (smooth cut-off) | Upper bound | $2.6$ | $320$ | $640$ | $|\mathcal{C}_{3}^{(6)}|/{\rm\Lambda}^{2}=0.18$
Full (hard cut-off) | Saturated | $1.9$ | $500$ | $>1250$ |  | $|\mathcal{C}_{3}^{(6)}|/{\rm\Lambda}^{2}=0.047$
---
$\sqrt{(\mathcal{C}_{2}^{(6)})^{2}+(\mathcal{C}_{4}^{(6)})^{2}}/{\rm\Lambda}^{2}=0.15$
Full (smooth cut-off) | Saturated | 2.0 | 420 | 840 |  | $|\mathcal{C}_{3}^{(6)}|/{\rm\Lambda}^{2}=0.052$
---
$\sqrt{(\mathcal{C}_{2}^{(6)})^{2}+(\mathcal{C}_{4}^{(6)})^{2}}/{\rm\Lambda}^{2}=0.23$
Table 5: Best-fit points from our various scans involving dimension-6
operators with restricted parameter ranges ($5\,\mathrm{GeV}\leq m_{\chi}\leq
500\,\mathrm{GeV}$ and $20\,\mathrm{GeV}\leq{\rm\Lambda}\leq
2\,\mathrm{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.

#### 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 ${\rm\Lambda}\sim
80\,\mathrm{GeV}$ together with a rather soft cut-off $a\approx 1.7$ of the
$\not{E}_{T}$ spectrum above ${\rm\Lambda}$. The resulting best-fit point has
$2\Delta\ln\mathcal{L}=2.9$, which is the highest likelihood found in any of
our scans.

Figure 11: Profile likelihood in the $m_{\chi}$–${\rm\Lambda}$ parameter plane
when considering dimension-6 and dimension-7 operators and including the full
LHC likelihood. In the left (right) panel, we impose a hard (smooth) cut-off
in the predicted missing energy spectrum for $\not{E}_{T}>{\rm\Lambda}$ (see
text for more details).

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 $\mathcal{Q}_{4}^{(7)}$ is essentially
unconstrained by direct detection and can induce sizeable LHC signals if
$\mathcal{C}_{4}^{(7)}$ 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
$\mathcal{Q}_{5\text{--}8,q}^{(7)}$. As a result, the allowed regions in
$m_{\chi}$–${\rm\Lambda}$ parameter space look very similar to the ones shown