Global fits of simplified models for dark matter with GAMBIT

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.


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 [1][2][3].The Weakly-Interacting Massive Particle (WIMP) hypothesis, in which DM is assumed to consist of a new species that interacts with a strength at least as weak as the weak nuclear force, is amongst the leading DM explanations due to its ability to explain the observed cosmological relic abundance of DM [4] at the same time as potentially being very tightly constrained in the near future by current experimental technologies [5].
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 EFT is valid up to some scale Λ, 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 [6][7][8][9][10][11] and direct [12][13][14][15][16] DM detection experiments, but difficult for higher-energy probes such as the ATLAS and CMS experiments at the Large Hadron Collider (LHC) [17][18][19][20][21][22][23].A detailed study of EFTs for Dirac DM has recently been performed by the GAMBIT collaboration [24].
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.[5,18,[25][26][27][28][29][30].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 [31][32][33] on the landscape of possible DM-SM interactions.
There are a plethora of studies of simplified DM models in the literature, including most notably a global scan of an s-channel vector-mediated Dirac fermion DM model performed by the MasterCode Collaboration [34].In this paper, we perform global fits of s-channel vector-mediated scalar, Dirac and Majorana fermion DM models, using GAMBIT v2.3 [35].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 [36].We complement the work of Ref. [34] 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 [24].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 DM 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 GAM-BIT input files and plotting scripts can be downloaded from Zenodo [37].

Models
We consider three models in this work, distinguished by the nature of the DM candidate: 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 [38].To enforce absolute stability of the DM candidate, we assume it to be odd under a new 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 DM masses.A detailed study of the case where this assumption breaks down for one model can be found in [39].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 couplings1 to avoid restrictive dilepton searches [42].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 q ) to the mediator.

Scalar DM
The Lagrangian density of our simplified model with a complex scalar DM candidate is where V µ is the mediator field, ϕ is the scalar DM candidate and F ′ µν is the mediator field strength tensor.The model has four free parameters: the DM mass m DM , the mediator mass m M , the mediator-WIMP coupling g V DM and the mediator-quark coupling g q .The DM candidate must be a complex scalar to avoid vanishing g V DM .This model is similar to one introduced in Ref. [43], differing by the absence of lepton couplings.
The decay width of the mediator to quark q, for all three models is while the decay width to the scalar DM candidate is

Dirac Fermion DM
The simplified model with Dirac fermion DM has the Lagrangian density where χ 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 V DM ) and axial-vector (g A DM ) couplings.We vary g V DM and g A 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 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. [42] for this: We apply this as an additional constraint, rejecting as theoretically invalid any parameter points that do not satisfy Eq. 6.

Majorana Fermion DM
The model Lagrangian of the BSM additions in the case of Majorana fermion DM is: 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 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 AT-LAS and CMS experiments.For this we employ GUM [36], 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 DM 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 L 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 DM, we simply give the best-case likelihood that results from predictions or parameters exactly matching their centrally measured values.
The details for these likelihoods, and their implementation in GAMBIT, are given in the following subsections.

Relic Abundance
The number density n DM,eq of DM particles in thermal equilibrium in the early Universe changes with time according to the Boltzmann equation [68] where section times the relative velocity and is given by where K 1,2 are the modified Bessel functions and v lab is the velocity of one of the annihilating DM particles or antiparticles in the rest frame of the other (see Ref. [69] for further discussion).Dirac fermion WIMPs give a total contribution to the observed DM density of n DM + n DM = 2n DM (where a possible initial asymmetry has been neglected [70]).We use GUM to generate tree-level annihilation crosssections for each model using CalcHEP v3.6.27 [71, 72].We obtain the WIMP relic density by using the Dark-Bit interface to DarkSUSY v6.2.5 [73, 74] to numerically solve Eq. ( 9) at each parameter point, assuming no departures from the standard cosmological history.For alternative scenarios, see Ref. [75].We compare the resulting value to the Planck 2018 measurement of Ω DM,obs h 2 = 0.120 ± 0.001 [67].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.[35,76].
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. [35].
We rescale all predicted direct detection signals by the fraction where Similarly, we rescale all indirect detection signals from DM-DM annihilation by f2 DM .Note that this encodes the two assumptions that f 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.

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 where ρ 0 is the local density of DM at our point in the Galactic halo, m 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 with µ = m T m DM /(m T + m DM ) 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 dσ/dE 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-relativistic EFT written as where the higher-dimensional operators O N i depend only on the spins of the WIMP and the nucleon ( ⃗ S DM and ⃗ S N ), the momentum transfer ⃗ q and the DM-nucleon relative velocity ⃗ v [12,77,78].A systematic classification of the non-relativistic operators that can result from the reduction of simplified models can be found in Refs [77,79], 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 N i (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 [76, 80] (through the DirectDM code, following Appendix A in [24]).We use DDCalc to compute the differential cross-section for each operator O N 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 p at each experiment, by convolving the differential rate with the product ϕ(E R ) of the energy resolution and detector acceptance, where M is the detector mass and T exp is the exposure time.

Indirect Detection
The processes that led to DM being in thermal equilibrium in the early Universe also allow for annihilation

Interaction
Effective Operator Relevant models Couplings to each operator follow [81].
in the present day.Annihilation products can thus potentially be seen originating from regions of high DM 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 DM.
An especially strong set of constraints on annihilating DM comes from observations of dwarf spheroidal galaxies [82].For a binned histogram of the DM-induced γ-ray flux, the flux from a target k in bin i can be written in the form Φ i • J k , where Φ i includes the relevant particle physics, and J k includes the relevant astrophysics (see Ref [76] for more details).
The velocity dependence of the annihilation crosssection 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 ̸ = 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 DM > m 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 DM > m M [5].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 mod-els that leads to significant indirect detection signals when m DM < m M .
The s-wave contribution gives a particle physics factor of where (σv) 0,j denotes the zero-velocity limit of the crosssection for DM, DM → j, N γ,j is the number of photons that results from the final state channel j (per annihilation), and f DM is the DM fraction defined in Eq. 11.The prefactor of 1/N DM reflects the nature of the DM candidate: N DM = 2 for self-conjugate particles and We use CalcHEP to compute the annihilation crosssections for each model, using the GUM interface.Each annihilation channel can produce either primary or secondary photons; the yields dN γ,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 Ω, 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 [66], computing the likelihood with the DarkBit interface to gamLike v1.0.1.The loglikelihood ln L exp is constructed from the product Φ i • J k summed over all targets and energy bins, An additional likelihood contribution comes from treating the J k factors of each dwarf spheroidal galaxy as nuisance parameters, giving ln The full likelihood, profiled over the J factors, is then given by ln L prof.dSphs = max An alternative to dwarf spheroidal measurements is to look for evidence of DM 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 [83].We briefly consider the future impact of CTA observations on the viable parameter space of our simplified models in Section 4.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 ( / E T ).In this study, we include CMS and ATLAS searches for monojet events based respectively on 137 fb −1 [65] and 139 fb −1 [64] 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 [84][85][86].
For a given parameter point of a simplified model, 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 where L is the integrated luminosity, σ is the total production cross-section, and ϵA is the product of the efficiency and acceptance for passing the kinematic selections that define the analysis.
The quantity ϵ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 [87].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 [88].We have therefore performed a more detailed simulation using MadGraph_aMC@NLO [89] (v3.1.1),interfaced to Pythia v8.3 [90] 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 [91], generated with FeynRules [92] and employing a 5-flavour scheme.We used MadAnalysis 5 [93] to perform detector simulation and implement each of the ATLAS and CMS monojet analyses in order to compute ϵA.As this set of simulations is too computationally expensive to run during the global fit, we precomputed grids of the cross sections (σ) and ϵA factors for each LHC experiment in advance, and interpolated them at runtime using Col-liderBit 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: - 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 [94], which approximates 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: 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 γ i quantifies the deviation from the nominal expected yields due to systematic uncertainties in signal region i.The set {γ i } thus gives 66 nuisance parameters in total for this particular analysis.The covariance matrix Σ provided by CMS encodes the correlations between the various γ 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 γ i parameters so that the final CMS likelihood is defined solely in terms of the simplified model signal estimates s: where γ denotes the combination of background nuisance parameters resulting in the highest value of the likelihood for the given signal s.
The ATLAS analysis does not come with a published covariance matrix, nor with a published likelihood in the HistFactory format of Ref. [95].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 / E T > 1000 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 where L ATLAS (s i , γ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 L LHC = ln L CMS + ln L 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 Col-liderBit solution is to instead set the LHC log-likelihood contribution to the total difference in log-likelihood between the signal and background-only (s = 0) cases: 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σ and 2σ confidence regions defined using the likelihood ratio L(Θ)/L(Θ best-fit ) around the best-fit parameters Θ 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 More detail on this procedure can be found in Ref. [96].
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, q , 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. [34].We interpolate these published limits in m 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 q,excl ) provided by a broad range of LHC dijet searches [55-63], using a likelihood of the form ln The factor of −2 arises because g q,excl is taken from the 95% confidence limit on the coupling provided by each analysis, corresponding to ∆χ 2 = 4.This produces coupling upper limits as shown in Figure 1.
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.

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.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).
Our treatment of the local DM density ρ 0 follows the default prescription in DarkBit, which assumes that ρ 0 is log-normally distributed with central value ρ 0 = 0.40 GeV cm −3 and error σ ρ0 = 0.15 GeV cm −3 .The asymmetric scan range for ρ 0 in Table 3 reflects this log-normal distribution.All other nuisance parameters are scanned over their 3σ 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 DM EFTs [80,97].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 peak = 240 ± 8 km s −1 [98] and v esc = 528 ± 25 km s −1 (based on Gaia data [99]), respectively.

Results
We have performed comprehensive scans of each simplified model parameter space using the differential evolution sampler Diver v1.0.4 [100] 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.

Scalar DM
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σ 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 DM > m M , where t-channel annihilation into a pair of mediators is possible.In regions either side of the diagonal resonance region, DM is overproduced and the model is inconsistent with relic density measurements.
Highest likelihood points lie along the resonance, where 2m DM ≈ m 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 lies (see Table 4).
When requiring that the scalar DM candidate explains all of DM (Figure 2), mediator masses up to approximately 2 TeV are excluded.Toward lower mediator masses the strength of the effective coupling used to calculate direct detection signals is increased.To escape direct detection bounds, the highest likelihood moves toward regions that underproduce the relic abundance.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 model was allowed to underproduce the DM abundance.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 q in the scan, and reducing this lower bound will expand the allowed region.
Direct detection limits give the lower bound on DM masses, along the resonance region in both scans and also when m DM > m M in Figure 2 (top) since these experiments are more constraining for lighter DM masses (provided that the mass is well above the energy threshold of the experiment).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.
Indirect detection limits give a slight preference to the region along the resonance.This is because when m DM > m 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 M .
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.

Dirac Fermion DM
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 DM ≈ m M .
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 bestfit point that would contribute to the CMS monojet search [65].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  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 A 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 DM > m M .The model overproduces DM below the resonance region, where m M > 2m DM .Requiring that the DM candidate constitute all of the observed DM excludes regions along the resonance where m DM ≲ 300 GeV and regions off the resonance where m M ≲ 600 GeV, where the model cannot escape direct detection bounds without underproducing the DM abundance.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σ 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 offresonance region seen in the capped collider results (Figure 4 right) and the uncapped collider results when requiring a saturated relic abundance (Figure 4 bottom left).To avoid unitarity violation, the model is excluded for larger g A 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 V DM .Since the strength of the direct detection signals is primarily from the g V DM coupling, and the unitarity bound is from the g A DM coupling, this shape in parameter space would differ in either the pure vector or pure axial-vector coupling cases (as studied in Ref. [34]).If the g 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. [34] 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 [34] without presenting duplication of their results.

Majorana Fermion DM
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σ contours, with little exclusion when m DM > m 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 V 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 such that abundance measurements exclude regions with a high mediator mass and low DM mass. 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 DM > m M region in the fits with capped L 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 DM < m M .

Future prospects
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 [83].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 DM > m M .
We show the predicted number of signal events in DARWIN, a next-generation direct detection experi-ment [101], as a function of the DM mass in Figure 11.Each model has the potential to produce tens or even more than a hundred events.DARWIN will therefore prove useful for constraining 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.

Conclusions
In this work we performed global scans with GAMBIT of three simplified DM 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 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

Fig. 1 :
Fig.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).

Fig. 2 :
Fig. 2: Scalar DM profile likelihood, profiling over couplings.The observed relic density of DM is taken as an upper limit (top) or to consist entirely of the scalar DM candidate (bottom).1σ and 2σ contours are shown in white and grey respectively.

Fig. 3 :
Fig. 3: Relic Density of the scalar DM model as a function of the mediator mass (left) and DM mass (right).The results are shown for the scan where the observed relic density of DM is allowed to exceed the abundance of the scalar DM candidate.1σ and 2σ contours are shown in white and grey respectively.

Fig. 4 :
Fig. 4: Profile likelihood the Dirac fermion DM model.The observed relic density of DM 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σ and 2σ contours are shown in white and grey respectively.
capped (Figure4, right), then the surviving parameter space opens up.

Fig. 5 :
Fig. 5: Missing energy spectra for the CMS monojet search.The SM background prediction (purple) and the observed event counts (black) are taken from Ref [65].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.

Fig. 6 : 4 Fig. 7 :
Fig. 6: Profile likelihood as a function of relic density of the Dirac fermion DM model, allowing the model to underproduce the relic abundance.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σ and 2σ contours are shown in white and grey respectively.

Fig. 8 :
Fig. 8: Profile likelihood for the Majorana fermion DM model.The observed relic density of DM 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σ and 2σ contours are shown in white and grey respectively.

Experiment ln L bg ln L max
(t)is the Hubble rate and n DM n DM ≡ n 2 DM for Majorana DM. ⟨σv rel ⟩ is the thermally averaged cross-

Table 2 :
Effective Operator matching to each model.The axial-vector couplings are momentum or velocity suppressed.

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 µ = 2 GeV.

Table 4 :
Approximate best-fit points for each model.∆lnL values are defined as ln L − ln L bg , where the background-only likelihood is detailed in Table1.