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 Mathematics and Physics, The University of
Queensland, St. Lucia, Brisbane, QLD 4072, Australia22institutetext:  Quantum
Brilliance Pty Ltd, The Australian National University, Daley Road, Acton ACT
2601, Australia33institutetext:  Institute for Theoretical Particle Physics
and Cosmology (TTK), RWTH Aachen University, D-52056 Aachen,
Germany44institutetext:  Institute for Theoretical Particle Physics (TTP),
Karlsruhe Institute of Technology (KIT), 76128 Karlsruhe,
Germany55institutetext:  ARC Centre of Excellence for Dark Matter Particle
Physics & CSSM, Department of Physics, University of Adelaide, Adelaide, SA
500566institutetext:  Department of Physics, University of Oslo, N-0316 Oslo,
Norway

<EMAIL_ADDRESS>
# Global fits of simplified models for dark matter with GAMBIT

I. Scalar and fermionic models with $s$-channel vector mediators

Christopher Changuq,a  Pat Scottqb  Tomás E. Gonzaloaachen,kitTTP  Felix
Kahlhoeferaachen,kitTTP  Anders Kvellestadoslo  Martin Whiteadelaide

(Received: date / Accepted: date)

###### Abstract

Simplified models provide a useful way to study the impacts of a small number
of new particles on experimental observables and the interplay of those
observables, without the need to construct an underlying theory. In this
study, we perform global fits of simplified dark matter models with GAMBIT
using an up-to-date set of likelihoods for indirect detection, direct
detection and collider searches. We investigate models in which a scalar or
fermionic dark matter candidate couples to quarks via an $s$-channel vector
mediator. Large parts of parameter space survive for each model. In the case
of Dirac or Majorana fermion dark matter, excesses in LHC monojet searches and
relic density limits tend to prefer the resonance region, where the dark
matter has approximately half the mass of the mediator. A combination of
vector and axial-vector couplings to the Dirac candidate also leads to
competing constraints from direct detection and unitarity violation.

††journal: Eur. Phys. J. C

###### Contents

  1. 1 Introduction
  2. 2 Models
    1. 2.1 Scalar DM
    2. 2.2 Dirac Fermion DM
    3. 2.3 Majorana Fermion DM
  3. 3 Constraints
    1. 3.1 Relic Abundance
    2. 3.2 Direct Detection
    3. 3.3 Indirect Detection
    4. 3.4 Collider searches for WIMPs using monojet events
    5. 3.5 Collider searches for the mediator using dijet events
    6. 3.6 Nuisance parameter likelihoods
  4. 4 Results
    1. 4.1 Scalar DM
    2. 4.2 Dirac Fermion DM
    3. 4.3 Majorana Fermion DM
    4. 4.4 Future prospects
  5. 5 Conclusions
  6. 6 Acknowledgements

##  1 Introduction

The Standard Model (SM) remains enormously successful as a theory of particle
physics, but is widely thought to be incomplete and expected to be superseded
by a more complete theory. One of the many motivations for searching for
beyond-Standard Model (BSM) physics is to explain the dark matter (DM) evident
in a number of astrophysical and cosmological observations Zwicky33 ; bullet ;
wmap3year . The WIMP hypothesis, in which dark matter is assumed to consist of
a new species of Weakly-Interacting Massive Particle, is amongst the leading
DM explanations due to its ability to explain the observed cosmological relic
abundance of dark matter Lee:1977ua  at the same time as potentially being
very tightly constrained in the near future by current experimental
technologies Arcadi:2017kky .

Whilst there are plenty of UV-complete theories that include viable WIMP
candidates, it is advantageous to take a model-independent approach and
construct a low-energy effective theory that includes the relevant
phenomenology for our current experimental probes, whilst remaining agnostic
about the high energy physics that we cannot currently observe. The simplest
way to construct such a theory is to write down an _effective field theory_
(EFT), in which the SM Lagrangian density is extended with a number of
effective operators that encode possible DM-SM interactions. An effective
field theory is valid up to some scale ${\rm\Lambda}$, at which point one
would start to resolve the physics that generates the operators in the low
energy description. This typically makes EFTs useful for studying the impact
of low energy experimental probes, such as indirect Goodman:2010qn ;
Beltran:2008xg ; Cheung:2011nt ; Harnik:2008uu ; DeSimone:2013gj ;
Karwin:2016tsw  and direct Fan:2010gt ; Agrawal:2010fh ; Crivellin:2014gpa ;
Hoferichter:2016nvd ; Kahlhoefer:2016eds  dark matter detection experiments,
but difficult for higher-energy probes such as the ATLAS and CMS experiments
at the Large Hadron Collider (LHC) Buchmueller:2013dya ; Abdallah:2015ter ;
Abdallah:2014hon ; Malik:2014ggr ; Busoni:2013lha ; Busoni:2014sya ;
Busoni:2014haa . A detailed study of EFTs for Dirac DM has recently been
performed by the GAMBIT collaboration GAMBIT:2021rlp .

An alternative, and only slightly more complicated, approach is to explicitly
add a particle that mediates interactions between the DM candidate and SM
species, giving what is popularly known as a _simplified model_. Most such
models are themselves also EFTs. Previous reviews of such models can be found
in Refs. Abdallah:2015ter ; Arina:2018zcq ; DeSimone:2016fbz ; Albert:2017onk
; Boveia:2016mrp ; Kahlhoefer:2017dnp ; Arcadi:2017kky ; Morgante:2018tiq . In
the limit of large mediator masses, the EFT approach may be recovered, but
simplified models remain the preferred way to investigate the simultaneous
impact of high energy collider and low energy DM probes DEramo:2016gos ;
Carpenter:2016thc ; Abercrombie:2015wmb  on the landscape of possible DM-SM
interactions.

There are a plethora of studies of simplified dark matter models in the
literature, including most notably a global scan of an $s$-channel vector-
mediated Dirac fermion dark matter model performed by the MasterCode
Collaboration Bagnaschi:2019djj . In this paper, we perform global fits of
$s$-channel vector-mediated scalar, Dirac and Majorana fermion dark matter
models, using GAMBIT v2.3 gambit . Setting up global fits of multiple
simplified models in a single study has recently become more accessible with
the tool for automatically generating GAMBIT code, GUM GUM . We complement the
work of Ref. Bagnaschi:2019djj  by also considering scalar and Majorana
fermion DM, and Dirac fermion DM with both vector and axial vector DM
couplings. The latter requires thorough treatment of direct detection
signatures to allow concurrent contribution from spin-dependent and spin-
indepent interactions.

Performing global fits of these models is of interest, due to the high degree
of interplay between different experiments. In the case of Dirac DM, it has
been shown in the EFT limit that requiring the relic abundance to be saturated
can be incompatible with constraints from laboratory experiments for low DM
masses GAMBIT:2021rlp . We find that this incompatibility remains when the EFT
is promoted to a simplified model. This lower bound on the DM mass is also
present in the complex scalar DM model, where collider searches have very
little influence on this bound. Suppressed direct detection signals in the
Majorana DM model give no such lower bound on the DM mass, and result in much
weaker exclusion as compared to the Dirac DM model. We also find that the
combination of both vector and axial-vector couplings in the Dirac DM model
leads to an interplay of constraints from direct detection and unitarity
violation. Exclusion from strong spin-independent direct detection signals,
given a vector coupling, overlaps with unitarity violation from the axial-
vector coupling to give an exclusion that would not be present when fixing
either to zero.

This paper is structured as follows. In Section 2, we describe the different
simplified dark matter models that we study. Section 3 contains a detailed
description of the experimental likelihoods that we include in our analysis.
Section 4 contains the details of our global fit technique and our results,
and Section 5 presents our conclusions. The samples from our scans and the
corresponding GAMBIT input files and plotting scripts can be downloaded from
Zenodo Zenodo_DMsimp .

##  2 Models

We consider three models in this work, distinguished by the nature of the DM
candidate:

  * •

spin $0$ (complex scalar)

  * •

spin $\frac{1}{2}$ (Dirac)

  * •

spin $\frac{1}{2}$ (Majorana)

In all cases, we assume that DM is a singlet under the SM gauge group. All
three models have a spin $1$ (vector) mediator. A thorough treatment of vector
DM with a vector mediator requires the derivation of new unitarity limits, and
we thus defer this option to Paper II of this series DMSimpII . To enforce
absolute stability of the DM candidate, we assume it to be odd under a new
$\mathbb{Z}_{2}$ symmetry under which the SM fields and the mediator are even.
We take the mediator to be a real vector. This could typically arise in an
extension of the SM gauge group by an abelian symmetry group such as U(1). We
assume that none of the models gives rise to observable mixing between the
mediator and any SM particles; this is required in order to reinterpret most
of the experimental results that we consider. The models that we consider also
assume that the mass generation mechanism that applies to the new particles
has no observable effect on any of the experimental results that we consider.
This could be achieved by, for example, the existence of a dark Higgs
mechanism where the mass scale of the dark Higgs is well above both the
mediator and dark matter masses. A detailed study of the case where this
assumption breaks down for one model can be found in Duerr:2016tmh . The
models that we consider remain valid provided that the details of the
encompassing higher order theory remain sufficiently decoupled at energy
scales probed by current experiments.

For each model, we assume purely vector couplings to quarks (no axial-vector
couplings) to prevent severe constraints from electroweak precision tests, and
no lepton couplings111Leptophobic mediators could be accommodated in some
extensions of the SM gauge group, such as gauged baryon number Duerr:2014wra
or anomaly-free extensions that require multiple new singlet fermions
Ellis:2017tkh . to avoid restrictive dilepton searches Kahlhoefer:2015bea . We
also assume Minimal Flavour Violation to avoid constraints from flavour
physics, and for simplicity, we assume equal couplings between up and down-
type quarks. These together require universal quark couplings
($g_{\mathrm{q}}$) to the mediator.

###  2.1 Scalar DM

The Lagrangian density of our simplified model with a complex scalar DM
candidate is

$\displaystyle\begin{split}\mathcal{L}_{\mathrm{BSM}}=&\partial_{\mu}\phi^{\dagger}\partial^{\mu}\phi-
m_{\mathrm{DM}}^{2}\phi^{\dagger}\phi\\\
-&\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{1}{2}m_{\mathrm{M}}^{2}V_{\mu}V^{\mu}+g_{\mathrm{q}}V_{\mu}\bar{q}\gamma^{\mu}q\\\
+&ig^{\mathrm{V}}_{\mathrm{DM}}V_{\mu}\Big{(}\phi^{\dagger}(\partial^{\mu}\phi)-(\partial^{\mu}\phi^{\dagger})\phi\Big{)}\,,\end{split}$
(1)

where $V_{\mu}$ is the mediator field, $\phi$ is the scalar DM candidate and
$F_{\mu\nu}$ is the mediator field strength tensor. The model has four free
parameters: the DM mass $m_{\mathrm{DM}}$, the mediator mass $m_{\mathrm{M}}$,
the mediator-WIMP coupling $g_{\mathrm{DM}}^{\mathrm{V}}$ and the mediator-
quark coupling $g_{\mathrm{q}}$. The DM candidate must be a complex scalar to
avoid vanishing $g^{\mathrm{V}}_{\mathrm{DM}}$.

The decay width of the mediator to quark $q$, for all three models is

$\displaystyle\Gamma(V_{\mu}\rightarrow\bar{q}q)=\frac{|g_{\mathrm{q}}|^{2}m_{\mathrm{M}}}{4\pi}\sqrt{1-\frac{4m_{\mathrm{q}}^{2}}{m_{\mathrm{M}}^{2}}}\Bigg{(}1+\frac{2m_{\mathrm{q}}^{2}}{m_{\mathrm{M}}^{2}}\Bigg{)}\,,$
(2)

while the decay width to the scalar DM candidate is

$\displaystyle\Gamma(V_{\mu}\rightarrow\phi^{\dagger}\phi)=\frac{|g_{\mathrm{DM}}^{\mathrm{V}}|^{2}m_{\mathrm{M}}}{48\pi}\Bigg{(}1-\frac{4m_{\mathrm{DM}}^{2}}{m_{\mathrm{M}}^{2}}\Bigg{)}^{3/2}.$
(3)

###  2.2 Dirac Fermion DM

The simplified model with Dirac fermion DM has the Lagrangian density

$\displaystyle\begin{split}\mathcal{L}_{\mathrm{BSM}}=&i\bar{\chi}\gamma^{\mu}\partial_{\mu}\chi-
m_{\mathrm{DM}}\bar{\chi}\chi\\\
-&\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{1}{2}m_{\mathrm{M}}^{2}V_{\mu}V^{\mu}+g_{\mathrm{q}}V_{\mu}\bar{q}\gamma^{\mu}q\\\
+&V_{\mu}\bar{\chi}(g^{\mathrm{V}}_{\mathrm{DM}}+g^{\mathrm{A}}_{\mathrm{DM}}\gamma^{5})\gamma^{\mu}\chi\,,\end{split}$
(4)

where $\chi$ represents the fermion DM candidate and the mediator and quark
fields are given as before. The model has five free parameters, including the
WIMP and mediator masses and the mediator-quark coupling, defined as before.
Note, however, that this time there are two possible mediator-WIMP couplings,
representing vector ($g^{\mathrm{V}}_{\mathrm{DM}}$) and axial-vector
($g^{\mathrm{A}}_{\mathrm{DM}}$) couplings. We vary
$g^{\mathrm{V}}_{\mathrm{DM}}$ and $g^{\mathrm{A}}_{\mathrm{DM}}$
independently in our scans, allowing for possible interference effects between
interactions mediated by the two.

The decay width to a given quark is given by Eq. 2 while the decay width to
the Dirac fermion DM candidate is

$\displaystyle\begin{split}\Gamma(V_{\mu}&\rightarrow\bar{\chi}\chi)=\frac{m_{\mathrm{M}}}{12\pi}\sqrt{1-\frac{4m_{\mathrm{DM}}^{2}}{m_{\mathrm{M}}^{2}}}\\\
&\times\Bigg{(}|g^{\mathrm{V}}_{\mathrm{DM}}|^{2}\Big{(}1+\frac{2m_{\mathrm{DM}}^{2}}{m_{\mathrm{M}}^{2}}\Big{)}+|g^{\mathrm{A}}_{\mathrm{DM}}|^{2}\Big{(}1-\frac{4m_{\mathrm{DM}}^{2}}{m_{\mathrm{M}}^{2}}\Big{)}\Bigg{)}.\end{split}$
(5)

The axial-vector coupling in the fermion DM model implies that perturbative
unitarity may be violated. The violation of unitarity stems from the omission
of a dark Higgs boson in the simplified model and is analagous to the
unitarity violation in the Standard Model without a Higgs boson. We adopt the
bound of Ref. Kahlhoefer:2015bea  for this:

$\displaystyle
m_{\mathrm{DM}}\leq\sqrt{\frac{\pi}{2}}\frac{m_{\mathrm{M}}}{g_{\mathrm{DM}}^{\mathrm{A}}}.$
(6)

We apply this as an additional constraint, rejecting as theoretically invalid
any parameter points that do not satisfy Eq. 6.

###  2.3 Majorana Fermion DM

The model Lagrangian of the BSM additions in the case of Majorana fermion DM
is:

$\displaystyle\begin{split}\mathcal{L}_{\mathrm{BSM}}=&\frac{1}{2}i\bar{\psi}\gamma^{\mu}\partial_{\mu}\psi-\frac{1}{2}m_{\mathrm{DM}}\bar{\psi}\psi\\\
-&\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{1}{2}m_{\mathrm{M}}^{2}V_{\mu}V^{\mu}+g_{\mathrm{q}}V_{\mu}\bar{q}\gamma^{\mu}q\\\
+&\frac{1}{2}g^{\mathrm{A}}_{\mathrm{DM}}V_{\mu}\bar{\psi}\gamma^{5}\gamma^{\mu}\psi\,,\end{split}$
(7)

Unlike the Dirac Fermion case, Majorana DM prevents pure vector couplings to
the mediator, and so this model has four free parameters. Like the Dirac
fermion case, perturbative unitarity will be violated due to the presence of
axial-vector couplings. The same bound (Eq. 6) on parameters is applicable in
this case, to exclude regions where unitarity is violated.

The decay width to a given quark is given by Eq. 2 while the decay width to
the Majorana fermion DM candidate is

$\displaystyle\Gamma(V_{\mu}\rightarrow\bar{\psi}\psi)=\frac{|g^{\mathrm{A}}_{\mathrm{DM}}|^{2}m_{\mathrm{M}}}{12\pi}\Bigg{(}1-\frac{4m_{\mathrm{DM}}^{2}}{m_{\mathrm{M}}^{2}}\Bigg{)}^{3/2}.$
(8)

##  3 Constraints

DM-SM interactions are constrained by a broad range of astrophysical,
cosmological and particle physics measurements.

We have implemented likelihoods for DM direct and indirect detection
experiments, the measurement of the DM relic abundance and collider searches
with the ATLAS and CMS experiments. For this we employ GUM GUM , the GAMBIT
Universal Model Machine. GUM generates the necessary model-specific module
functions for GAMBIT 2.3, along with interfaces to the relevant backend codes
that contain calculations for each dark matter observable. For the fermionic
DM models, we supplement these likelihoods with the unitarity bound of Eq. 6.

In Table 1 we provide a summary of the likelihoods for different measurements
sensitive to BSM physics included in our scans. For searches where BSM effects
are sought above some existing SM-induced event rate, we also give the value
$\ln\mathcal{L}^{\text{bg}}$ that each likelihood takes in the pure SM case,
where BSM physics is absent and only background events are observed. In other
cases, i.e. likelihoods corresponding to values of parameters of the SM
itself, or of the relic density of dark matter, we simply give the best-case
likelihood that results from predictions or parameters exactly matching their
centrally measured values.

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$ |
LZ 2022 LZ:2022ufs  | 0 |
LHC Dijets CMS:2019gwf ; ATLAS:2019fgd ; ATLAS:2018qto ; CDF:2008ieg ; ATLAS:2018hbc ; ATLAS:2018hzj ; CMS:2019emo ; ATLAS:2019itm ; CMS:2019xai  | $0$ |
ATLAS monojet Aad:2021egl  | $0$ |
CMS monojet CMS:2021snz  | $0$ |
_Fermi_ -LAT LATdwarfP8  | $-33.245$ |
_Planck_ 2018: $\Omega h^{2}$ Aghanim:2018eyx  |  | $5.989$
Nuisances (see Table 3) |  | $6.728$
Table 1: A summary of all likelihoods included in our scans. For likelihoods
that search for DM events above an SM background, we also give the SM-only
(i.e. background-only) log-likelihood $\ln\mathcal{L}^{\text{bg}}$ for
comparison. For the remaining likelihoods, where there is no ‘SM background‘
to sensibly speak of, we give the maximum achievable value of the log-
likelihood $\ln\mathcal{L}^{\text{max}}$, which corresponds to the case of a
‘perfect fit’ where the predicted value of the relic density or the chosen
value of a nuisance parameter is exactly equal to its measured value.

The details for these likelihoods, and their implementation in GAMBIT, are
given in the following subsections.

###  3.1 Relic Abundance

The number density $n_{\scaleto{\mathrm{DM}}{3pt},\textrm{eq}}$ of DM
particles in thermal equilibrium in the early Universe changes with time
according to the Boltzmann equation Gondolo:1990dk

$\displaystyle\frac{dn_{\scaleto{\mathrm{DM}}{3pt}}}{dt}+3Hn_{\scaleto{\mathrm{DM}}{3pt}}=-\langle\sigma
v_{\textrm{rel}}\rangle\left(\left(n_{\scaleto{\mathrm{DM}}{3pt}}n_{\bar{\scaleto{\mathrm{DM}}{3pt}}}\right)-\left(n_{\scaleto{\mathrm{DM}}{3pt}}n_{\bar{\scaleto{\mathrm{DM}}{3pt}}}\right)_{\textrm{eq}}\right)\,,$
(9)

where $H(t)$ is the Hubble rate and
$n_{\scaleto{\mathrm{DM}}{3pt}}n_{\bar{\scaleto{\mathrm{DM}}{3pt}}}\equiv
n_{\scaleto{\mathrm{DM}}{3pt}}^{2}$ for Majorana DM. $\langle\sigma
v_{\textrm{{rel}}}\rangle$ is the thermally averaged cross-section times the
relative velocity and is given by

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

where $K_{1,2}$ are the modified Bessel functions and $v_{\rm lab}$ is the
velocity of one of the annihilating DM particles or antiparticles in the rest
frame of the other (see Ref. Binder:2017rgn  for further discussion). Dirac
fermion WIMPs give a total contribution to the observed DM density of
$n_{\scaleto{\mathrm{DM}}{3pt}}+n_{\bar{\scaleto{\mathrm{DM}}{3pt}}}=2n_{\scaleto{\mathrm{DM}}{3pt}}$
(where a possible initial asymmetry has been neglected Kaplan:2009ag ).

We use GUM to generate tree-level annihilation cross-sections for each model
using CalcHEP v3.6.27 Pukhov:2004ca ; Belyaev:2012qa . We obtain the WIMP
relic density by using the DarkBit interface to DarkSUSY v6.2.5 darksusy ;
darksusy4  to numerically solve Eq. (9) at each parameter point, assuming no
departures from the standard cosmological history. For alternative scenarios,
see Ref. Arbey:2021gdg . We compare the resulting value to the _Planck_ 2018
measurement of $\Omega_{\textrm{DM,obs}}\,h^{2}=0.120\pm 0.001$
Aghanim:2018eyx . We include a $1\%$ theoretical error on our computed relic
density values, added in quadrature with the quoted _Planck_ uncertainty. More
details on this prescription can be found in Refs. gambit ; DarkBit .

Given the Planck measurement of the DM relic abundance, it is interesting to
consider both the case where our hypothesised WIMPs constitute all of DM, and
the case where they form only a subcomponent. In the former case, we use the
Planck measurement to define a Gaussian likelihood based on the predicted WIMP
abundance. In the latter case, we modify the likelihood such that it is flat
if the predicted value is smaller than the observed one; details can be found
in Ref. gambit .

We rescale all predicted direct detection signals by the fraction

$\displaystyle
f_{\textrm{DM}}\equiv\Omega_{\textrm{DM}}/\Omega_{\mathrm{DM,obs}}.$ (11)

where
$\Omega_{\textrm{DM}}\equiv\Omega_{\phi}+\Omega_{\phi^{\dagger}}=2\Omega_{\phi}$
for scalar DM,
$\Omega_{\textrm{DM}}\equiv\Omega_{\chi}+\Omega_{\bar{\chi}}=2\Omega_{\chi}$
for Dirac DM and $\Omega_{\textrm{DM}}\equiv\Omega_{\psi}$ for Majorana DM.
Similarly, we rescale all indirect detection signals from DM-DM annihilation
by $f_{\mathrm{DM}}^{2}$. Note that this encodes the two assumptions that
$f_{\mathrm{DM}}$ is the same in all astrophysical systems, and the
conservative assumption that any extra DM components are not visible to any of
the astrophysical experiments that we consider.

###  3.2 Direct Detection

Direct detection experiments search for the scattering of DM particles off
nuclei in some highly pure detector material, by measuring the nuclear recoil
energy spectrum. The differential event rate with respect to the recoil energy
$E_{R}$ is

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

where $\rho_{0}$ is the local density of DM at our point in the Galactic halo,
$m_{\mathrm{T}}$ is mass of each target nucleus and $f(v)$ describes the local
velocity distribution of the DM. Nuclear recoils of energy $E_{R}$ occur above
some minimum DM velocity, which is given by

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

with
$\mu=m_{\mathrm{T}}\,m_{\scaleto{\mathrm{DM}}{3pt}}/(m_{\mathrm{T}}+m_{\scaleto{\mathrm{DM}}{3pt}})$
the reduced mass of the DM-nucleus system.

Both the local density and the local velocity distribution of DM are subject
to sizeable uncertainties. We therefore include the relevant quantities as
nuisance parameters in our scans (see Section 3.6). The remaining part of Eq.
12 is the differential scattering cross-section
$\mathrm{d}\sigma/\mathrm{d}E_{\mathrm{R}}$. In calculating this, one can make
use of the knowledge that the scattering of WIMPs from the Galactic halo on
nuclei is non-relativistic, typically involving momentum transfers below 200
MeV. If the mediator mass is much heavier than this scale (as assumed in our
study), the mediator can be integrated out of the theory, and its contribution
to scattering processes can instead be modelled using a non-relativisitic EFT
written as

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

where the higher-dimensional operators $\mathcal{O}^{\mathrm{N}}_{\mathrm{i}}$
depend only on the spins of the WIMP and the nucleon
($\vec{S}_{\scaleto{\mathrm{DM}}{3pt}}$ and $\vec{S}_{\mathrm{N}}$), the
momentum transfer $\vec{q}$ and the DM-nucleon relative velocity $\vec{v}$
Fitzpatrick:2012ix ; Anand:2013yka ; Fan:2010gt . A systematic classification
of the non-relativistic operators that can result from the reduction of
simplified models can be found in Refs Dent:2015zpa ; Fitzpatrick:2012ix ,
giving 16 operators in total.

One can then proceed by translating the parameters of the simplified model to
the coefficients of the relevant operators
$c_{\mathrm{i}}^{\mathrm{N}}(q^{2})$ in the non-relativistic EFT. The relevant
operators for our simplified models are provided in Table 2. The associated
prefactors are passed to DDCalc v2.2.0 DarkBit ; HP  (through the DirectDM
code, following Appendix A in GAMBIT:2021rlp ). We use DDCalc to compute the
differential cross-section for each operator
$\mathcal{O}^{\mathrm{N}}_{\mathrm{i}}$ and target element of interest, and to
perform the velocity integral in Eq. 12 in order to obtain the differential
event rate. DDCalc then implements detector effects and computes the final
predicted number of events $N_{\mathrm{p}}$ at each experiment, by convolving
the differential rate with the product $\phi(E_{\mathrm{R}})$ of the energy
resolution and detector acceptance,

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

where $M$ is the detector mass and $T_{\text{exp}}$ is the exposure time.

In this work, we use direct detection data from the most recent XENON1T
analysis Aprile:2018dbl , LUX 2016 LUXrun2 , PandaX 2016 and 2017 Tan:2016zwf
; Cui:2017nnn , CDMSlite Agnese:2015nto , CRESST-II and CRESST-III
Angloher:2015ewa ; Abdelhameed:2019hmk , PICO-60 2017 and 2019 Amole:2017dex ;
Amole:2019fdf , DarkSide-50 Agnes:2018fwg  and LZ LZ:2022ufs .

Interaction | Effective Operator | Relevant models
---|---|---
scalar | $1_{\mathrm{DM}}1_{\mathrm{N}}$ | Scalar
vector | $1_{\mathrm{DM}}1_{\mathrm{N}}$ | Dirac
axial-vector |  $i\hat{\textbf{S}}\cdot\left(\hat{\textbf{S}}_{\mathrm{N}}\times\frac{\hat{q}}{m_{\mathrm{N}}}\right)$, $\hat{\textbf{S}}\cdot\hat{\textbf{v}}^{\perp}1_{\mathrm{N}}$ | Dirac, Majorana
Table 2: Effective Operator matching to each model. The axial-vector couplings
are momentum or velocity suppressed. Couplings to each operator follow
Baum:2017kfa .

###  3.3 Indirect Detection

The processes that led to DM being in thermal equilibrium in the early
Universe also allow for annihilation in the present day. Annihilation products
can thus potentially be seen originating from regions of high dark matter
density. Of the possible products, gamma rays are particularly useful to
search for, given that they should point back to the source without
deflection. In recent years, both satellite and ground-based gamma ray
observations have been used to constrain the possible interactions of dark
matter.

An especially strong set of constraints on annihilating dark matter comes from
observations of dwarf spheroidal galaxies Bringmann:2012ez . For a binned
histogram of the DM-induced $\gamma$-ray flux, the flux from a target $k$ in
bin $i$ can be written in the form $\Phi_{i}\cdot J_{k}$, where $\Phi_{i}$
includes the relevant particle physics, and $J_{k}$ includes the relevant
astrophysics (see Ref DarkBit  for more details).

The velocity dependence of the annihilation cross-section

$\displaystyle\langle\sigma v\rangle=a+bv^{2}+\mathcal{O}(v^{4}),$ (16)

plays a strong role in the sensitivity of indirect detection searches to
different models. If $a=0$, the leading-order term in the annihilation cross-
section is proportional to $v^{2}$ ($p$-wave annihilation), and the low
average DM velocities in the dwarf spheroidal galaxies will result in a
suppression of the predicted gamma-ray signals, compared to velocity-
independent $s$-wave annihilation (where $a\neq 0$). The two primary channels
of annihilation are through the s-channel to a pair of quarks, and the
t-channel to a pair of mediator particles.

In the case of scalar DM, annihilation to quarks is $p$-wave and annihilation
to mediators is $s$-wave. When the latter channel is open
($m_{\mathrm{DM}}>m_{\mathrm{M}}$), it will dominate gamma-ray signals. The
Dirac DM model has dominant $s$-wave annihilations to both quarks and
mediators. Without vector couplings to DM or axial-vector quark couplings, the
Majorana DM model’s annihilation into quarks has no $s$-wave contribution.
Like the scalar DM model, annihilation is $s$-wave to mediators when
$m_{\mathrm{DM}}>m_{\mathrm{M}}$ Arcadi:2017kky . Only $s$-wave contributions
are large enough to impact searches towards dwarf spheroidals in the models
that we consider, so we do not include the $p$-wave contributions in our
gamma-ray flux predictions. The Dirac fermion is therefore the only one of our
models that leads to significant indirect detection signals when
$m_{\mathrm{DM}}<m_{\mathrm{M}}$.

The $s$-wave contribution gives a particle physics factor of

$\displaystyle\Phi_{i}$
$\displaystyle=\frac{f_{\mathrm{DM}}^{2}}{N_{\mathrm{DM}}}\sum_{j}\frac{(\sigma
v)_{0,j}}{4\pi m_{\scaleto{\mathrm{DM}}{3pt}}^{2}}\int_{\Delta
E_{i}}dE\,\frac{dN_{\gamma,j}}{dE}\,,$ (17)

where $(\sigma v)_{0,j}$ denotes the zero-velocity limit of the cross-section
for $\mathrm{DM},\bar{\mathrm{DM}}\to j$, $N_{\gamma,j}$ is the number of
photons that results from the final state channel $j$ (per annihilation), and
$f_{\mathrm{DM}}$ is the DM fraction defined in Eq. 11. The prefactor of
$1/N_{\mathrm{DM}}$ reflects the nature of the DM candidate:
$N_{\mathrm{DM}}=2$ for self-conjugate particles and $N_{\mathrm{DM}}=4$
otherwise (assuming
$n_{\scaleto{\mathrm{DM}}{3pt}}=n_{\bar{\scaleto{\mathrm{DM}}{3pt}}}$).

We use CalcHEP to compute the annihilation cross-sections for each model,
using the GUM interface. Each annihilation channel can produce either primary
or secondary photons; the yields ${dN_{\gamma,j}}/{dE}$ are provided by
DarkBit based on tabulated Pythia runs provided by DarkSUSY.

The astrophysics factor $J_{k}$ for each dwarf spheroidal galaxy $k$ is given
by the line-of-sight integral over the DM distribution, assuming an NFW DM
halo 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}\,,$ (18)

where $D_{k}$ is the distance to the galaxy. We use the Pass-8 combined
analysis of 15 dwarf spheroidal galaxies performed by the _Fermi_ -LAT
Collaboration using 6 years of LAT data LATdwarfP8 , computing the likelihood
with the DarkBit interface to gamLike v1.0.1. The log-likelihood
$\ln\mathcal{L}_{\rm{exp}}$ is constructed from the product $\Phi_{i}\cdot
J_{k}$ 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)\,.$ (19)

An additional likelihood contribution comes from treating the $J_{k}$ factors
of each dwarf spheroidal galaxy as nuisance parameters, giving
$\ln\mathcal{L}_{J}=\sum_{k}\ln\mathcal{L}(J_{k})$ DarkBit ; LATdwarfP8 . The
full likelihood, profiled over the J factors, is then given by

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

An alternative to dwarf spheroidal measurements is to look for evidence of
dark matter annihilation in the centre of our Galaxy. Although _Fermi_ -LAT
Galactic Centre limits are not nearly as robust as those from dwarf spheroidal
galaxies, the forthcoming Cherenkov Telescope Array (CTA) is expected to probe
thermally-produced WIMPs up to particle masses of several TeV Acharyya:2020sbj
. We briefly consider the future impact of CTA observations on the viable
parameter space of our simplified models in Section 4.4.

###  3.4 Collider searches for WIMPs using monojet events

The simplified models defined in Section 2 allow for pair production of WIMPs
in proton-proton collisions at the LHC. This process becomes visible if one of
the incoming partons radiates a jet through initial state radiation, giving a
potential signal of events with a single jet plus missing transverse energy
($\not{E}_{\mathrm{T}}$). In this study, we include CMS and ATLAS searches for
monojet events based respectively on $137\,\rm{fb}^{-1}$ CMS:2021snz  and
$139\,\rm{fb}^{-1}$ Aad:2021egl  of Run II data. We neglect other signatures
such as monophoton events, which are known to give weaker constraints on our
simplified models than monojet searches Bauer:2017fsw ; Zhou:2013fla ;
Brennan:2016xjh .

For a given parameter point of a simplified model, the the key theory input to
the monojet search likelihood for an LHC experiment is the set of predicted
event yields in each bin of the missing transverse energy distribution. In
each bin, the yield is given by

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

where $L$ is the integrated luminosity, $\sigma$ is the total production
cross-section, and $\epsilon A$ is the product of the efficiency and
acceptance for passing the kinematic selections that define the analysis.

The quantity $\epsilon A$ can be obtained by combining a Monte Carlo
simulation of DM production with a simulation of the ATLAS/CMS detector. The
standard approach to this in GAMBIT is to run the Pythia Monte Carlo generator
at each point in the global fit to simulate events, followed by a fast
detector simulation based on four-vector smearing with typical resolution
functions ColliderBit . The problem for monojet events, however, is that the
expected signatures crucially depend on the ISR model used to simulated jet
radiation in order to correctly predict the transverse missing energy
distribution Buckley:2014fba . We have therefore performed a more detailed
simulation using MadGraph_aMC@NLO Alwall:2011uj  (v3.1.1), interfaced to
Pythia v8.3 Sjostrand:2007gs  for parton showering and hadronization. We use
the CKKW prescription to perform the matching between MadGraph and Pythia. We
computed matrix elements for MadGraph starting from Universal FeynRules Output
(UFO) files Degrande:2011ua , generated with FeynRules Alloul:2013bka  and
employing a 5-flavour scheme. We used MadAnalysis 5 Conte:2012fm  to perform
detector simulation and implement each of the ATLAS and CMS monojet analyses
in order to compute $\epsilon A$. As this set of simulations is too
computationally expensive to run during the global fit, we precomputed grids
of the cross sections ($\sigma$) and $\epsilon A$ factors for each LHC
experiment in advance, and interpolated them at runtime using ColliderBit in
order to obtain predicted LHC signal yields. We then fed the predicted yields
to the likelihood functions contained in ColliderBit in order to obtain the
final constraints.

Our interpolation grids were defined as follows:

  * •

mediator mass: 20 values, 50 GeV–10 TeV

  * •

DM/mediator mass ratio: 31 values, 0.1–40

  * •

quark-mediator coupling: 6 values, 0.01–1.0

  * •

DM-mediator couplings: 7 values (each), 0.01–3.0

We used the ratio of DM and mediator masses, rather than the masses
themselves, so as to be able to include a higher density of points across the
resonance region, where rapid changes in signal prediction occur. In order to
avoid simulating points unnecessarily, we did not simulate parameter
combinations that violate Eq. 6. We limit the DM mass/mediator mass ratio to
0.1, as below this one can safely extrapolate to small DM masses. After
imposing the unitarity requirement on the grid points, the total numbers of
points are 26040 (Scalar DM), 8880 (Majorana DM) and 62160 (Dirac DM).

The CMS analysis that we include has 66 exclusive signal regions. The
Collaboration have published a covariance matrix that allows all of these to
be used simultaneously in constructing the likelihood function. We use the
“simplified likelihood” method Collaboration:2242860 , which approximates the
the full experimental likelihood function by a Poisson counting term convolved
with a multi-dimensional Gaussian likelihood describing the correlated
systematic uncertainties on the background predictions:

$\displaystyle\begin{split}\mathcal{L}_{\text{CMS}}(\bm{s},\bm{\gamma})=&\prod_{i=1}^{66}\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}$
(22)

For each signal region $i$, the observed counts, expected signal and expected
background yields are represented by $n_{i}$, $s_{i}$ and $b_{i}$,
respectively. The term $\gamma_{i}$ quantifies the deviation from the nominal
expected yields due to systematic uncertainties in signal region $i$. The set
$\\{\gamma_{i}\\}$ thus gives 66 nuisance parameters in total for this
particular analysis. The covariance matrix $\bm{\Sigma}$ provided by CMS
encodes the correlations between the various $\gamma_{i}$ factors. We
supplement these by adding the signal yield uncertainties in quadrature along
the diagonal. For every parameter point, we profile out the 66 $\gamma_{i}$
parameters so that the final CMS likelihood is defined solely in terms of the
simplified model signal estimates $\bm{s}$:

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

where $\hat{\hat{\bm{\gamma}}}$ denotes the combination of background nuisance
parameters resulting in the highest value of the likelihood for the given
signal $\bm{s}$.

The ATLAS analysis does not come with a published covariance matrix, nor with
a published likelihood in the HistFactory format of Ref. ATL-PHYS-PUB-2019-029
. The conservative course of action is therefore to calculate a likelihood
using the signal region with the best expected sensitivity. To maximize the
sensitivity of this procedure, we combine the three highest missing energy
bins so that $\not{E}_{\mathrm{T}}>1000\,\mathrm{GeV}$ is the highest bin in
the analysis. This is justified by the fact that the systematic uncertainties
in the background estimations for these bins (and hence their correlations)
are negligible. The ATLAS likelihood is then given by

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

where $\mathcal{L}_{\text{ATLAS}}(s_{i},\hat{\hat{\gamma_{i}}})$ is the
single-bin equivalent of Eq. (22), whilst $i$ labels the signal region that
would give the lowest likelihood in the case $n_{i}=b_{i}$.

Assuming that the ATLAS and CMS searches are independent, we then calculate a
total log-likelihood for LHC monojet searches as
$\ln\mathcal{L}_{\text{LHC}}=\ln\mathcal{L}_{\text{CMS}}+\ln\mathcal{L}_{\text{ATLAS}}$.

Note, however, that the choice of different ATLAS signal regions for different
parameter points leads to a large variation in the effective likelihood
normalisation between different parameter points. The standard ColliderBit
solution is to instead set the LHC log-likelihood contribution to the total
difference in log-likelihood between the signal and background-only
($\bm{s}=\bm{0}$) cases:

$\displaystyle\Delta\ln\mathcal{L}_{\text{LHC}}=\ln\mathcal{L}_{\text{LHC}}(\bm{s})-\ln\mathcal{L}_{\text{LHC}}(\bm{s}=\bm{0}).$
(25)

Positive values of this quantity indicate that a DM model fits the observed
data better than the background-only hypothesis. In cases where there are
small excesses in the LHC data, this can lead to regions of the simplified
model parameter space fitting the data better than other regions that might be
indistinguishable from the SM. As our global fit results are presented as
$1\sigma$ and $2\sigma$ confidence regions defined using the likelihood ratio
$\mathcal{L}(\bf{\Theta})/\mathcal{L}(\bf{\Theta}_{\text{best-fit}})$ around
the best-fit parameters $\bf{\Theta}_{\text{best-fit}}$, this can exclude
parameter regions that exhibit only a little worse agreement with the data
than the SM. Whilst this is of course the correct result in the case where one
takes excesses at face value and attempts to fit them, it can also be
instructive to consider LHC results under the conservative assumption that any
excesses simply arise from statistical fluctuations rather than BSM physics.
We therefore run separate scans for this latter case, where we “cap” the LHC
likelihood as

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

More detail on this procedure can be found in Ref. EWMSSM .

###  3.5 Collider searches for the mediator using dijet events

Because our simplified models explicitly include a mediator particle that
interacts with SM quarks, it is possible for the LHC to produce a mediator
that decays to quarks, rather than WIMPs. This should generate an excess of
dijet events, each with a dijet invariant mass approximately equal to the mass
of the mediator. Searches for dijet resonances thus provide powerful
constraints on DM simplified models, though analyses have to employ various
clever tricks to increase sensitivity given the extremely large multijet
background at the LHC.

Assuming that a narrow width approximation holds,
$\sigma(\mathrm{pp}\rightarrow
V_{\mu}\rightarrow\mathrm{qq})\approx\sigma(\mathrm{pp}\rightarrow
V_{\mu})\times BR(V_{\mu}\rightarrow\mathrm{qq})$. In this case,
$\sigma\propto g_{\mathrm{q}}^{2}$, so we implement ATLAS and CMS dijet limits
by appropriately scaling the published limits by the mediator-quark coupling
and the branching ratio into quarks, following the same approach as Ref.
Bagnaschi:2019djj . We interpolate these published limits in $m_{\mathrm{M}}$
during our scans and select the most constraining search for a given mediator
mass. This way, we are able to recast the published limits without using Monte
Carlo simulation.

We then compare our results to the coupling upper limits
($g_{\mathrm{q,excl}}$) provided by a broad range of LHC dijet searches
CMS:2019gwf ; ATLAS:2019fgd ; ATLAS:2018qto ; CDF:2008ieg ; ATLAS:2018hbc ;
ATLAS:2018hzj ; CMS:2019emo ; ATLAS:2019itm ; CMS:2019xai , using a likelihood
of the form

$\displaystyle\ln{\mathcal{L}_{\mathrm{dijet}}}=-2\Big{(}\frac{g_{\mathrm{q,pred}}^{4}\times
BR(V_{\mu}\rightarrow\mathrm{qq})^{2}}{g_{\mathrm{q,excl}}^{4}}\Big{)}\,.$
(27)

The factor of $-2$ arises because $g_{\mathrm{q,excl}}$ is taken from the
$95\%$ confidence limit on the coupling provided by each analysis,
corresponding to $\Delta\chi^{2}=4$. This produces coupling upper limits as
shown in Figure 1.

Figure 1: Quark coupling upper limits from each dijet search included in our
likelihood Eq. 27. For a given mediator mass, the $95\%$ confidence dijet
limit that we use in our likelihood is the one that is the most constraining
(i.e. closest to the bottom of this plot).

There is a high degree of variation with mediator mass in the limits presented
from each analysis, due to the signal prediction moving in and out of
different experimental analysis bins. The effect of these variations will not
show up strongly in our later results, as we profile over the model couplings
and the highest likelihood points will tend to be where dijet constraints are
unconstraining.

###  3.6 Nuisance parameter likelihoods

The model parameters for each of our simplified models are supplemented by a
series of nuisance parameters, which contribute to our astrophysical
likelihoods. We give a complete list of nuisance parameters in Table 3.

Our treatment of the local DM density $\rho_{0}$ follows the default
prescription in DarkBit, which assumes that $\rho_{0}$ is log-normally
distributed with central value $\rho_{0}=0.40$ GeV cm-3 and error
$\sigma_{\rho_{0}}=0.15$ GeV cm-3. The asymmetric scan range for $\rho_{0}$ in
Table 3 reflects this log-normal distribution. All other nuisance parameters
are scanned over their $3\sigma$ range as provided in Table 3.

We treat the Milky Way halo in the same way as in our previous studies of
Higgs portal and dark matter effective field theories HP ; DMEFT . The DM
velocity follows a Maxwell-Boltzmann distribution, with uncertainties on the
peak of the distribution and the Galactic escape velocity described by
Gaussian likelihoods with $v_{\rm{peak}}=240\,\pm\,8$ km s-1 Reid:2014boa  and
$v_{\rm{esc}}=528\pm 25$ km s-1 (based on _Gaia_ data Deason:2019kgj ),
respectively.

The expected yields in direct DM detection experiments are strongly affected
by our knowledge of various hadronic quantities needed for calculating form
factors at the nuclear scale. We use a Gaussian likelihood for the quantity
$r_{\rm s}^{2}$ Djukanovic:2019jtp ; Sufian:2016pex , using the mean and
standard deviation listed in Table 3. The fixed values of other hadronic
quantities assumed by DirectDM v2.2.0 can be found in Ref. DMEFT .

##  4 Results

We have performed comprehensive scans of each simplified model parameter space
using the differential evolution sampler Diver v1.0.4 ScannerBit  with a
population of $10\,000$ and a convergence threshold of $10^{-6}$. We present
results as profile likelihood maps in planes of the parameters and/or derived
quantities. We carried out scans for the cases where the observed DM relic
density is taken as an upper limit or as a two-sided measurement, and also for
the cases where the LHC likelihood is capped or uncapped. This gives four
combinations of scans for each model. The parameter points with the highest
likelihoods for each model are given in Table 4.

Table 3 outlines our complete list of parameters and their associated scan
ranges. We do not consider models with large fine-tuning or large hierarchies
between the different couplings, as this may be challenging to achieve when
considering plausible UV embeddings of the simplified models. In order to
equally explore each order of magnitude in the Lagrangian parameters, we
sampled them using log-uniform priors. For nuisance parameters, we sampled
using flat priors. We stress, however, that these ‘prior’ distributions play
no formal role in the final statistical analysis of the profile likelihood
maps that we present, but merely provide efficiently distributed starting
guesses from which to hunt for best-fit points and map likelihood functions.

Parameters | Range
---|---
DM mass, $m_{\mathrm{DM}}$ |  $[50,10000]$ GeV
Mediator mass, $m_{\mathrm{M}}$ |  $[50,10000]$ GeV
quark-mediator coupling, $g_{\mathrm{q}}$ | $[0.01,1.0]$
mediator-DM coupling (vector), $g^{\rm V}_{\mathrm{DM}}$ | $[0.01,3.0]$
mediator-DM coupling (axial vector), $g^{\rm A}_{\mathrm{DM}}$ | $[0.01,3.0]$
Nuisance Parameters |  Value ($\pm 3\sigma$ range)
strange quark proton charge radius, $r^{2}_{\rm s}$ |  $-0.115(0.105)$ GeV-2
Local DM density, $\rho_{0}$ |  $[0.2,0.8]$ GeV cm-3
Most probable speed, $v_{\mathrm{peak}}$ |  $240(24)$ km s-1
Galactic escape speed, $v_{\mathrm{esc}}$ |  $528(75)$ km s-1
Table 3: Ranges scanned over for model and nuisance parameters. The axial-vector coupling is present only in the Dirac fermion model. Hadronic input parameters are given at $\mu=2$ GeV. LHC | Relic Density | Best Fit $m_{\mathrm{DM}}$ | Best Fit $m_{\mathrm{M}}$ | Best Fit $g_{\mathrm{q}}$ | Best Fit $g^{\rm V}_{\mathrm{DM}}$ | Best Fit $g^{\rm A}_{\mathrm{DM}}$ | $\Delta\ln\mathcal{L}$
---|---|---|---|---|---|---|---
Scalar DM |  |  |  |  |  |
Full | Upper limit | 4965 | 10000 | 0.0100 | 2.333 | – | -0.019
Full | All DM | 4532 | 9203 | 0.0101 | 0.825 | – | -0.469
Dirac fermion DM |  |  |  |  |  |
Full | Upper limit | 262 | 537 | 0.0301 | 0.0100 | 0.990 | 4.48
Capped | Upper limit | 146 | 300 | 0.0108 | 0.0103 | 2.525 | -0.089
Full | All DM | 588 | 1320 | 0.0369 | 0.0100 | 0.754 | 0.881
Capped | All DM | 3762 | 7744 | 0.0151 | 0.0118 | 0.536 | -0.559
Majorana fermion DM |  |  |  |  |  |
Full | Upper limit | 50 | 114 | 0.0130 | – | 0.243 | 4.779
Capped | Upper limit | 668 | 1400 | 0.0139 | – | 0.763 | 0.001
Full | All DM | 69.1 | 166 | 0.0104 | – | 0.307 | 3.12
Capped | All DM | 2423 | 5079 | 0.0279 | – | 0.125 | -0.449
Table 4: Approximate best-fit points for each model. $\Delta\ln\mathcal{L}$
values are defined as $\ln\mathcal{L}-\ln\mathcal{L}^{\mathrm{bg}}$, where the
background-only likelihood is detailed in Table 1.

###  4.1 Scalar DM

Figure 2: Scalar DM profile likelihood, profiling over couplings. The observed
relic density of dark matter is taken as an upper limit (top) or to consist
entirely of the scalar DM candidate (bottom). $1\sigma$ and $2\sigma$ contours
are shown in white, with the star representing the best-fit point.

Figure 3: Relic Density of the scalar DM model as a function of the mediator
mass (left) and DM mass (right). The observed relic density of dark matter is
taken as an upper limit on the abundance of the scalar DM candidate $1\sigma$
and $2\sigma$ contours are shown in white, with the star representing the
best-fit point.

Results from global scans of the scalar DM model are shown in Figure 2. Any
excesses present in the mono-jet likelihoods can only be fitted by models that
are already robustly excluded by other searches, so we do not show any results
for this model based on the capped collider likelihood.

The results show two separate regions allowed at the $2\sigma$ level. The
shape of the surviving parameter space is defined primarily by the relic
abundance. The viable parameter space is split in two corresponding to the two
DM annihilation channels, along the resonance for s-channel annihlation into a
pair of quarks, and $m_{\mathrm{DM}}>m_{\mathrm{M}}$, where $t$-channel
annihilation into a pair of mediators is possible. Highest likelihood points
lie along the resonance, where $2m_{\mathrm{DM}}\approx m_{\mathrm{M}}$ and
annihilation is enhanced, bringing the DM relic density down to or below the
observed value. This region is most preferred toward the upper mass limits of
the scan, where the best-fit point can be seen indicated by white star in
Figure 2.

When requiring that the scalar DM candidate explains all of DM (Figure 2),
mediator masses up to approximately $2$ TeV are excluded. This is because for
lower mediator masses, there are no allowed regions where the relic abundance
is saturated. This can be seen clearly in Figure 3, where we plot the relic
density as a function of the mediator mass (left) and DM mass (right), in the
scan where the relic density is taken as an upper limit. We can see here that
relic abundance estimates produce strong limits on the dark matter mass.
Whether the relic density is taken as an upper limit or a two-sided
constraint, DM masses below approximately $1$ TeV are excluded. This exclusion
is, of course, dependent on the limits of the quark coupling $g_{\rm q}$ in
the scan, and reducing this lower bound will expand the allowed region.

Indirect detection limits give a slight preference to the region along the
resonance. This is because when $m_{\mathrm{DM}}>m_{\mathrm{M}}$, $t$-channel
annihilation of DM to mediator particles (and subsequent decays to SM
products) would produce an observable effect on gamma-ray searches. This
signal would be in weak tension with the absence of a positive detection thus
far. This effect is reduced toward the lower mediator masses, as gamma ray
predictions are scaled by the relic abundance which significantly
underpredicts the DM abundance for lower $m_{\mathrm{M}}$.

Direct detection limits give the lower bound on DM masses, both along the
resonance region and when $m_{\mathrm{DM}}>m_{\mathrm{M}}$. These experiments
also drive the best-fit point towards the border of the scan, where the
predicted signal decreases. Near the boundaries of the scan, this likelihood
estimate is close enough to the zero signal likelihood that the magnitude of
the profile likelihood ratios would not be noticeably changed by extending the
scan limits to higher mediator masses.

Dijet searches contribute by giving preference to models with lower quark
couplings, but monojet searches do not have a strong influence on the overall
profile likelihood at all.

Given that the likelihood is weakly dependent on the couplings, and dependent
primarily on the ratio of the mediator and DM masses, a reasonable estimate
for the number of effective degrees of freedom would be 1 or 2. We compute an
approximate $p$-value of the best-fit likelihood against the background only
scenario described in Table 1 whenever the background only scenario is
preferred. At the best-fit point, for 1–2 effective degrees of freedom, the
$p$-value is between $0.85$ and $0.98$ when allowing the scalar DM to
underpredict the relic abundance, and between $0.33$ and $0.63$ when
saturating the abundance. Neither of these are statistically distinct from the
Standard Model.

###  4.2 Dirac Fermion DM

Figure 4: Profile likelihood for the Dirac fermion DM model. The observed
relic density of dark matter is taken as an upper limit (top) or to consist
entirely of the Dirac DM candidate (bottom). The collider likelihood is either
uncapped (left), or capped to prevent preference over the Standard Model
(right). $1\sigma$ and $2\sigma$ contours are shown in white, with the star
representing the best-fit point.

The constraints on the Dirac fermion DM model are shown in Figure 4. As with
the scalar DM, a large portion of the non-excluded parameter space lies on the
resonance region where $2m_{\mathrm{DM}}\approx m_{\mathrm{M}}$.

Figure 5: Missing energy spectra for the CMS monojet search. The SM background
prediction (purple) and the observed event counts (black) are taken from Ref
CMS:2021snz . The green distribution shows the signal prediction for the best-
fit Dirac DM model. The bottom panel shows the residuals, defined as (data -
prediction)/uncertainty for both the SM and the SM + DM predictions.

Excesses in the monojet collider searches are partially fit by the model. This
leaves a highly-preferred region around a mediator mass of $500$ GeV (Figure
4, top left). In Figure 5, we show the signal for the best-fit point that
would contribute to the CMS monojet search CMS:2021snz . The signal regions in
this search are split into three years of data taking, with a strong
difference between SM prediction and data in the last year. This difference is
present in the simplified likelihood only, but was not present in their joint
fit of control and signal regions. The model cannot entirely fit the 2018
excess without strongly over predicting the signal in the two other years. If
these excesses are assumed to simply be statistical fluctuations, and the
monojet likelihood is capped (Figure 4, right), then the surviving parameter
space opens up.

For the sake of illustration, we also show the bounds from unitarity that
apply to this model in the top right panel of Figure 4. There is no strong
preference toward the unitarity bound as the profile likelihood in this region
is in fact very flat in $g^{\rm A}_{\mathrm{DM}}$ — a fact reflected in the
range of best-fit couplings shown in Table 4.

Relic density limits have a strong influence on the exclusion contours, such
that the allowed regions lie either on the resonance or, in the case where
collider likelihoods are capped, where $m_{\mathrm{DM}}>m_{\mathrm{M}}$.
Requiring that the DM candidate constitute all of the observed DM excludes
lower mass regions along the resonance. This shifts the best-fit point to
higher masses around DM mass of $1300$,GeV, and reduces the ability of the
best fit to fit the monojet excesses. This suppression of the likelihood of
the best fit (compared to that found when imposing the relic density as an
upper limit) opens up some of the parameter space off resonance at the very
limits of the 2$\sigma$ region. For the results where we cap the collider
likelihood, the best-fit likelihood is also reduced by requiring the DM
abundance to be saturated, which broadens the surviving parameter space along
the resonance region. Figure 6 shows that if the monojet excesses were
explained by this model, rather than by statistical fluctuations in the
experiments, this particle would not saturate the DM relic abundance.
Additional DM candidates would be required to explain the observed relic
abundance.

The combination of direct detection, relic abundance and unitarity constraints
provide the shape of the off-resonance region for the capped collider results.
To avoid unitarity violation, the model is excluded for larger $g^{\rm
A}_{\mathrm{DM}}$. This in turn may prevent sufficient annihilation and cause
the predicted relic abundance to exceed the observed value. Direct detection
experiments provide a lower bound on the mediator mass for a given DM mass and
$g^{\rm V}_{\mathrm{DM}}$. Since the strength of the direct detection signals
is primarily from the $g^{\rm V}_{\mathrm{DM}}$ coupling, and the unitarity
bound is from the $g^{\rm A}_{\mathrm{DM}}$ coupling, this shape in parameter
space would differ in either the pure vector or pure axial-vector coupling
cases (as studied in Ref. Bagnaschi:2019djj ). If the $g_{\rm q}$ limit was
lowered, this allowed region would expand such that the two off-resonance
regions would become one. Indirect detection searches do not play a strong
role in the overall exclusion contours for this model.

Since the model is preferred over the Standard Model when the collider
likelihood is uncapped, we only calculate a $p$-value for the capped collider
scans. At the best-fit point, for 1–2 effective degrees of freedom, the
$p$-value is between $0.67$ and $0.91$ when allowing the Dirac fermion to
underpredict the relic abundance, and between $0.29$ and $0.57$ when
saturating the abundance. These are not statistically distinct from the
Standard Model.

For the narrow width approximation to hold, the ratio of mediator decay width
to mediator mass must remain low. Figure 7 shows that this ratio increases
with higher mediator masses and lower DM masses. In the surviving parameter
space of the scan, this ratio can reach at most roughly 0.4–0.5. As this is
close to the limit that would prevent accurate recasting of dijet search
limits, doubt could be raised about the validity of applying these limits.
However, this occurs in regions where dijet limits are unconstraining. In all
regions where collider limits contribute noticeably, the model safely
satisfies the narrow width approximation.

The results differ from those in Ref. Bagnaschi:2019djj  as we present
combined fits of all 5 model parameters varying concurrently, whereas they
separate the model into pure vector/axial-vector cases. We also allow the
model to fit monojet excesses and give an overall preference over the Standard
Model, where they do not. In this way, this study is complementary to
Bagnaschi:2019djj  without presenting duplication of their results.

Figure 6: Profile likelihood as a function of relic density of the Dirac
fermion DM model, taking the relic abundance as an upper limit. The capped
collider result is shown on the right. We do not show the dependence on
mediator mass as it does not differ greatly from the dependence on DM mass.
$1\sigma$ and $2\sigma$ contours are shown in white, with the star
representing the best-fit point. Figure 7: Ratio of mediator decay width to
mass for all regions allowed in Figure 4 (top right).

###  4.3 Majorana Fermion DM

Figure 8: Profile likelihood for the Majorana fermion DM model. The observed
relic density of dark matter is taken as an upper limit (top) or a two-sided
measurement that the model must match (bottom). The collider likelihood is
either uncapped (left), or capped to prevent preference over the Standard
Model (right). $1\sigma$ and $2\sigma$ contours are shown in white, with the
star representing the best-fit point.

Figure 9: Relic density of the Majorana fermion DM model, taking the relic
abundance as an upper limit. Capped collider results are shown on the right.
As with the Dirac model, we do not show the dependence on mediator mass as it
does not differ greatly from the dependence on DM mass. $1\sigma$ and
$2\sigma$ contours are shown in white, with the star representing the best-fit
point.

Results from global scans of the Majorana fermion DM model are shown in Figure
8. Like the Dirac fermion model, there is a strong preference over the
background along the resonance. The collider excess is predominantly fit along
the resonance region.

Figure 8 (right) shows how capping the collider likelihood expands the
$2\sigma$ contours, with little exclusion when
$m_{\mathrm{DM}}>m_{\mathrm{M}}$. The reason this region is much larger in the
Majorana models compared to the Dirac model, is that the direct detection
signals from axial-vector couplings are suppressed in the non-relativistic
limit. The presence of vector couplings in the Dirac model give rise to strong
enough direct detection signals such that this region is not allowed for the
scanned coupling ranges, i.e. even with $g^{\rm V}_{\mathrm{DM}}=0.01$.

As with the Dirac fermion and scalar models, relic density limits play a
strong role in determining the shape of the allowed regions. Figure 9 shows
the relic abundance as a function of the DM mass. The best-fit point appears
to under-predict the relic abundance, although there is little preference for
the best-fit point over a region that saturates the abundance. Requiring that
the DM candidate saturates the relic abundance, in Figure 8 (bottom), shifts
the location of the best-fit point slightly towards higher mediator/DM masses,
however the capped collider results are very similar regardless of whether the
abundance is saturated or not. The approximate $p$-value for 1–2 effective
degrees of freedom is between $0.34$ and $0.64$ in the case of a saturated
relic abundance and a capped collider likelihood. The best-fit points for the
other three scans gave preferences over the Standard Model.

The preference for the resonance regions over the
$m_{\mathrm{DM}}>m_{\mathrm{M}}$ region in the fits with capped
$\mathcal{L}_{\text{LHC}}$ is because in the latter region of parameter space,
the gamma-ray flux is not negligible when the annihilation channel into
mediators is open. This increase in the annihilation cross-section at late
times is not matched by a drop in the relic density fraction, as $p$-wave
annihilation dominates at early times. The result is that the _Fermi_ -LAT
likelihood gives a preference to regions where
$m_{\mathrm{DM}}<m_{\mathrm{M}}$.

###  4.4 Future prospects

Figure 10: Allowed parameter space as a function of the DM mass and rescaled
annihilation cross section. Capped collider likelihoods are shown for Dirac
(middle) and Majorana (bottom) DM as these are the least constraining. The
non-capped likelihood is shown for the scalar DM model (top). The solid red
line indicates the projected “initial construction” sensitivity of the
Cherenkov Telescope Array (CTA) towards the Galactic Centre (GC).

Figure 11: Mediator mass as a function of the predicted number of signal
events in the DARWIN experiment, when the DM candidate accounts for all of the
observed DM abundance. Capped collider likelihoods are shown for Dirac
(middle) and Majorana (bottom) DM as these are the least constraining. The
non-capped likelihood is shown for the scalar DM model (top).

The upcoming Cherenkov Telescope Array (CTA) has the potential to probe the
surviving parameter space of these models. Figure 10 shows the annihilation
cross section, rescaled by the DM fraction (see section 3.1 for details),
overlayed with projections of the sensitivity after initial construction for
the CTA experiment Acharyya:2020sbj . The scalar DM model may only be affected
at the upper limits of DM masses in the scan, above $5$ TeV. The strongest
projected limits will occur for the Dirac fermion model. CTA would seem to
have minimal impact on the Majorana model. All of the parameter space that
would be constrained corresponds to the off-resonance regions, where
$m_{\mathrm{DM}}>m_{\mathrm{M}}$.

We show the predicted number of signal events in DARWIN, a next-generation
direct detection experiment DARWIN , as a function of the DM mass in Figure
11. We restrict these to the case where the DM relic abundance is saturated.
Each model has the potential to produce tens or even more than a hundred
events. DARWIN will therefore prove useful for constraining for all the models
that we consider here. This should strongly constrain both the scalar and
Dirac fermion DM models, including much of the off-resonance regions. For the
Majorana model, the DARWIN experiment may impact the lowest DM masses, but
much of the high mediator mass regions would remain untouched.

##  5 Conclusions

In this work we performed global scans with GAMBIT of three simplified dark
matter models that interact with quarks via an $s$-channel vector mediator. We
found that in the case of scalar DM, a great deal of parameter space with high
DM mass and low mediator mass survives, provided that it does not constitute
all of the observed relic abundance of DM. Requiring that the model saturates
the relic abundance restricts the parameter space to lie either towards the
upper mass bounds in this study, or along the region with resonant enhancement
in the cross-section. Global constraints on Dirac fermion DM are driven by
small excesses in monojet searches, which prefer dark matter masses of around
260 GeV and mediator masses of around 540 GeV. In these regions, the candidate
could only constitute a subcomponent of the total relic abundance; requiring
that it saturates the DM abundance raises the preferred region to mediator
masses of roughly $1.3$ TeV and DM masses of $580$ GeV. Assuming that the
collider excesses are only background fluctuations, we also perform scans
where the LHC likelihood can only exclude the parameter regions that fit the
data worse than the SM. We then see the parameter space open up much more,
allowing regions with non-resonant production to survive. Similarly Majorana
fermion DM also fits these monojet excesses, and can do so whilst saturating
the DM relic abundance. For all of the simplified models that we studied,
direct detection, dijet searches and relic abundance measurements provide
complementary limits. A lower bound on the DM mass was found for the scalar
and Dirac DM models from competing constraints from relic abundance and
laboratory experiments, as was found for the EFT in Ref. GAMBIT:2021rlp . This
bound is strongly dependent on the scan limits on the couplings. We showed
with projections for the CTA and DARWIN experiments that these will be capable
of exploring some of the remaining parameter space of the scalar and Dirac
fermion models. Also, improvements to dijet limits from Run 3 of the LHC
should lower the upper bound on mediator-quark couplings, which will further
restrict the surviving parameter space when combined with other constraints.

##  6 Acknowledgements

We thank PRACE for awarding us access to Joliot-Curie at CEA. This project was
also undertaken with assistance of resources and services from National
Computational Infrastructure, which is supported by the Australian Government.
This work was in part performed using the Cambridge Service for Data Driven
Discovery (CSD3), part of which is operated by the University of Cambridge
Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk).
The DiRAC component of CSD3 was funded by BEIS capital funding via STFC
capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant
ST/R00689X/1. DiRAC is part of the National e-Infrastructure. PS acknowledges
funding support from the Australian Research Council under Future Fellowship
FT190100814. TEG and FK were funded by the Deutsche Forschungsgemeinschaft
(DFG) through the Emmy Noether Grant No. KA 4662/1-1 and grant 396021762 - TRR
257. This article made use of pippi v2.2 pippi . MJW is supported by the ARC
Centre of Excellence for Dark Matter Particle Physics (CE200100008).

## References

  * (1) F. Zwicky, Die Rotverschiebung von extragalaktischen Nebeln, Helvetica Physica Acta 6 (1933) 110–127.
  * (2) D. Clowe, M. Bradač, et. al., A Direct Empirical Proof of the Existence of Dark Matter, ApJ 648 (2006) L109–L113, [astro-ph/0608407].
  * (3) D. N. Spergel, R. Bean, et. al., Three-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Implications for Cosmology, ApJS 170 (2007) 377–408, [astro-ph/0603449].
  * (4) B. W. Lee and S. Weinberg, Cosmological Lower Bound on Heavy Neutrino Masses, Phys. Rev. Lett. 39 (1977) 165–168.
  * (5) G. Arcadi, M. Dutra, et. al., The waning of the WIMP? A review of models, searches, and constraints, Eur. Phys. J. C78 (2018) 203, [arXiv:1703.07364].
  * (6) J. Goodman, M. Ibe, et. al., Gamma Ray Line Constraints on Effective Theories of Dark Matter, Nucl. Phys. B 844 (2011) 55–68, [arXiv:1009.0008].
  * (7) M. Beltran, D. Hooper, E. W. Kolb, and Z. C. Krusberg, Deducing the nature of dark matter from direct and indirect detection experiments in the absence of collider signatures of new physics, Phys. Rev. D 80 (2009) 043509, [arXiv:0808.3384].
  * (8) K. Cheung, P.-Y. Tseng, and T.-C. Yuan, Gamma-ray Constraints on Effective Interactions of the Dark Matter, JCAP 06 (2011) 023, [arXiv:1104.5329].
  * (9) R. Harnik and G. D. Kribs, An Effective Theory of Dirac Dark Matter, Phys. Rev. D 79 (2009) 095007, [arXiv:0810.5557].
  * (10) A. De Simone, A. Monin, A. Thamm, and A. Urbano, On the effective operators for Dark Matter annihilations, JCAP 02 (2013) 039, [arXiv:1301.1486].
  * (11) C. Karwin, S. Murgia, T. M. P. Tait, T. A. Porter, and P. Tanedo, Dark Matter Interpretation of the Fermi-LAT Observation Toward the Galactic Center, Phys. Rev. D 95 (2017) 103005, [arXiv:1612.05687].
  * (12) J. Fan, M. Reece, and L.-T. Wang, Non-relativistic effective theory of dark matter direct detection, JCAP 1011 (2010) 042, [arXiv:1008.1591].
  * (13) P. Agrawal, Z. Chacko, C. Kilic, and R. K. Mishra, A Classification of Dark Matter Candidates with Primarily Spin-Dependent Interactions with Matter, arXiv:1003.1912.
  * (14) A. Crivellin and U. Haisch, Dark matter direct detection constraints from gauge bosons loops, Phys. Rev. D 90 (2014) 115011, [arXiv:1408.5046].
  * (15) M. Hoferichter, P. Klos, J. Menéndez, and A. Schwenk, Analysis strategies for general spin-independent WIMP-nucleus scattering, Phys. Rev. D 94 (2016) 063505, [arXiv:1605.08043].
  * (16) F. Kahlhoefer and S. Wild, Studying generalised dark matter interactions with extended halo-independent methods, JCAP 10 (2016) 032, [arXiv:1607.04418].
  * (17) O. Buchmueller, M. J. Dolan, and C. McCabe, Beyond Effective Field Theory for Dark Matter Searches at the LHC, JHEP 01 (2014) 025, [arXiv:1308.6799].
  * (18) J. Abdallah et. al., Simplified Models for Dark Matter Searches at the LHC, Phys. Dark Univ. 9-10 (2015) 8–23, [arXiv:1506.03116].
  * (19) J. Abdallah et. al., Simplified Models for Dark Matter and Missing Energy Searches at the LHC, arXiv:1409.2893.
  * (20) S. A. Malik et. al., Interplay and Characterization of Dark Matter Searches at Colliders and in Direct Detection Experiments, Phys. Dark Univ. 9-10 (2015) 51–58, [arXiv:1409.4075].
  * (21) G. Busoni, A. De Simone, E. Morgante, and A. Riotto, On the Validity of the Effective Field Theory for Dark Matter Searches at the LHC, Phys. Lett. B 728 (2014) 412–421, [arXiv:1307.2253].
  * (22) G. Busoni, A. De Simone, J. Gramling, E. Morgante, and A. Riotto, On the Validity of the Effective Field Theory for Dark Matter Searches at the LHC, Part II: Complete Analysis for the $s$-channel, JCAP 06 (2014) 060, [arXiv:1402.1275].
  * (23) G. Busoni, A. De Simone, T. Jacques, E. Morgante, and A. Riotto, On the Validity of the Effective Field Theory for Dark Matter Searches at the LHC Part III: Analysis for the $t$-channel, JCAP 09 (2014) 022, [arXiv:1405.3101].
  * (24) GAMBIT: P. Athron et. al., Thermal WIMPs and the scale of new physics: global fits of Dirac dark matter effective field theories, Eur. Phys. J. C 81 (2021) 992, [arXiv:2106.02056].
  * (25) C. Arina, Impact of cosmological and astrophysical constraints on dark matter simplified models, Front. Astron. Space Sci. 5 (2018) 30, [arXiv:1805.04290].
  * (26) A. De Simone and T. Jacques, Simplified models vs. effective field theory approaches in dark matter searches, Eur. Phys. J. C 76 (2016) 367, [arXiv:1603.08002].
  * (27) A. Albert et. al., Recommendations of the LHC Dark Matter Working Group: Comparing LHC searches for dark matter mediators in visible and invisible decay channels and calculations of the thermal relic density, Phys. Dark Univ. 26 (2019) 100377, [arXiv:1703.05703].
  * (28) A. Boveia et. al., Recommendations on presenting LHC searches for missing transverse energy signals using simplified $s$-channel models of dark matter, Phys. Dark Univ. 27 (2020) 100365, [arXiv:1603.04156].
  * (29) F. Kahlhoefer, Review of LHC Dark Matter Searches, Int. J. Mod. Phys. A 32 (2017) 1730006, [arXiv:1702.02430].
  * (30) E. Morgante, Simplified Dark Matter Models, Adv. High Energy Phys. 2018 (2018) 5012043, [arXiv:1804.01245].
  * (31) F. D’Eramo, B. J. Kavanagh, and P. Panci, You can hide but you have to run: direct detection with vector mediators, JHEP 08 (2016) 111, [arXiv:1605.04917].
  * (32) L. M. Carpenter, R. Colburn, J. Goodman, and T. Linden, Indirect Detection Constraints on s and t Channel Simplified Models of Dark Matter, Phys. Rev. D 94 (2016) 055027, [arXiv:1606.04138].
  * (33) D. Abercrombie et. al., Dark Matter Benchmark Models for Early LHC Run-2 Searches: Report of the ATLAS/CMS Dark Matter Forum, Phys. Dark Univ. 26 (2019) 100371, [arXiv:1507.00966].
  * (34) E. Bagnaschi et. al., Global Analysis of Dark Matter Simplified Models with Leptophobic Spin-One Mediators using MasterCode, Eur. Phys. J. C 79 (2019) 895, [arXiv:1905.00892].
  * (35) GAMBIT Collaboration: P. Athron, C. Balázs, et. al., GAMBIT: The Global and Modular Beyond-the-Standard-Model Inference Tool, Eur. Phys. J. C 77 (2017) 784, [arXiv:1705.07908]. Addendum in gambit_addendum.
  * (36) GAMBIT Collaboration: S. Bloor, T. E. Gonzalo, et. al., The GAMBIT Universal Model Machine: from Lagrangians to likelihoods, Eur. Phys. J. C 81 (2021) 1103, [arXiv:2107.00030].
  * (37) GAMBIT Collaboration, Supplementary Data: Global Fits of simplified models for dark matter with GAMBIT. I. Scalar and fermionic models with s-channel vector mediators., (2022), https://zenodo.org/record/6615830.
  * (38) C. Chang, P. Scott, et. al., Global fits of simplified models for dark matter with GAMBIT II. Vector dark matter with an $s$-channel vector mediator, in prep (2022).
  * (39) M. Duerr, F. Kahlhoefer, K. Schmidt-Hoberg, T. Schwetz, and S. Vogl, How to save the WIMP: global analysis of a dark matter model with two s-channel mediators, JHEP 09 (2016) 042, [arXiv:1606.07609].
  * (40) M. Duerr and P. Fileviez Perez, Theory for Baryon Number and Dark Matter at the LHC, Phys. Rev. D 91 (2015) 095001, [arXiv:1409.8165].
  * (41) J. Ellis, M. Fairbairn, and P. Tunney, Anomaly-Free Dark Matter Models are not so Simple, JHEP 08 (2017) 053, [arXiv:1704.03850].
  * (42) F. Kahlhoefer, K. Schmidt-Hoberg, T. Schwetz, and S. Vogl, Implications of unitarity and gauge invariance for simplified dark matter models, JHEP 02 (2016) 016, [arXiv:1510.02110].
  * (43) SuperCDMS: R. Agnese et. al., New Results from the Search for Low-Mass Weakly Interacting Massive Particles with the CDMS Low Ionization Threshold Experiment, Phys. Rev. Lett. 116 (2016) 071301, [arXiv:1509.02448].
  * (44) CRESST: G. Angloher et. al., Results on light dark matter particles with a low-threshold CRESST-II detector, Eur. Phys. J. C76 (2016) 25, [arXiv:1509.01515].
  * (45) CRESST: A. H. Abdelhameed et. al., First results from the CRESST-III low-mass dark matter program, Phys. Rev. D 100 (2019) 102002, [arXiv:1904.00498].
  * (46) DarkSide: P. Agnes et. al., DarkSide-50 532-day Dark Matter Search with Low-Radioactivity Argon, Phys. Rev. D 98 (2018) 102006, [arXiv:1802.07198].
  * (47) LUX: D. S. Akerib et. al., Results from a search for dark matter in the complete LUX exposure, Phys. Rev. Lett. 118 (2017) 021303, [arXiv:1608.07648].
  * (48) PICO: C. Amole et. al., Dark Matter Search Results from the PICO-60 C3F8 Bubble Chamber, Phys. Rev. Lett. 118 (2017) 251301, [arXiv:1702.07666].
  * (49) PICO: C. Amole et. al., Dark Matter Search Results from the Complete Exposure of the PICO-60 C3F8 Bubble Chamber, Phys. Rev. D 100 (2019) 022001, [arXiv:1902.04031].
  * (50) PandaX-II: A. Tan et. al., Dark Matter Results from First 98.7 Days of Data from the PandaX-II Experiment, Phys. Rev. Lett. 117 (2016) 121303, [arXiv:1607.07400].
  * (51) PandaX-II: X. Cui et. al., Dark Matter Results From 54-Ton-Day Exposure of PandaX-II Experiment, Phys. Rev. Lett. 119 (2017) 181302, [arXiv:1708.06917].
  * (52) XENON: E. Aprile et. al., Dark Matter Search Results from a One Ton-Year Exposure of XENON1T, Phys. Rev. Lett. 121 (2018) 111302, [arXiv:1805.12562].
  * (53) LZ: J. Aalbers et. al., First Dark Matter Search Results from the LUX-ZEPLIN (LZ) Experiment, arXiv:2207.03764.
  * (54) CMS: A. M. Sirunyan et. al., Search for high mass dijet resonances with a new background prediction method in proton-proton collisions at $\sqrt{s}=$ 13 TeV, JHEP 05 (2020) 033, [arXiv:1911.03947].
  * (55) ATLAS: G. Aad et. al., Search for new resonances in mass distributions of jet pairs using 139 fb-1 of $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector, JHEP 03 (2020) 145, [arXiv:1910.08447].
  * (56) ATLAS: M. Aaboud et. al., Search for low-mass dijet resonances using trigger-level jets with the ATLAS detector in $pp$ collisions at $\sqrt{s}=13$ TeV, Phys. Rev. Lett. 121 (2018) 081801, [arXiv:1804.03496].
  * (57) CDF: T. Aaltonen et. al., Search for new particles decaying into dijets in proton-antiproton collisions at $\sqrt{s}=1.96$ TeV, Phys. Rev. D 79 (2009) 112002, [arXiv:0812.4036].
  * (58) ATLAS: M. Aaboud et. al., Search for light resonances decaying to boosted quark pairs and produced in association with a photon or a jet in proton-proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector, Phys. Lett. B 788 (2019) 316–335, [arXiv:1801.08769].
  * (59) ATLAS Collaboration, Search for boosted resonances decaying to two b-quarks and produced in association with a jet at $\sqrt{s}=13$ TeV with the ATLAS detector, 2018.
  * (60) CMS: A. M. Sirunyan et. al., Search for low mass vector resonances decaying into quark-antiquark pairs in proton-proton collisions at $\sqrt{s}=$ 13 TeV, Phys. Rev. D 100 (2019) 112007, [arXiv:1909.04114].
  * (61) ATLAS: M. Aaboud et. al., Search for low-mass resonances decaying into two jets and produced in association with a photon using $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector, Phys. Lett. B 795 (2019) 56–75, [arXiv:1901.10917].
  * (62) CMS: A. M. Sirunyan et. al., Search for Low-Mass Quark-Antiquark Resonances Produced in Association with a Photon at $\sqrt{s}$ =13 TeV, Phys. Rev. Lett. 123 (2019) 231803, [arXiv:1905.10331].
  * (63) ATLAS: G. Aad et. al., Search for new phenomena in events with an energetic jet and missing transverse momentum in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector, arXiv:2102.10874.
  * (64) CMS collaboration, Search for new particles in events with energetic jets and large missing transverse momentum in proton-proton collisions at $\sqrt{s}=13~{}\mathrm{TeV}$, CMS-PAS-EXO-20-004 (2021).
  * (65) Fermi-LAT: M. Ackermann et. al., Searching for Dark Matter Annihilation from Milky Way Dwarf Spheroidal Galaxies with Six Years of Fermi Large Area Telescope Data, Phys. Rev. Lett. 115 (2015) 231301, [arXiv:1503.02641].
  * (66) Planck: N. Aghanim et. al., Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6, [arXiv:1807.06209].
  * (67) P. Gondolo and G. Gelmini, Cosmic abundances of stable particles: Improved analysis, Nucl. Phys. A 360 (1991) 145–179.
  * (68) T. Binder, T. Bringmann, M. Gustafsson, and A. Hryczuk, Early kinetic decoupling of dark matter: when the standard way of calculating the thermal relic density fails, Phys. Rev. D 96 (2017) 115010, [arXiv:1706.07433].
  * (69) D. E. Kaplan, M. A. Luty, and K. M. Zurek, Asymmetric Dark Matter, Phys. Rev. D 79 (2009) 115016, [arXiv:0901.4117].
  * (70) A. Pukhov, CalcHEP 2.3: MSSM, structure functions, event generation, batchs, and generation of matrix elements for other packages, hep-ph/0412191.
  * (71) A. Belyaev, N. D. Christensen, and A. Pukhov, CalcHEP 3.4 for collider physics within and beyond the Standard Model, Comp. Phys. Comm. 184 (2013) 1729–1769, [arXiv:1207.6082].
  * (72) T. Bringmann, J. Edsjö, P. Gondolo, P. Ullio, and L. Bergström, DarkSUSY 6 : An Advanced Tool to Compute Dark Matter Properties Numerically, JCAP 1807 (2018) 033, [arXiv:1802.03399].
  * (73) P. Gondolo, J. Edsjo, et. al., DarkSUSY: Computing supersymmetric dark matter properties numerically, JCAP 0407 (2004) 008, [astro-ph/0406204].
  * (74) A. Arbey and F. Mahmoudi, Dark matter and the early Universe: A review, Prog. Part. Nucl. Phys. 119 (2021) 103865, [arXiv:2104.11488].
  * (75) GAMBIT Dark Matter Workgroup: T. Bringmann, J. Conrad, et. al., DarkBit: A GAMBIT module for computing dark matter observables and likelihoods, Eur. Phys. J. C 77 (2017) 831, [arXiv:1705.07920].
  * (76) A. L. Fitzpatrick, W. Haxton, E. Katz, N. Lubbers, and Y. Xu, The Effective Field Theory of Dark Matter Direct Detection, JCAP 1302 (2013) 004, [arXiv:1203.3542].
  * (77) N. Anand, A. L. Fitzpatrick, and W. C. Haxton, Weakly interacting massive particle-nucleus elastic scattering response, Phys. Rev. C 89 (2014) 065501, [arXiv:1308.6288].
  * (78) J. B. Dent, L. M. Krauss, J. L. Newstead, and S. Sabharwal, General analysis of direct dark matter detection: From microphysics to observational signatures, Phys. Rev. D 92 (2015) 063515, [arXiv:1505.03117].
  * (79) GAMBIT Collaboration: P. Athron et. al., Global analyses of Higgs portal singlet dark matter models using GAMBIT, Eur. Phys. J. C 79 (2019) 38, [arXiv:1808.10465].
  * (80) S. Baum, R. Catena, J. Conrad, K. Freese, and M. B. Krauss, Determining dark matter properties with a XENONnT/LZ signal and LHC Run 3 monojet searches, Phys. Rev. D 97 (2018) 083002, [arXiv:1709.06051].
  * (81) T. Bringmann and C. Weniger, Gamma Ray Signals from Dark Matter: Concepts, Status and Prospects, Phys. Dark Univ. 1 (2012) 194–217, [arXiv:1208.5481].
  * (82) CTA: A. Acharyya et. al., Sensitivity of the Cherenkov Telescope Array to a dark matter signal from the Galactic centre, JCAP 01 (2021) 057, [arXiv:2007.16129].
  * (83) M. Bauer, M. Klassen, and V. Tenorth, Universal properties of pseudoscalar mediators in dark matter extensions of 2HDMs, JHEP 07 (2018) 107, [arXiv:1712.06597].
  * (84) N. Zhou, D. Berge, and D. Whiteson, Mono-everything: combined limits on dark matter production at colliders from multiple final states, Phys. Rev. D 87 (2013) 095013, [arXiv:1302.3619].
  * (85) A. J. Brennan, M. F. McDonald, J. Gramling, and T. D. Jacques, Collide and Conquer: Constraints on Simplified Dark Matter Models using Mono-X Collider Searches, JHEP 05 (2016) 112, [arXiv:1603.01366].
  * (86) GAMBIT Collider Workgroup: C. Balázs, A. Buckley, et. al., ColliderBit: a GAMBIT module for the calculation of high-energy collider observables and likelihoods, Eur. Phys. J. C 77 (2017) 795, [arXiv:1705.07919].
  * (87) M. R. Buckley, D. Feld, and D. Goncalves, Scalar Simplified Models for Dark Matter, Phys. Rev. D 91 (2015) 015017, [arXiv:1410.6497].
  * (88) J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer, and T. Stelzer, MadGraph 5 : Going Beyond, JHEP 06 (2011) 128, [arXiv:1106.0522].
  * (89) T. Sjostrand, S. Mrenna, and P. Z. Skands, A Brief Introduction to PYTHIA 8.1, Comput. Phys. Commun. 178 (2008) 852–867, [arXiv:0710.3820].
  * (90) C. Degrande, C. Duhr, et. al., UFO - The Universal FeynRules Output, Comput. Phys. Commun. 183 (2012) 1201–1214, [arXiv:1108.2040].
  * (91) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, FeynRules 2.0 - A complete toolbox for tree-level phenomenology, Comput. Phys. Commun. 185 (2014) 2250–2300, [arXiv:1310.1921].
  * (92) E. Conte, B. Fuks, and G. Serret, MadAnalysis 5, A User-Friendly Framework for Collider Phenomenology, Comput. Phys. Commun. 184 (2013) 222–256, [arXiv:1206.1599].
  * (93) CMS Collaboration, Simplified likelihood for the re-interpretation of public CMS results, CMS-NOTE-2017-001, 2017.
  * (94) ATLAS Collaboration, Reproducing searches for new physics with the ATLAS experiment through publication of full statistical likelihoods, ATL-PHYS-PUB-2019-029, 2019.
  * (95) GAMBIT Collaboration: P. Athron et. al., Combined collider constraints on neutralinos and charginos, Eur. Phys. J. C 79 (2019) 395, [arXiv:1809.02097].
  * (96) GAMBIT: P. Athron et. al., Thermal WIMPs and the scale of new physics: global fits of Dirac dark matter effective field theories, Eur. Phys. J. C 81 (2021) 992, [arXiv:2106.02056].
  * (97) M. J. Reid et. al., Trigonometric Parallaxes of High Mass Star Forming Regions: the Structure and Kinematics of the Milky Way, Astrophys. J. 783 (2014) 130, [arXiv:1401.5377].
  * (98) A. J. Deason, A. Fattahi, et. al., The local high-velocity tail and the galactic escape speed, MNRAS 485 (2019) 3514–3526, [arXiv:1901.02016].
  * (99) D. Djukanovic, K. Ottnad, J. Wilhelm, and H. Wittig, Strange electromagnetic form factors of the nucleon with $N_{f}=2+1$ $\mathcal{O}(a)$-improved Wilson fermions, Phys. Rev. Lett. 123 (2019) 212001, [arXiv:1903.12566].
  * (100) R. S. Sufian, Y.-B. Yang, et. al., Strange Quark Magnetic Moment of the Nucleon at the Physical Point, Phys. Rev. Lett. 118 (2017) 042001, [arXiv:1606.07075].
  * (101) GAMBIT Scanner Workgroup: G. D. Martinez, J. McKay, et. al., Comparison of statistical sampling methods with ScannerBit, the GAMBIT scanning module, Eur. Phys. J. C 77 (2017) 761, [arXiv:1705.07959].
  * (102) DARWIN: J. Aalbers et. al., DARWIN: towards the ultimate dark matter detector, JCAP 1611 (2016) 017, [arXiv:1606.07001].
  * (103) P. Scott, Pippi – painless parsing, post-processing and plotting of posterior and likelihood samples, Eur. Phys. J. Plus 127 (2012) 138, [arXiv:1206.2245].