Global Analysis of Dark Matter Simplified Models with Leptophobic Spin-One Mediators using MasterCode

TL;DR

This paper conducts a comprehensive global analysis of leptophobic spin-one dark matter models, integrating cosmological, direct detection, and collider constraints to identify viable parameter regions and predict experimental signatures.

Contribution

It provides the first full parameter space scan of leptophobic spin-one mediators, combining multiple experimental constraints within the MasterCode framework.

Findings

Two main allowed regions in parameter space identified.

Mediator and dark matter masses could be within LHC reach.

Predictions made for spin-independent and -dependent scattering cross sections.

Abstract

We report the results of a global analysis of dark matter simplified models (DMSMs) with leptophobic mediator particles of spin one, considering the cases of both vector and axial-vector interactions with dark matter (DM) particles and quarks. We require the DMSMs to provide all the cosmological DM density indicated by Planck and other observations, and we impose the upper limits on spin-independent and -dependent scattering from direct DM search experiments. We also impose all relevant LHC constraints from searches for monojet events and measurements of the dijet mass spectrum. We model the likelihood functions for all the constraints and combine them within the MasterCode framework, and probe the full DMSM parameter spaces by scanning over the mediator and DM masses and couplings, not fixing any of the model parameters. We find, in general, two allowed regions of the parameter spaces: one in which the mediator couplings to Standard Model (SM) and DM particles may be comparable to those in the SM and the cosmological DM density is reached via resonant annihilation, and one in which the mediator couplings to quarks are $\lesssim 10^{-3}$ and DM annihilation is non-resonant. We find that the DM and mediator masses may well lie within the ranges accessible to LHC experiments. We also present predictions for spin-independent and -dependent DM scattering, and present specific results for ranges of the DM couplings that may be favoured in ultraviolet completions of the DMSMs.

Figures

Tables

Table 1. Table 1: The ranges of the DMSM parameters sampled, together with the numbers of segments into which they were divided during the sampling.

Parameter	Range	# of Segments
$m_Y$	$(0.1,5)\mathrm{TeV}$	10
$m_{\chi }$	$(0,2.5)\mathrm{TeV}$	8
$g_{\mathrm{SM}}$	$(0,\sqrt{4\pi })$	2
$g_{\mathrm{DM}}$	$(0,\sqrt{4\pi })$	2
Total # of segments		320

Equations

$${\mathcal{L}}_{\chi }^1=\bar{\chi }{\gamma }_{\mu }\left(g_{\chi }^V+g_{\chi }^A{\gamma }_5\right)\chi Y^{\mu },$$

$${\mathcal{L}}_q^1=\sum _i\overline{q_i}{\gamma }_{\mu }\left(g_{q_i}^V+g_{q_i}^A{\gamma }_5\right)q_i{Y}^{\mu }.$$

$$g_{\chi }^V\equiv g_{\mathrm{DM}}\ne 0,$$

$$g_q^V\equiv g_{\mathrm{SM}}\ne 0,$$

$$g_{\chi }^A\equiv g_{\mathrm{DM}}\ne 0,$$

$$g_q^A\equiv g_{\mathrm{SM}}\ne 0,$$

$${\Omega }_{\mathrm{CDM}}{h}^2=0.120\pm 0.001.$$

$$⟨\sigma v{⟩}_s\simeq ⟨\frac{c_Y}{32\pi }\frac{g_{\mathrm{SM}}^2{g}_{\mathrm{DM}}^2{m}_{\chi }^2}{(m_{\chi }^2{v}^2+4m_{\chi }^2-m_Y^2{)}^2+m_Y^2{\Gamma }_Y^2}⟩,$$

$${\Gamma }_Y^{XX}=\frac{g_X^2{m}_Y}{12\pi }\cdot f\left(\frac{m_X^2}{m_Y^2}\right),$$

$$f(x)=\{\begin{array}{ll}(1+2x)(1-4x{)}^{1/2}& (\mathrm{vector})\\ (1-4x{)}^{3/2}& (\mathrm{axial}\mathrm{vector})\end{array},$$

$${\Omega }_{\chi }{h}^2\text{mbox}\underset{\text{vbox}\sim }{\text{vbox}\propto }\frac{1}{⟨\sigma v{⟩}_s}\text{mbox}\underset{\text{vbox}\sim }{\text{vbox}\propto }\frac{m_Y^2}{g_{\mathrm{SM}}^2{g}_{\mathrm{DM}}^2},$$

$$1/[(m_Y^2-4m_{\chi }^2{)}^2+m_Y^2{\Gamma }_Y^2]$$

$$\to \delta (m_Y^2-4m_{\chi }^2)/(m_Y{\Gamma }_Y),$$

$${\Omega }_{\chi }{h}^2\text{mbox}\underset{\text{vbox}\sim }{\text{vbox}\propto }\frac{1}{⟨\sigma v{⟩}_s}\text{mbox}\underset{\text{vbox}\sim }{\text{vbox}\propto }\frac{m_Y^2}{g_{\mathrm{DM}}^2}$$

$${\Omega }_{\chi }{h}^2\text{mbox}\underset{\text{vbox}\sim }{\text{vbox}\propto }\frac{m_Y^2{g}_{\mathrm{SM}}^2}{g_{\mathrm{DM}}^2}$$

$${\Omega }_{\chi }{h}^2\text{mbox}\underset{\text{vbox}\sim }{\text{vbox}\propto }\frac{m_Y^2}{g_{\mathrm{SM}}^2}$$

$$⟨\sigma v{⟩}_t\simeq \frac{g_{\mathrm{DM}}^4}{32\pi }\frac{(m_{\chi }^2-m_Y^2{)}^{3/2}}{m_{\chi }(2m_{\chi }^2-m_Y^2{)}^2},$$

$${\sigma }_p^{\mathrm{SI}}\simeq {\left(\frac{g_{\mathrm{DM}}}{0.1}\right)}^2{\left(\frac{g_{\mathrm{SM}}}{0.1}\right)}^2{\left(\frac{m_Y}{1\mathrm{TeV}}\right)}^{-4}10^{-43}[{\mathrm{cm}}^{-2}]$$

$$\Delta {\chi }^2=5.99\times {\left(\frac{1}{R_i^{\mathrm{UL}}(\mathbf{m})}\frac{{\sigma }_{\mathrm{MG}}(\mathbf{m})}{{\sigma }_{\mathrm{MG}}^{\mathrm{fix}}(\mathbf{m})}\right)}^2,$$

$${\sigma }_{Y\to qq}\le \sigma (g_q^{\ast })$$

$${\sigma }_{Y\to qq}=c\frac{g_{\mathrm{SM}}^4}{{\Gamma }_q+{\Gamma }_{\chi }},\sigma (g_q^{\ast })=c\frac{(g_q^{\ast }{)}^4}{{\Gamma }_q^{\ast }}$$

$$g_{\mathrm{SM}}^{\ast }\equiv \sqrt{\frac{(g_q^{\ast }{)}^2}{2}\left(1+\sqrt{1+\frac{4{\Gamma }_{\chi }}{{\Gamma }_q^{\ast }}}\right)},$$

$$\Delta {\chi }^2=4\times {\left(\left[\frac{g_{\mathrm{SM}}^4}{{\Gamma }_q+{\Gamma }_{\chi }}\right]/\left[\frac{(g_{\mathrm{SM}}^{\ast }{)}^4}{{\Gamma }_q(g_{\mathrm{SM}}^{\ast })+{\Gamma }_{\chi }}\right]\right)}^2$$

$${\Omega }_{\chi }{h}^2\sim 0.1{\left(\frac{m_Y}{2\mathrm{TeV}}\right)}^2{\left(\frac{g_{\mathrm{DM}}}{0.01}\right)}^{-2}$$

$$g_{\mathrm{SM}}\lesssim 0.1\cdot \left(\frac{m_{\chi }}{1\mathrm{TeV}}\right){\left[\frac{{\sigma }_{\mathrm{SI}}^{\mathrm{UL}}(m_{\chi })}{10^{-45}[{\mathrm{cm}}^{-2}]}\right]}^{\frac{1}{2}}$$

$$\Big{(}\frac{m_{Y}}{2\,\rm TeV}\Big{)}^{2}\Big{(}\frac{g_{\rm DM}}{0.01}\Big{)}^{-2}\lesssim 1\,,$$

Full text

KCL-PH-TH/2019-10, CERN-TH-2019-007, DESY-19-071, PSI-PR-19-06, IFT-UAM/CSIC-18-120

FTPI-MINN-19/05, UMN-TH-3814/19

1 Introduction

The nature of Dark Matter (DM) is one of the most pressing issues in contemporary physics. For over 80 years, astrophysical and cosmological observations have indicated its existence indirectly [1, 2]. They have also suggested that DM is mainly cold, i.e., it has been non-relativistic since at least the epoch of cosmological recombination [3]. Beyond this, we know very little, and candidates for DM range in mass from black holes heavier than the Sun, to various species of (meta)stable neutral particles, to very light particles with large occupation numbers best described by classical fields, see Ref. [4] for a review. If the DM is composed of weakly-interacting massive particles (WIMPs) that were in thermal equilibrium with Standard Model (SM) particles in the early Universe, freeze-out calculations suggest that the WIMP is likely to weigh ${\cal O}$(TeV), in which case it could be produced at accelerators, notably the LHC. Thus, the search for WIMP DM particles has been one of the principal research objectives of the LHC experiments, particularly ATLAS and CMS [5, 6, 7, 8, 9].

Several approaches have been taken to DM searches at the LHC, in both predictions for possible signals and interpretations of the search limits. Initially, many analyses were based on specific models that predict WIMPs capable of providing the cold DM, such as supersymmetry (SUSY) [10] and some other models with new physics at the TeV scale [11]. These searches were typically either for the production of heavier new particles followed by cascade decays to the dark matter particle, or direct production of dark matter particles in association with a single SM particle used for tagging purposes (the ’mono-X’ signature). DM models giving rise to long-lived signatures have also been studied recently [12].

We have used the MasterCode framework [13] previously to study several SUSY scenarios with various experimental signatures [14], and other groups have made parallel analyses [15]. In particular, in Ref. [14] we provided detailed predictions for the preferred parameter spaces in various concrete SUSY models. However, there are many other scenarios for possible physics beyond the SM, and it is desirable to fashion search strategies that facilitate the interpretation of results with differing characteristics, minimizing theoretical biases.

Recently, growing attention is being paid to more general model frameworks and analyses motivated by more generic experimental signatures. One of the first such approaches was to construct an effective field theory (EFT) for the interactions between DM and SM particles [16, 17, 18, 19], with a focus on the mono-X signatures. A shortcoming of this EFT approach is that, in order to obtain the appropriate cosmological DM density, the DM/SM interactions are likely to be mediated by particles in the TeV mass range, for which the approximation of a contact interaction may be inadequate. In particular, the LHC might also be capable of producing the mediator particle directly, providing an additional signature beyond the EFT framework [20, 21]. For this reason, interest has been growing in Dark Matter Simplified Models (DMSMs) [21, 22, 23, 24, 25, 26, 27, 28, 29], which postulate effective Lagrangians that include explicitly the mediator particle and its interactions with both DM and SM particles. Recommendations for the formulation of DMSMs for use by the LHC experiments have been made by the LHC Dark Matter Working Group (DMWG) [5, 6, 7, 8, 9].

In this paper we report the results of global likelihood analyses within the MasterCode framework of some of the DMSMs suggested by the DMWG and analyzed by ATLAS and CMS.

We take a bottom-up approach, postulating a mediator particle in the LHC mass range, without requiring an explicit UV completion, commenting later how the models investigated here could arise in UV-complete scenarios. We consider only $s$-channel mediators of spin one. Avoiding the strong LHC constraints on mediators interacting with charged leptons [30, 31], we assume that their interactions are leptophobic. We assume that the DM particle is a neutral Dirac fermion $\chi$. The Lagrangian for the spin-one mediator $Y$ interactions with the DM particle $\chi$ is

[TABLE]

and that for its interactions with quarks is

[TABLE]

For simplicity, following the DMWG [8] we assume that the interactions between the $Y$-boson and quarks are generation-independent, and we also assume that they are the same for charge-2/3 and -1/3 quarks ($g_{q}^{V,A}:=g_{q_{i}}^{V,A}\;\forall\;i$) 111We do not discuss here the cancellation of gauge anomalies, which is possible with a suitable non-minimal ‘dark sector’, whose experimental signatures are more model-dependent and not discussed in detail here [32, 33, 34].. Being leptophobic, the $Y$ boson has no interactions with leptons. We consider two scenarios for the $Y$ couplings, one in which they are purely vectorial:

[TABLE]

and another in which they are purely axial:

[TABLE]

Note that the axial vector model could in principle violate unitarity at high energies. However, the couplings required for this to occur are larger than the range we consider in our scans, and would even call the particle interpretation of the mediator into question [26, 35].

UV completions of the DMSMs we study could be achieved in a number of different ways. These could favour specific ranges of $g_{\rm SM}$, $g_{\rm DM}$ and their ratio, as we explore below. Since our spin-1 models involve massive vector bosons, they will likely involve a dark Higgs responsible for giving them masses, or else a Stückelberg mechanism [36]. However, our intention here is not to fit the parameters of a complete model, but rather to study correlations between the model parameters defined above, which are the minimal sets leading to an interplay between LHC and DM phenomenology that is of interest.

Each of these DMSMs has four free parameters, the masses of the DM particle and the mediator, $m_{\chi}$ and $m_{Y}$, and the mediator couplings to the DM and SM particles, $g_{\rm DM}$ and $g_{\rm SM}$. Sampling these 4-dim parameter spaces is computationally tractable. However, experimental results are often interpreted for fixed values of the couplings, and thus previous investigations have generally fixed two out of the four parameters. A more general approach is desirable for combining the direct DM constraints with those from the LHC, particularly because obtaining the preferred cosmological value of $\Omega_{\chi}h^{2}$ requires values of $g_{\rm DM}$ and $g_{\rm SM}$ that depend on $m_{\chi}$ and $m_{Y}$. Consequently, without fixing any parameter of the model, we implement the predictions of the models using micrOMEGAs [37] and Madgraph5_aMC@NLO [38] using DMSIMP model files [39] to build a likelihood function for each of the astrophysical and accelerator constraints, and combine them using the MasterCode framework. We use MultiNest[40] to identify the regions where the global likelihood is minimized. We do not include in our analysis indirect searches for astrophysical DM particles via their annihilations into SM particles, which were found in [34] to be unimportant in simple DM models when $m_{\chi}\gtrsim 50\,\,\mathrm{GeV}$.

The layout of this paper is as follows. In Section 2 we set out the constraints we use, and in Section 3 we describe our analysis framework. Our results are described in Section 4, Possible UV completions are discussed in Section 5, and our conclusions and some perspectives are set out in Section 6.

2 Astrophysical and LHC Constraints

2.1 Dark Matter Constraints

2.1.1 Dark Matter density

The mean density of cold DM in the Universe is tightly constrained by Planck measurements of the cosmic microwave background and other observations [41]:

[TABLE]

We assume here that the dominant source of this cold DM is the WIMP $\chi$, so that $\Omega_{\chi}h^{2}\simeq\Omega_{\rm CDM}h^{2}$, with interactions that are described in the Lagrangians (3) and (4) in the DMSMs we study.

The relic abundance is approximately proportional to the inverse of the annihilation cross-section, $\Omega_{\chi}h^{2}\propto 1/\langle\sigma v\rangle$. If the annihilation is dominated by the $s$-channel process $\chi\chi\to Y^{(*)}\to$ Standard Model particles, the leading expression for the cross-section takes the form [42, 43]

[TABLE]

where $c_{Y}$ is an ${\cal O}(1)$ constant for the vector mediator, whereas for the axial mediator it is suppressed by the quark mass and the relative velocity, $v$, of the annihilating DM particles: $c_{Y}\sim a\frac{m_{q}^{2}}{m_{\chi}^{2}}+bv^{2}$. The total width of $Y$ is given by $\Gamma_{Y}=\Gamma_{Y}^{\chi\chi}+N_{c}\sum_{q}\Gamma_{Y}^{qq}$, where $N_{c}=3$ is a colour factor and

[TABLE]

where $X=\chi\,({\rm DM}),q\,({\rm SM})$, and

[TABLE]

where $x\equiv(m_{X}/m_{Y})^{2}$. The width of $Y$ becomes particularly important in the resonant region, where $m_{Y}\simeq 2m_{\chi}$ and $\Gamma_{Y}^{\chi\chi}\ll\Gamma_{Y}$ because of phase-space suppression.

In the axial-vector case, $\langle\sigma v\rangle_{s}$ is approximately proportional to $v^{2}$ through $c_{Y}$ and, if $g_{\rm SM}\ll 1$, the velocity term dominates the denominator of Eq. (6) at the resonance. This leads to the approximate scaling

[TABLE]

where \mathrel{\raisebox{3.01385pt}{\mbox{\underset{\vbox{\hbox{$\sim$}}}{\vbox{\hbox{$\propto$}}}}}} denotes approximate proportionality. On the other hand, in the vector case the annihilation is $s$-wave and

the narrow-width approximation is applicable for $g_{\rm SM}\ll 1$, and one finds that in the limit $g_{\rm SM}\to 0$

[TABLE]

so the scaling of the relic density is expected to be

[TABLE]

near the $Y$ peak. However, we recall that $\chi\chi$ annihilation takes place in a thermal bath, with collisions taking place at a range of centre-of-mass energies $\geq 2m_{\chi}$. For this reason, annihilations at an effective centre-of-mass energy within $\mathcal{O}(\Gamma_{Y})$ of $m_{Y}$ do not dominate when $\Gamma_{Y}$ is very small. We find numerically that (11) captures the $g_{\rm DM}$ dependence of $\Omega_{\chi}h^{2}$, and that it depends only weakly on $g_{\rm SM}$ over an intermediate range of $g_{\rm SM}$, but that $\Omega_{\chi}h^{2}$ increases at both large and small $g_{\rm SM}$ :

[TABLE]

for $g_{\rm SM}\gg g_{\rm DM}$ and

[TABLE]

for $g_{\rm SM}\ll g_{\rm DM}$. The dependence on $g_{\rm DM}$ in this case depends on the relative sizes of the width and mass difference between $m_{Y}$ and $2m_{\chi}$ [43]. If $\Gamma_{Y}/m_{Y}>(1-4m_{\chi}^{2}/m_{Y}^{2})$ and $m_{Y}>2m_{\chi}$ then $\Omega_{\chi}h^{2}$ is proportional to $g_{\rm DM}^{2}$. However, if the width is very small (for example if $g_{\rm SM}$ is very small and $m_{Y}<2m_{\chi}$), then the denominator in (6) is dominated by the mass difference and $\Omega_{\chi}h^{2}$ is inversely proportional to $g_{\rm DM}^{2}$, as seen in Eq. (9) for the axial-vector case.

Another important process is the $t$-channel process $\chi\chi\to YY$ for $m_{\chi}>m_{Y}$. The leading contribution to the $t$-channel cross-section is given by [42]

[TABLE]

and is kinematically forbidden when $m_{\chi}<m_{Y}$.

We use micrOMEGAs [37] to evaluate $\Omega_{\chi}h^{2}$ numerically in each model parameter set, including a ${\pm 10}$% range around Eq. (5). This enables us to obtain a substantial but unbiased sample of model points.

As a first validation step, we show in Fig. 1 a comparison between our micrOMEGAs calculations of $\Omega_{\chi}h^{2}$ (colours and solid blue line) with the $\Omega_{\chi}h^{2}=0.12$ curve shown by the CMS Collaboration in [44] (dashed line) for the case of a spin-one mediator with vector-like couplings $g_{\rm SM}^{V}=0.25,g_{\rm DM}^{V}=1$. As seen in Fig. 1, the agreement is excellent, providing a valuable cross-check on our calculations. We note that the red shaded region at larger values of $m_{Y}$ is excluded because $\Omega_{\chi}h^{2}>0.12$, and that the DM particle is underdense in the blue region at large $m_{\chi}$.

2.1.2 Spin-independent DM scattering

We evaluate the constraint on the spin-independent DM scattering cross-section $\sigma^{\rm SI}_{p}$ from a combination of the LUX [45], PandaX-II [46] and XENON1T [47] experiments using a joint likelihood function, evaluating the theoretical prediction for $\sigma^{\rm SI}_{p}$ using micrOMEGAs. This constraint is important for the vector mediator case, where the spin-independent cross-section has the following scaling [7] $\sigma^{\rm SI}_{p}\propto g_{\rm DM}^{2}g_{\rm SM}^{2}/m_{Y}^{4}$, with the following approximate numerical value:

[TABLE]

in the $s$-channel region where $m_{Y}\simeq 2m_{\chi}$. The $\sigma^{\rm SI}_{p}$ constraint provides important limits on the vector mediator mass $m_{Y}$ and the product $g_{\rm DM}g_{\rm SM}$.

We allow for an overall uncertainty of 10% in the spin-independent cross section.

2.1.3 Spin-dependent DM scattering

We also evaluate the constraint imposed by the upper limit on the spin-dependent DM-proton scattering cross-section $\sigma^{\rm SD}_{p}$ from the PICO-60 experiment [48], and that on the spin-dependent DM-neutron scattering cross-section $\sigma^{\rm SD}_{n}$ from the XENON1T experiment [49]. We evaluate theoretical predictions for $\sigma^{\rm SD}_{p}$ and $\sigma^{\rm SD}_{n}$ using micrOMEGAs, which uses estimates of the spin-dependent scattering matrix elements consistent with [50].

We also allow for an overall uncertainty of 10% in the spin-dependent cross section.

2.1.4 Halo model

The majority of direct searches for DM scattering assume a standard halo model (SHM) in which the local density of cold DM is 0.3 GeV/cm3, with negligible uncertainty. However, several recent analyses favour a larger local density, and here we follow [51] in assuming a central value for the local DM density of 0.55 GeV/cm3. For this reason, the direct DM scattering limits, estimated future sensitivity curves and neutrino ‘floors’ we use are rescaled relative to the published curves.

It has been suggested in a recent analysis of Gaia data that there may be a local debris flow that modifies the local velocity distribution from that given by the SHM [52]. However, this modification of the the local velocity distribution was shown to have only a small effect on the interpretation of the XENON1T experiment for $m_{\chi}$ in the range of interest in this paper, and a similar conclusion was reached in [51]. For another estimate of the local density of DM based on Gaia data, see [53].

When implementing the experimental constraints on $\sigma^{\rm SI}_{p}$ and $\sigma^{\rm SD}_{p},\sigma^{\rm SD}_{n}$, we assume an overall uncertainty of 30% [54, 51], in the local density.

2.2 LHC Experimental Constraints

2.2.1 Monojet Constraints

We use the constraints provided by the CMS Collaboration in [44] using 35.9/fb of data from collisions at 13 TeV, based on an analysis employing signal regions targeting monojet final states. In particular, we use the 95 % CL upper bounds $R_{i}^{\rm UL}({\bf m})$, where $i$ labels the vector and axial-vector mediators, imposed on $\sigma({\bf m})/\sigma^{\rm fix}({\bf m})$ evaluated at each point in the ${\bf m}=(m_{\chi},m_{Y})$ plane, where $\sigma^{\rm fix}$ is the cross-section evaluated with $(g_{\rm SM},g_{\rm DM})=(0.25,1)$.

We model the likelihood function from this search using the procedure outlined in Fastlim [55]:

[TABLE]

where $\sigma_{\rm MG}$ is the monojet cross-section calculated using MG5$\underline{~{}}$aMC@(N)LO [38].

As a second step in our validation, we display in Fig. 2 our implementations (colours and solid red lines) of the CMS monojet constraints on the spin-1 mediator with vector-like interactions (left panel) and axial couplings (right panel). We find excellent agreement with the CMS limits [44], which are indicated by dashed red lines.

2.2.2 Dijet Constraints

The ATLAS and CMS Collaborations have both published constraints from dijet invariant-mass distributions on $Z^{\prime}$ resonances [56, 57, 58, 59, 60, 61] with up to 139/fb of data at 13 TeV [62]. We use a combination of the limits presented in the plane of the $Z^{\prime}$ mass and its ${\bar{q}}q$ coupling, $g_{q}$. In recasting this constraint for the $pp\to Y\to qq$ process in our DM models, we demand

[TABLE]

at $m_{Y}=m_{Z^{\prime}}$, where $g_{q}^{*}$ is the upper limit placed on the coupling of the $Z^{\prime}$ to quarks at a given $m_{Z^{\prime}}$. We note that $\sigma_{Y\to qq}$ and $\sigma({g_{q}^{*}})$ can be written as (using $\Gamma_{x}\equiv\Gamma(Y\to xx)$)

[TABLE]

with the same constant $c$ at $m_{Y}=m_{Z^{\prime}}$, and \Gamma_{q}=\big{(}\frac{g_{\rm SM}}{g_{q}^{*}}\big{)}^{2}\Gamma_{q}^{*}. Using these relations, the equality in Eq. (17) occurs when $g_{\rm SM}$ takes the value

[TABLE]

which serves as the upper bound on $g_{\rm SM}$ for given $m_{Y}$, $\Gamma_{q}$ and $\Gamma_{\chi}$. We model the corresponding likelihood function using

[TABLE]

where $\Gamma_{q}$ and $\Gamma_{\chi}$ are obtained with MG5$\underline{~{}}$aMC@(N)LO.

As a further step in our validation, we have cross-checked our recasting with the CMS and ATLAS limits [56, 57, 58, 59, 60, 61, 62], as shown in Fig. 3 222The discontinuity in the limit at $m_{Y}=2\,\,\mathrm{TeV}$ is due to the improved sensitivity of the recent ATLAS limit using 139/fb of data at 13 TeV [62].. The coloured lines correspond to the various experimental analyses, which can directly compared to the MasterCode $\chi^{2}$ evaluation indicated by colours in the plane. Note that we use the most sensitive upper limit for each mass value, and do not attempt to combine different experiments. Since the width of the band in Fig. 3 where the variation in the likelihood function is significant (as indicated by the colour transition) is typically smaller than the differences between the various experimental limits, this should be a good approximation. The excellent agreement between our 95% CL upper limit on $g_{\rm SM}$ and that quoted by the experiments is a valuable cross-check of our implementation of the dijet invariant-mass constraints.

In addition to these invariant-mass searches, the CMS Collaboration has also published a limit on $g_{\rm SM}$ as a function of $m_{Z^{\prime}}$ from an analysis of dijet angular distributions [63]. However, this limit is weaker than the limits from the dijet invariant mass distribution for the range of masses we study. Accordingly, we do not include it in our analysis.

3 Analysis Framework

Global-fitting frameworks such as MasterCode [14] have been used to study the constraints on different possible extensions of the SM, including SUSY models in particular, but can also be applied to other models with WIMP candidates for DM. Since DMSMs are not complete models, we do not consider it meaningful to compute $p$-values. Instead, here we use the MasterCode framework to correlate the impacts of constraints from different sectors, explore correlations and find regions of the model parameters that are still allowed.

We recall that MasterCode is a frequentist framework for combining constraints, written in Python/Cython and C++. As already mentioned, we use MG5$\underline{~{}}$aMC@(N)LO [38] to calculate mediator properties and the collider constraints on models implemented using DMSIMP [39], and we use micrOMEGAs [37] for DM density and scattering calculations. We use the MultiNest [40] algorithm for efficient sampling of the model parameter spaces. Since each of the DMSMs (3) and (4) that we study has a parameter space of only 4 dimensions, and since the constraint set is not large, sampling the model parameter spaces is not computationally onerous. In studies for our analysis we have made use of udocker software framework [64] to automatize the deployment of MasterCode inside Linux containers. This is a middleware suite developed in the context of the INDIGO datacloud project [65] to run docker containers in userspace, without requiring root privileges for installation or for execution.

The ranges of DMSM parameters that we study are shown in Table 1, together with the numbers of segments we use for our basic sampling of the parameter space. The range of $m_{Y}$ was chosen to avoid the low-mass region where mixing with the $Z$ could be subject to important constraints from precision electroweak data 333See [34] for a treatment of these constraints, which we discuss in more detail below. and indirect searches for astrophysical DM annihilations should be considered 444Searches for $\gamma$-rays from hadronic DM annihilations are generally insensitive to the cross-section required to obtain the correct cosmological DM density for $m_{\chi}\gtrsim 50\,\,\mathrm{GeV}$ [34], see also [66] and references therein. Indirect constraints from searches for energetic solar neutrinos are not competitive with direct searches for spin-dependent DM scattering, as discussed below., but include all the masses for which LHC searches are sensitive. The couplings $g_{\rm SM}$ and $g_{\rm DM}$ were restricted to perturbative ranges $<\sqrt{4\pi}$.

In addition to the generic sampling ranges shown in Table 1, we have made dedicated scans of certain regions in order to clarify certain DMSM features. In particular, we gathered dedicated samples of the regions where $2m_{\chi}/m_{Y}$ deviates from unity by $<10^{-3}$, so as to sample adequately annihilations near the $Y$ peak when $\Gamma_{Y}\ll m_{Y}$. We used similar sampling procedure for both the vector and axial-vector DMSMs, generating $\sim 100$ million parameter sets in each case.

4 Results

4.1 Vector DM Couplings

Fig. 4 displays the $(m_{Y},m_{\chi})$ plane in the vector-like model (3) after application of the constraints discussed above. The parameter regions with $\Delta\chi^{2}<2.30$ (5.99), which are favoured at the 68% (95%) CL and regarded as proxies for 1- (2-)$\sigma$ regions, are delineated by red (blue) contours, respectively. In this and subsequent figures, we illustrate the dominant mechanism bringing the DM density into the allowed range at the point that minimizes $\chi^{2}$ in the displayed 2-dimensional projection of the four-dimensional parameter space using colour coding:

$\bullet$ Green: annihilation via $t$-channel $\chi$ exchange into pairs of mediator particles $Y$ that subsequently decay into SM particles, in the region where $m_{\chi}\geq m_{Y}$;

$\bullet$ Yellow: rapid annihilation directly into SM particles via the $s$-channel $Y$ resonance, in the region where $0.9<m_{Y}/(2m_{\chi})<1.1$.

We see clearly two separated favoured regions, one with $m_{\chi}\simeq m_{Y}/2$ where rapid annihilation via the $Y$ funnel dominates, and another with $m_{\chi}>m_{Y}$ where annihilation into pairs of mediator particles $Y$ dominates. The boundaries of both allowed regions are very sharp, reflecting the steepness of the resonance peak in the $s$-channel case and the phase-space limit in the $t$-channel case 555In other projections, the $s$- and $t$-channel regions may overlap, and have quite similar $\chi^{2}$, so that points with different colours appear interspersed..

We see that the $t$-channel region is located at smaller $m_{Y}$ than the $s$-channel region for any value of $m_{\chi}$, subject to the kinematic constraint $m_{Y}<m_{\chi}$. In this region the annihilation cross-section of this channel, $\chi\chi\to YY$, is proportional to $g_{\rm DM}^{4}$ and independent of $g_{\rm SM}$. Therefore, the relic abundance can be brought to the observed value by tuning $g_{\rm DM}$, unless $m_{\chi}\lesssim{100}$ GeV, whilst the experimental constraints from the LHC and the direct DM detection can be avoided by taking $g_{\rm SM}$ small enough.

As can be seen in Fig. 4, the lower limit $m_{Y}>100\,\,\mathrm{GeV}$ in our sample enforces a corresponding lower limit $m_{\chi}\gtrsim 100\,\,\mathrm{GeV}$ in the $t$-channel region. This is because the annihilation process $\chi\chi\to YY$ is kinematically blocked for $m_{\chi}<m_{Y}$, so the cross-section for $\chi\chi$ annihilation becomes very small, resulting in DM overdensity. On the other hand, in the $s$-channel region the relic density constraint requires $g_{\rm DM}$ to be very small for $m_{Y}\sim 100\,\,\mathrm{GeV}$, in which case $\Gamma_{Y}/m_{Y}\ll 1$ and $\chi\chi$ annihilations are rapid enough only if $m_{\chi}\simeq m_{Y}/2\gtrsim 50\,\,\mathrm{GeV}$, the inequality being due to our scanning limit on $m_{Y}$. The uncoloured vertical band at small mass in this and subsequent figures is a reflection of this restriction on the range of $m_{Y}$ sampled.

We find that the $\Delta\chi^{2}$ function is negligible for all values of $m_{Y}$ above the 100 GeV cut, and we also find that the $\Delta\chi^{2}$ function is also very small and essentially featureless for $m_{\chi}>50\,\,\mathrm{GeV}$. Moreover, we find no upper limits on $m_{\chi,Y}$ within the ranges sampled. Our numerical results indicate that

[TABLE]

for intermediate values of $g_{\rm SM}$, with larger values of $\Omega_{\chi}h^{2}$ when $g_{\rm SM}$ is either large or small. Extrapolation of this approximation suggests that the upper bound on $m_{Y}$ in the $s$-channel region is $m_{Y}\sim\mathcal{O}(700)\,\,\mathrm{TeV}$ for $g_{\rm DM}\lesssim\sqrt{4\pi}$.

Fig. 5 uses logarithmic scales to show the $(g_{\rm SM},g_{\rm DM})$ plane in the vector-like model. Again, we distinguish immediately two distinct preferred regions, separated where $g_{\rm SM}\lesssim 3\times 10^{-4}$ and $g_{\rm DM}\gtrsim 0.3$. As can be seen from the colour coding, the region at small $g_{\rm SM}$ and large $g_{\rm DM}$ corresponds to the triangular $t$-channel region in Fig. 4, whereas the larger values of $g_{\rm SM}$  and $g_{\rm DM}$ outside this region appear in the $s$-channel funnel region. This region is bounded at large values of the product g_{\rm SM}$$g_{\rm DM} by the upper limit on $\sigma^{\rm SI}_{p}$ and the limited scanning range $m_{Y}<5\,\,\mathrm{TeV}$, as indicated [see Eq. (15)]: the upper right portion of the figure with larger values of g_{\rm SM}$$g_{\rm DM} would be allowed for larger $m_{Y}$. In the region $0.07>g_{\rm DM}>3\times 10^{-4}$, $\Omega h^{2}\propto m_{Y}^{2}/g_{\rm DM}^{2}$ and can be sufficiently small if $m_{\chi}\simeq m_{Y}/2$ and $m_{Y}<5\,\,\mathrm{TeV}$. In this case, Eqs. (15) and (21) imply

[TABLE]

where $\sigma_{\rm SI}^{\rm UL}(m_{\chi})$ is the mass-dependent experimental upper bound on $\sigma^{\rm SI}_{p}$. We note also that the dijet constraint is the strongest at the largest value of $g_{\rm SM}\sim 0.3$. Thus the lower right part of the figure is excluded by the upper limit on the elastic scattering cross section. Finally, the apparent exclusion in the lower region in the figure, where $g_{\rm DM}\lesssim 3\times 10^{-4}$, rising at smaller $g_{\rm SM}$, is also an artefact of the scanning limit $m_{Y}>100\,\,\mathrm{GeV}$ combined with the relic density constraint, as indicated.

We have indicated by a diagonal dotted line where $g_{\rm DM}=g_{\rm SM}$. As discussed later, we might expect $g_{\rm DM}$ and $g_{\rm SM}$ to be similar in magnitude in many UV completions of DMSMs, and we study later the restriction to the band between the dashed lines in Fig. 5 that is shaded darker yellow, where $1/3<g_{\rm DM}/g_{\rm SM}<3$, see Section 5.

We do not show the one-dimensional likelihood function for $g_{\rm SM}$, which is quite featureless apart from a sharp rise for $g_{\rm SM}\gtrsim 0.3$, nor that for $g_{\rm DM}$, which is also featureless apart from a steep rise for g_{\rm DM}\;\raisebox{-3.00003pt}{\stackrel{{\scriptstyle\displaystyle<}}{{\sim}}}\;3\times 10^{-4} that is an artefact of the limit on our scanning range for $m_{Y}$.

Fig. 6 displays the planes of $m_{Y}$ and the two couplings $g_{\rm SM},\ g_{\rm DM}$. We see in the left panel that within the $t$-channel region there is a strong upper limit on $g_{\rm SM}\lesssim 10^{-3}$, which is enforced by the combination of the upper limit on $\sigma^{\rm SI}_{p}$, which constrains the product $g_{\rm DM}\,g_{\rm SM}$, and the lower limit on $g_{\rm DM}$ visible in the right panel, which is another result of our limited scan to $m_{Y}>100\,\,\mathrm{GeV}$. Since $\sigma^{\rm SI}_{p}$ scales as $g_{\rm DM}^{2}g_{\rm SM}^{2}/m_{Y}^{4}$, the limit is stronger for smaller $m_{Y}$, and $g_{\rm SM}$ must be smaller than $10^{-4}$ for $m_{Y}\lesssim 100$ GeV. The upper limit on $m_{Y}$ in the $t$-channel region visible for small $g_{\rm SM}$ is an artefact of the limited scanning range $m_{\chi}<2.5\,\,\mathrm{TeV}$. The upper limit on $g_{\rm SM}$ in the $s$-channel region comes mainly from the LHC dijet constraint, with the $\sigma^{\rm SI}_{p}$ constraint also playing a role when $m_{Y}\lesssim 500\,\,\mathrm{GeV}$.

We see in the left panel of Fig. 6 that $g_{\rm SM}\lesssim 10^{-2}$ for $m_{Y}\sim 100\,\,\mathrm{GeV}$ and $\lesssim 0.1$ for $m_{Y}>1\,\,\mathrm{TeV}$. Comparing with Fig. 11 of [34], where leptophobic models with universal quark couplings were analyzed 666We recall that in such models the $Z$ and $Y$ have only kinetic mixing, which is loop-induced and hence suppressed., and identifying $g_{\rm SM}$ with the combination $gY^{\prime}_{q}$ of the parameters defined in those models, we see that the precision electroweak data do not constrain our model sample.

In the right panel of Fig. 6, the lower bound on $g_{\rm DM}$ in the $t$-channel region at low $m_{Y}$ is due to the DM density constraint. In the $s$-channel region the relic constraint imposes a weaker bound on $g_{\rm DM}$, which is given roughly by

[TABLE]

as follows from Eq. (21) for $2m_{\chi}\sim m_{Y}$.

Fig. 7 displays the $(m_{\chi},g_{\rm SM})$ plane in the left panel and the $(m_{\chi},g_{\rm DM})$ plane in the right panel. The upper bounds on $g_{\rm SM}$ and the lower bounds on $g_{\rm DM}$ are the same as in Fig. 6, all increasing with $m_{\chi}$.

Fig. 8 shows the likelihood function in the $(m_{\chi},\sigma^{\rm SI}_{p})$ plane for the vector-like model. We see that, as already visible in the left panel of Fig. 7, only values of $m_{\chi}\gtrsim 50$ GeV are allowed, for the reason mentioned previously, namely the interplay of the cut $m_{Y}>100\,\,\mathrm{GeV}$ (as indicated) and the relic density constraint. Above this value, a large range of values of $\sigma^{\rm SI}_{p}$ is allowed in both the $t$- and $s$-channel regions. The upper limits on $\sigma^{\rm SI}_{p}$ at various confidence levels are determined by the combined experimental likelihood for the LUX [45], PANDAX-II [46] and XENON1T [47] experiments, which we have rescaled to account for the different local DM density that we assume. In the $s$-channel region, small values of $\sigma^{\rm SI}_{p}$ are allowed when $m_{\chi}\sim m_{Y}/2$ and small values of $g_{\rm DM}$ and/or $g_{\rm SM}$ are favoured by the relic density constraint, and small values of $\sigma^{\rm SI}_{p}$ are allowed in the $t$-channel region because small values of $g_{\rm SM}$ are allowed, as discussed previously. On the other hand, we see that $\sigma^{\rm SI}_{p}$ may well lie within the range to be probed by upcoming experiments such as LUX-ZEPLIN (LZ) [67] and XENONnT [68], though $\sigma^{\rm SI}_{p}$ may also be much smaller than the current experimental sensitivity, even below the neutrino ‘floor’ indicated by the dashed orange line in Fig. 8, which is based on [51], updating [69].

4.2 Axial-Vector DM Couplings

We now turn to the case of DM with axial-vector couplings. Fig. 9 displays the $(m_{Y},m_{\chi})$ plane in this case: it is also colour-coded according to the dominant DM annihilation mechanism using the same shading scheme as for the vector case, and the 1- and 2-$\sigma$ contours are indicated by red and blue lines, respectively. We see an $s$-channel funnel feature that is rather broader than in the case of DM with vector couplings shown in Fig. 4, and in this case we shade yellow the region where $0.6<m_{Y}/(2m_{\chi})<2$.

This broadening occurs because, whereas the direct detection constraint is very severe for the vector mediator, so that $g_{\rm DM}\,g_{\rm SM}$ cannot be large and the parameters should be near the peak of the resonance where $m_{Y}\sim 2m_{\chi}$, in order for the DM particles to annihilate sufficiently (see Fig. 2 of [26] and the accompanying text), the strong $\sigma^{\rm SI}_{p}$ constraint is absent in the axial-vector case, so that off-resonance regions of parameter space with larger values of $g_{\rm SM}\,g_{\rm DM}$ are allowed where annihilation is $p$-wave or $m_{q}$ suppressed (see the discussion around Eq. (7) of [27]). This opens up more parameter space in the $(m_{Y},m_{\chi})$ plane, with the deviation from $m_{Y}\sim 2m_{\chi}$ bounded only by the dijet constraint. We note that this constraint is weaker when $m_{Y}>2m_{\chi}$ because the decay channel $Y\to\chi\chi$ is open. We also note that, as in the vector case, at low masses the funnel region merges with the $t$-channel annihilation region where $m_{\chi}>m_{Y}$, so that the preferred parameter space is simply connected also in the axial-vector case.

Also as in the vector-like case, we see again in the axial case in Fig. 9 the lower bound $m_{\chi}\gtrsim 50\,\,\mathrm{GeV}$ due to the interplay of the sampling limit $m_{Y}>100\,\,\mathrm{GeV}$ and the relic density constraint. Above this value of $m_{\chi}$, the one-dimensional $\Delta\chi^{2}$ likelihood function for $m_{\chi}$ is featureless, as is that for $m_{Y}>100\,\,\mathrm{GeV}$.

In Fig. 10 we show the $(g_{\rm SM},g_{\rm DM})$ plane in the axial-vector model, using logarithmic scales. The $s$-channel region extends to larger values of $g_{\rm SM}$ and $g_{\rm DM}$ than in the vector case, because of the suppression of axial-vector-mediated $s$-channel annihilation discussed above. Values of $g_{\rm DM}$ as large as the sampling limit $\sqrt{4\pi}$ are allowed. We see a separation between the regions of this plane where the $s$- and $t$-channel mechanisms are dominant when $g_{\rm SM}\sim 10^{-2}$ and $g_{\rm DM}\sim 10^{-1}$. The upper bound on $g_{\rm SM}$ comes from the dijet constraint. As in the vector case, the dotted line is where $g_{\rm SM}=g_{\rm DM}$, and the deeper shading indicates where $1/3<g_{\rm SM}/g_{\rm DM}<3$, favouring the large-$g_{\rm SM}$ part of the $s$-channel annihilation region, see the discussion in Section 5. As in the vector case, there is a region at low $g_{\rm DM}$, rising at small $g_{\rm SM}$, whose exclusion by the relic density constraint is an artefact of the sampling restriction $m_{Y}>100\,\,\mathrm{GeV}$.

We display in Fig. 11 the $(m_{Y},g_{\rm SM})$ and $(m_{Y},g_{\rm DM})$ planes (left and right panels, respectively) in the scenario with axial-vector couplings. We see in the left panel that the $t$-channel DM mechanism is important for $m_{Y}\lesssim 2.6\,\,\mathrm{TeV}$. The limit visible in Fig. 9, which is due to the sampling limit $m_{\chi}<2.5\,\,\mathrm{TeV}$. Annihilation via the $s$-channel becomes more important as $m_{Y}$ increases, as also seen in Fig. 9, and is the only mechanism for $m_{Y}\gtrsim 2.6\,\,\mathrm{TeV}$ up to the sampling limit of $5\,\,\mathrm{TeV}$. The upper limit on $g_{\rm SM}$ for $m_{Y}<5\,\,\mathrm{TeV}$ is due to the dijet constraint. Unlike the vector case seen in the left panel of Fig. 6, the $s$-channel region is excluded for small $g_{\rm SM}$. This is due to the $p$-wave nature of the $s$-channel annihilation via the axial mediator, which leads to the scaling behaviour \Omega h^{2}\mathrel{\raisebox{3.01385pt}{\mbox{\underset{\vbox{\hbox{$\sim$}}}{\vbox{\hbox{$\propto$}}}}}}m_{Y}^{2}/(g_{\rm SM}^{2}g_{\rm DM}^{2}), as discussed around Eq. (9).

We see again in the right panel that the $t$-channel mechanism is relevant for $m_{Y}\lesssim 3\,\,\mathrm{TeV}$, whereas the $s$-channel mechanism is dominant for larger $m_{Y}$. The lower limit on $g_{\rm DM}$ is given by the relic density constraint since. even exactly at $m_{Y}=2m_{\chi}$, \Omega h^{2}\mathrel{\raisebox{3.01385pt}{\mbox{\underset{\vbox{\hbox{$\sim$}}}{\vbox{\hbox{$\propto$}}}}}}m_{Y}^{2}/(g_{\rm SM}^{2}g_{\rm DM}^{2}) and DM is overproduced for sufficiently small $g_{\rm DM}$.

Fig. 12 shows the $(m_{\chi},g_{\rm SM})$ and $(m_{\chi},g_{\rm DM})$ planes (left and right panels, respectively) in the axial-vector model. We see again in the left panel that values of $m_{\chi}$ up to the sampling limit of $2.5\,\,\mathrm{TeV}$ are allowed both at smaller $g_{\rm SM}$ where the $t$-channel mechanism dominates and at larger $g_{\rm SM}$ where the $s$-channel mechanism dominates. There is a region at low $m_{\chi}$ where the $s$-channel mechanism dominates. The lower limit on $g_{\rm DM}$ comes from the relic density constraint as can be seen in Eqs. (12), (13) and (14).

The one-dimensional $\Delta\chi^{2}$ function for $g_{\rm DM}$ in the axial case rises sharply below $\sim 4\times 10^{-3}$, and that for $g_{\rm SM}$ rises above $\sim 0.3$. The $\Delta\chi^{2}$ functions for $m_{\chi}$ and $m_{Y}$ are featureless above 50 and 100 GeV, respectively.

Finally, in Fig. 13 we compare the experimental constraints, ranges favoured in our axial-vector DMSM analysis and prospective experimental sensitivities to the cross section for spin-dependent scattering on a proton ($\sigma^{\rm SD}_{p}$, inferred from the PICO-60 search with a C3F8 target) [48] (left panel) and to that on a neutron ($\sigma^{\rm SD}_{n}$, inferred from a search with the XENON1T detector) [49] (right panel), again accounting for the different local DM density that we assume. Since we consider here leptophobic mediators, the constraints on $\sigma^{\rm SD}_{p}$ provided by Super-Kamiokande [70] and IceCube [71] limits on the annihilations into $\tau^{+}\tau^{-}$ of DM particles trapped in the Sun are not relevant, and the limits on hadronic annihilations are not competitive with the direct constraints on $\sigma^{\rm SD}_{p}$. We also show in the left panel the estimated neutrino ‘floor’ applicable to experiments using a C3F8 target (shaded blue), adapted from [72] using a similar factor as in [51]. We see that $\sigma^{\rm SD}_{p}$ approaches the PICO-60 limit most closely for a small range $100~{}{\rm GeV}\lesssim m_{\chi}\lesssim 200$ GeV, and that $\sigma^{\rm SD}_{p}$ may be accessible to the LZ experiment [67] or the PICO-500 experiment [72] for $m_{\chi}\lesssim 1\,\,\mathrm{TeV}$. However, we see in the right panel that the most sensitive limit on spin-dependent scattering is currently set by the XENON1T experiment, which is sensitive to $\sigma^{\rm SD}_{n}$, and that the LZ experiment may be able to increase this sensitivity significantly. We also display the neutrino ‘floor’ for an experiment using a Xenon target [72] (shaded blue). Encouragingly, we note that there are significant regions of the axial-vector parameter space where both $\sigma^{\rm SD}_{p}$ and $\sigma^{\rm SD}_{n}$ may be detectable above the corresponding neutrino ‘floors’, in both the $s$- and the $t$-channel regions. However, we note that, whereas the short-term advantage may lie with searches for $\sigma^{\rm SD}_{n}$, since that currently provides the stronger constraint, the longer-term advantage may lie with searches for $\sigma^{\rm SD}_{p}$, since the expected ‘floor’ is lower in that case.

.

5 Possible Ultraviolet Completions

The vector- and axial-like leptophobic DMSMs analyzed here were chosen following the recommendations of the LHC Dark Matter Working Group [5, 6, 7, 8, 9], on the basis of their phenomenological simplicity and without regard for their possible ultraviolet (UV) completions. In any such UV completion, the spin-one boson could be expected to have comparable couplings to SM and DM particles, modulo possible group-theoretical factors and mixing angles.

Looking beyond the requirements of specific grand unified or string scenarios, one should consider the important consistency conditions on gauge couplings that are imposed by the cancellation of anomalous triangle anomalies required for renormalizability. These entail, characteristically, that the gauge couplings to different particle species are related by rational algebraic factors that are ${\cal O}(1)$ before mixing. Supersymmetry is an important example of a framework where such mixing factors are important, but the couplings of supersymmetric WIMPs to the SU(2)$\times$U(1) gauge bosons of the SM are typically not much smaller than those of SM particles.

The construction of DMSMs with anomaly-free U(1)′ gauge bosons has been studied in [33], where it was found that in order to be leptophobic, such DMSMs would necessarily contain non-trivial dark sectors containing other particles besides the DM particle. Explicit examples have been given of leptophobic DMSMs with purely vector- or axial-like couplings to quarks, and DMSMs with purely vector- or axial-like couplings to the DM particle would be subject to further constraints. We do not discuss the construction of such models here, but expect on the basis of the argument in the previous paragraph that they would in general have $g_{\rm DM}/g_{\rm SM}={\cal O}(1)$. Additionally, one might expect that in a UV completion featuring unification in a non-Abelian gauge group the spin-1 mediator couplings would be $\gtrsim{\cal O}(0.1)$, favouring parts of the funnel regions away from the regions where $t$-channel annihilation is dominant.

We have included in Fig. 5 and 10 diagonal dotted lines where $g_{\rm DM}=g_{\rm SM}$ and shaded the strips where $1/3<g_{\rm SM}/g_{\rm DM}<3$, which are bounded by dashed lines. We see that in the vector case it traverses the funnel region where the DM particles annihilate via the $s$-channel, and does not approach the regions where DM particles annihilate mainly via $t$-channel exchanges into pairs of mediator particles, which are dominant when $g_{\rm DM}\gg g_{\rm SM}$. On the other hand, in the axial-vector case both $t$- and $s$-channel annihilations are possible in the strip with $1/3<g_{\rm SM}/g_{\rm DM}<3$, as we discuss later. We now discuss in more detail the implications of the constraint $1/3<g_{\rm SM}/g_{\rm DM}<3$ for the DMSM parameter spaces.

We start by displaying in Fig. 14 various parameter planes for the vector-coupling case with the selection $1/3<g_{\rm SM}/g_{\rm DM}<3$. We see that both $g_{\rm SM}$ and $g_{\rm DM}$ are $>10^{-3}$ when this selection is imposed, and that values of $g_{\rm SM}$ and $g_{\rm DM}\lesssim 0.3$ lie within the favoured range, consistent with unification scenarios. The ranges $m_{Y}\gtrsim 2\,\,\mathrm{TeV}$ and $m_{\chi}\gtrsim 1\,\,\mathrm{TeV}$ are compatible with $g_{\rm SM},g_{\rm DM}>0.1$. We note that the lower bound comes from the DM density condition, which requires $s$-channel annihilation to be fast enough. The upper bounds on $g_{\rm SM}$ and $g_{\rm DM}$ are, on the other hand, largely due to the spin-independent DM-nucleus scattering constraint, with the dijet constraint also playing a role in constraining $g_{\rm SM}$ for $m_{Y}\gtrsim 1.8$ TeV ($m_{\chi}\gtrsim 0.9$ TeV). These features are also visible in the region of Fig. 5 with the darker yellow band.

Fig. 15 shows the corresponding parameter planes for the axial-coupling case with the selection $1/3<g_{\rm SM}/g_{\rm DM}<3$. In the upper panels, we see that the $t$-channel region is confined to a small region where $m_{Y}\lesssim 500$ GeV and $g_{\rm SM}\gtrsim 0.1$. We see from the right panel of Fig. 11 that $g_{\rm DM}$ must be larger than $\sim 0.3$ (0.9) at $m_{Y}=500$ GeV (2 TeV) to account for the DM density. The $1/3<g_{\rm SM}/g_{\rm DM}<3$ condition then implies that $g_{\rm SM}\gtrsim 0.1$ (0.3) for $m_{Y}=500$ GeV (2 TeV). However, such a large value of $g_{\rm SM}$ is not compatible with the dijet constraint in this range of $m_{Y}$. Therefore, no region with $m_{Y}\gtrsim 500$ GeV is allowed in the $t$-channel annihilation region. As discussed above, $s$-channel annihilation undergoes a $p$-wave or $m_{q}$ suppression in the axial-vector mediator case, and $g_{\rm SM}\gtrsim 10^{-3}$ and $g_{\rm DM}\gtrsim 3\times 10^{-3}$ (see the darker yellow band in Fig. 10), are needed to compensate this suppression. However, there is a narrow band of $g_{\rm SM}$ values that is also compatible with the dijet constraint for all the sampled range of $m_{Y}$. We see that the $s$-channel region is allowed for all the sampled range of $m_{Y}$, for the most part also if $g_{\rm SM}>0.1$.

In the lower panels of Fig. 15 we see the corresponding $(m_{\chi},g_{\rm SM})$ and $(m_{\chi},g_{\rm DM})$ planes. We see in the left plot that in this projection the $t$-channel region appears in the coupling range $0.03\lesssim g_{\rm SM}\lesssim 0.3$, and in the right plot its range is $0.05\lesssim g_{\rm DM}\lesssim 0.5$, sandwiched by the DM density condition from below and the LHC dijet constraint from above. As seen in this projection, the $t$-channel region is less restricted in the $m_{\chi}$ direction, and any value of $m_{\chi}\gtrsim 300\,\,\mathrm{GeV}$ is in principle compatible with this mechanism. In the $s$-channel region all values of $m_{\chi}$ between 50 GeV and 2.5 TeV are allowed. We note also that $g_{\rm DM},g_{\rm SM}>0.1$ are possible for all the range of $m_{\chi}$.

In the vector-like case, the main changes in the one-dimensional $\Delta\chi^{2}$ functions for $g_{\rm DM}$ and $g_{\rm SM}$ due to the selection $1/3<g_{\rm SM}/g_{\rm DM}<3$ are reductions in their upper limits to $\sim 0.3$, and there are negligible changes in the one-dimensional $\Delta\chi^{2}$ functions for $m_{Y}$ and $m_{\chi}$. In the axial case the main changes when the selection is made are that $g_{\rm DM}\lesssim 0.5$ and $g_{\rm SM}\gtrsim 10^{-3}$, while all values of $m_{\chi}\gtrsim 50\,\,\mathrm{GeV}$ and $m_{Y}\gtrsim 100\,\,\mathrm{GeV}$ still have very low $\Delta\chi^{2}$.

Finally, we show in Fig. 16 the preferred ranges of $(m_{\chi},\sigma^{\rm SI}_{p})$ in the vector-like model (left panel) and $(m_{\chi},\sigma^{\rm SD}_{n})$ in the axial model (right panel) after the selection $1/3<g_{\rm SM}/g_{\rm DM}<3$. We see that there are good prospects for detecting spin-independent DM scattering in the vector-like case, since $\sigma^{\rm SI}_{p}$ lies partially above the neutrino ‘floor’ within the range of $m_{\chi}$ sampled, though smaller values of $\sigma^{\rm SI}_{p}$ could also occur for any value of $m_{\chi}>50\,\,\mathrm{GeV}$. In the case of spin-dependent scattering after the selection $1/3<g_{\rm SM}/g_{\rm DM}<3$, we see in the right panel of Fig. 16 that in the $t$-channel exchange region the predicted values of $\sigma^{\rm SD}_{n}$ are relatively close to the XENON1T limit [49], and hence may offer prospects for detection with the planned LZ experiment [67] within the sampled range of DM masses. However, values of $\sigma^{\rm SD}_{n}$ may be considerably lower in the $s$-channel exchange region 777Note that we have again not included in the right panel the upper limits and prospective detection ‘floors’ presented by the Super-Kamiokande [70] and IceCube [71] Collaborations, which are based on specific assumptions about the DM annihilation channels that are not applicable in the leptophobic DMSMs studied here.. This can be traced to features visible in Fig. 15: we see in the upper panels that the $t$-channel region corresponds to small values of $m_{Y}$ with larger values of $g_{\rm SM}$ and $g_{\rm DM}$ than in the $s$-channel region, and in the lower panels these projections show that the $t$-channel region may have larger values of $g_{\rm SM}$ and $g_{\rm DM}$ than in the $s$-channel region for a large range of values of $m_{\chi}$. The intermediate ‘hole’ is enforced by the LHC dijet constraint, in particular.

It is interesting that, with the selection $1/3<g_{\rm SM}/g_{\rm DM}<3$ motivated by the prospect of UV completion, searches for $\sigma^{\rm SI}_{p}$ and $\sigma^{\rm SD}_{n}$ are able to probe both the $s$- and $t$-channel regions of the DMSM parameter spaces

6 Conclusions and Perspectives

In this paper we have used MasterCode to make a global analysis of the parameter spaces of dark matter simplified models with leptophobic spin-one mediator particles $Y$ with either vectorial or axial couplings to SM particles and to the dark matter particle $\chi$. Each of these models is characterized by four free parameters: $m_{Y}$ and $m_{\chi}$, the coupling $g_{\rm DM}$ of the mediator to the dark matter particle, and the coupling $g_{\rm SM}$ of the mediator particle to quarks, which we assume to be independent of flavour.

We have implemented constraints on the model parameter spaces due to LHC searches for monojet events and measurements of the dijet invariant mass spectrum, as well as the cosmological constraint on the dark matter density and direct upper limits on spin-independent and -dependent scattering on nuclei. We have scanned mediator masses $m_{Y}\leq 5\,\,\mathrm{TeV}$ and dark matter particle masses $m_{\chi}\leq 2.5\,\,\mathrm{TeV}$, delineating the regions of the model parameters with $\Delta\chi^{2}<2.30(5.99)$, which are favoured at the 68% (95%) CL and regarded as proxies for 1- (2-)$\sigma$ regions, respectively. Within these regions we have identified two main mechanisms for bringing the cosmological dark matter density into the range allowed by cosmology, namely annihilation via $t$-channel $\chi$ exchange and annihilation via the $Y$ boson in the $s$ channel. With an eye to possible ultraviolet completions of the simplified models studied here, we have also explored the portions of the favoured parameter space where $1/3<g_{\rm DM}/g_{\rm SM}<3$, as discussed below.

In the vector-like case, we find a relatively clear separation between the regions where the $t$- and $s$-channel mechanisms dominate, with the former being more important at smaller mediator masses, small values of $g_{\rm SM}$ and relatively large values of $g_{\rm DM}$. The one-dimensional likelihood functions for both $m_{\chi}$ and $m_{Y}$ are quite small and flat above thresholds $\sim 50\,\,\mathrm{GeV}$ and $\sim 100\,\,\mathrm{GeV}$. Thus the LHC still has interesting prospects for discovering DM and mediator particles in these simplified models. Any value of $g_{\rm SM}$ between $\sim 10^{-6}$ and the dijet limit of $\sim 0.3$ is possible without any $\chi^{2}$ penalty, as is the case for $g_{\rm DM}\gtrsim 10^{-3}$, where the lower limit is due to the upper limit on the sampling range for $m_{Y}$. The spin-independent dark matter scattering cross section $\sigma^{\rm SI}_{p}$ may be very close to the present experimental upper limit, within the range accessible to the upcoming LZ and XENONnT experiments. However, in both the $s$- and the $t$-channel cases $\sigma^{\rm SI}_{p}$ may also be much below the neutrino ‘floor’.

In the axial case the $t$- and $s$-channel regions are more connected. The one-dimensional likelihood functions for $m_{\chi}$ and $m_{Y}$ are again featureless above thresholds $\sim 50\,\,\mathrm{GeV}$ and $100\,\,\mathrm{GeV}$, respectively, and that for $g_{\rm SM}$ is quite featureless, as is that for $g_{\rm DM}\gtrsim 10^{-2}$. We find that the spin-dependent dark matter scattering cross sections $\sigma^{\rm SD}_{p}$ and $\sigma^{\rm SD}_{n}$ may also be very close to the present experimental upper limit, within the range accessible to the upcoming PICO-500 and LZ experiments, though much lower values below the corresponding ‘floors’ are also possible.

Finally, we have also explored the possibility that $1/3<g_{\rm DM}/g_{\rm SM}<3$, as might be suggested in some ultraviolet completions of the dark matter simplified models considered here. In this case, values of $g_{\rm SM},g_{\rm DM}\sim 0.1$ are possible, as as might also be suggested in some scenarios with unified gauge interactions, and wide ranges of DM and mediator masses are again accessible to the LHC experiments. In the case of vector-like couplings, this reduced range of $g_{\rm DM}/g_{\rm SM}$ disfavours the $t$-channel mechanism. However, in the axial case the $t$-channel mechanism is still possible, yields values of $\sigma^{\rm SD}_{n}$ that are relatively close to the upper limit set by the PICO-60 experiment, and may well be within reach of the upcoming PICO-500 and LZ experiments, whereas lower values of $\sigma^{\rm SD}_{n}$ are possible if the $t$-channel mechanism dominates.

In this paper, we have performed a first global study - including all the relevant constraints from Run 2 of the LHC and direct scattering searches and not fixing any of the underlying parameters - of some DMSMs proposed by the LHC Dark Matter Working Group [5, 6, 7, 8, 9] using the MasterCode framework. Our analysis also provides a more detailed study of the vector and axial leptophobic models considered in [34]. These models are undoubtedly over-simplified, and it would be interesting and useful to extend this type of analysis to other models that may be more realistic, which generally contain more parameters. For example, one should study models with spin-one mediators that are not leptophobic, and also models with spin-zero mediators that may be either scalar or pseudoscalar. One could also study models that are not flavour-universal in the quark sector, which may be motivated by the anomalies reported in $B$ meson decays. Consistent ultraviolet completions of DMSMs should include a mechanism for anomaly cancellation, which typically include more ‘dark’ particles whose possible phenomenological signatures could also be considered. The MasterCode tool is a very suitable tool for such analyses, and we plan to use MasterCode for such analyses in the future.

Acknowledgements

We would like to thank Y. Mambrini and M. Voloshin for helpful discussions. G.W. acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306. The work of K.S. has been partially supported by the National Science Centre, Poland, under research grants DEC-2014/15/B/ST2/02157, DEC-2015/18/M/ST2/00054 and DEC- 2015/19/D/ST2/03136. The work of M.B. and D.M.S. has been supported by the European Research Council via Grant BSMFLEET 639068. The work of J.C.C. is supported by CNPq (Brazil). The work of M.J.D. is supported in part by the Australia Research Council. The work of J.E. is supported in part by STFC (UK) via the research grant ST/L000258/1 and in part via the Estonian Research Council via a Mobilitas Pluss grant, and the work of H.F. is also supported in part by STFC (UK) via grant ST/N000250/1. The work of S.H. is supported in part by the MEINCOP Spain under contract FPA2016-78022-P, in part by the Spanish Agencia Estatal de Investigación (AEI) and the EU Fondo Europeo de Desarrollo Regional (FEDER) through the project FPA2016-78645-P, in part by the AEI through the grant IFT Centro de Excelencia Severo Ochoa SEV-2016-0597, and by the ”Spanish Red Consolider Multidark” FPA2017-90566-REDC. The work of M.L. is supported by XuntaGal. The work of K.A.O. is supported in part by DOE grant desc0011842 at the University of Minnesota. During part of this work, to deploy MasterCode on clusters we used the middle-ware suite udocker [64], which was developed by the EC H2020 project INDIGO-Datacloud (RIA 653549). We also thank Imperial College London and the University of Bristol for making available to us cluster computing resources that have been used intensively to carry out this work.