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,3,4,5]. They have also suggested that DM is mainly cold, i.e., it has been non-relativistic since at least the epoch of cosmological recombination [6, 7]. 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. [8] 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 \(\mathcal{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  [9,10,11,12,13].

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) [14, 15] and some other models with new physics at the TeV scale [16, 17]. 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 [18,19,20].

We have used the MasterCode framework [21] previously to study several SUSY scenarios with various experimental signatures [22,23,24,25,26,27,28,29,30,31,32], and other groups have made parallel analyses [33,34,35]. In particular, in Refs. [22,23,24,25,26,27,28,29,30,31,32] 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 [36,37,38,39,40], 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 [41, 42]. For this reason, interest has been growing in Dark Matter Simplified Models (DMSMs) [42,43,44,45,46,47,48,49,50], 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) [9,10,11,12,13].

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 [51, 52], 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

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

and that for its interactions with quarks is

$$\begin{aligned} \mathcal{L}^1_q \; = \; \sum _{i} \bar{q_i} \gamma _\mu \left( g_{q_i}^V + g_{q_i}^A \gamma _5 \right) q_i Y^\mu \, . \end{aligned}$$
(2)

For simplicity, following the DMWG [12] 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\)).Footnote 1 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:

$$\begin{aligned} g_\chi ^V \; \equiv \; g_\mathrm{DM}\; \ne \; 0,&g_\chi ^A \; = \; 0, \nonumber \\ g_q^V \; \equiv \; g_\mathrm{SM}\; \ne \; 0,&g_q^A \; = \; 0, \end{aligned}$$
(3)

and another in which they are purely axial:

$$\begin{aligned} g_\chi ^A \; \equiv \; g_\mathrm{DM}\; \ne \; 0,&g_\chi ^V \; = \; 0, \nonumber \\ g_q^A \; \equiv \; g_\mathrm{SM}\; \ne \; 0,&g_q^V \; = \; 0. \end{aligned}$$
(4)

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 [47, 62].

UV completions of the DMSMs we study could be achieved in a number of different ways. These could favour specific ranges of \(g_\mathrm{SM}\), \(g_\mathrm{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 [63]. 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_\mathrm{DM}\) and \(g_\mathrm{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_\mathrm{DM}\) and \(g_\mathrm{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 5.0.2 [64] and MG5_aMC 2.4.3 [65], which we run at LO for performance reasons, using DMSIMP model files [66, 67] to build a likelihood function for each of the astrophysical and accelerator constraints, and combine them using the MasterCode framework. We use MultiNest 2.18 [68,69,70] 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 [55] to be unimportant in simple DM models when \(m_\chi > rsim 50 \,\, \mathrm {GeV}\).

The layout of this paper is as follows. In Sect. 2 we set out the constraints we use, and in Sect. 3 we describe our analysis framework. Our results are described in Sect. 4, Possible UV completions are discussed in Sect. 5, and our conclusions and some perspectives are set out in Sect. 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 [71]:

$$\begin{aligned} \Omega _\mathrm{CDM} h^2 \; = \; 0.120 \pm 0.001 \, . \end{aligned}$$
(5)

We assume here that the dominant source of this cold DM is the WIMP \(\chi \), so that \(\Omega _\chi h^2\simeq \Omega _\mathrm{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 \rightarrow Y^{(*)} \rightarrow \) Standard Model particles, the leading expression for the cross-section takes the form [72, 73]

$$\begin{aligned} \langle \sigma v \rangle _s \simeq \left\langle \frac{c_Y}{32 \pi } \frac{g_\mathrm{SM}^2 g_\mathrm{DM}^2 m^2_\chi }{(m_\chi ^2 v^2 + 4 m_\chi ^2 - m_Y^2)^2 + m_Y^2 \Gamma _Y^2}\right\rangle , \end{aligned}$$
(6)

where \(c_Y\) is an \(\mathcal{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} + b v^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

$$\begin{aligned} \Gamma _Y^{XX} \; = \; \frac{g_X^2 m_Y}{12 \pi } \cdot f \left( \frac{m_X^2}{m_Y^2} \right) \, , \end{aligned}$$
(7)

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

$$\begin{aligned} f(x) \; = \; {\left\{ \begin{array}{ll} (1 + 2 x) (1 - 4 x)^{1/2} &{}\quad (\mathrm{vector}) \\ (1 - 4 x)^{3/2} &{}\quad (\mathrm{axial\,vector}) \end{array}\right. }, \end{aligned}$$
(8)

where \(x \equiv (m_X/m_Y)^2\). The width of Y becomes particularly important in the resonant region, where \(m_Y \simeq 2 m_\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_\mathrm{SM}\ll 1\), the velocity term dominates the denominator of Eq. (6) at the resonance. This leads to the approximate scaling

(9)

where 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_\mathrm{SM}\ll 1\), and one finds that in the limit \(g_\mathrm{SM}\rightarrow 0\)

$$\begin{aligned}&1/[(m_Y^2 - 4 m_\chi ^2)^2 + m_Y^2 \Gamma _Y^2] \nonumber \\&\quad \rightarrow \delta (m_Y^2 - 4 m_\chi ^2)/(m_Y \Gamma _Y), \end{aligned}$$
(10)

so the scaling of the relic density is expected to be

(11)

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 \(\ge 2 m_\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_\mathrm{DM}\) dependence of \(\Omega _\chi h^2\), and that it depends only weakly on \(g_\mathrm{SM}\) over an intermediate range of \(g_\mathrm{SM}\), but that \(\Omega _\chi h^2\) increases at both large and small \(g_\mathrm{SM}\) :

(12)

for \(g_\mathrm{SM}\gg g_\mathrm{DM}\) and

(13)

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

Another important process is the t-channel process \(\chi \chi \rightarrow Y Y\) for \(m_\chi > m_Y\). The leading contribution to the t-channel cross-section is given by [72]

$$\begin{aligned} \langle \sigma v \rangle _t \simeq \frac{g_\mathrm{DM}^4}{32 \pi } \frac{(m_\chi ^2 - m_Y^2)^{3/2} }{m_\chi (2 m_\chi ^2 - m_Y^2 )^2}, \end{aligned}$$
(14)

and is kinematically forbidden when \(m_\chi < m_Y\).

We use MicrOMEGAs 5.0.2 [64] 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 5.0.2 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 [75] (dashed line) for the case of a spin-one mediator with vector-like couplings \(g_\mathrm{SM}^V = 0.25, g_\mathrm{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 \).

Fig. 1
figure 1

Comparison of our calculation 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 [75] (dashed white line) for the DMSM with a spin-one mediator having vector-like couplings \(g_\mathrm{SM}^V = 0.25\), \(g_\mathrm{DM}^V = 1\) to quarks and to DM particles \(\chi \), respectively

2.1.2 Spin-independent DM scattering

We evaluate the constraint on the spin-independent DM scattering cross-section \(\sigma ^\mathrm{SI}_p\) from a combination of the LUX [76], PandaX-II [77] and XENON1T [78] experiments using a joint likelihood function, evaluating the theoretical prediction for \(\sigma ^\mathrm{SI}_p\) using MicrOMEGAs 5.0.2. This constraint is important for the vector mediator case, where the spin-independent cross-section has the following scaling [11] \(\sigma ^\mathrm{SI}_p\propto g_\mathrm{DM}^2 g_\mathrm{SM}^2/m_Y^4\), with the following approximate numerical value:

$$\begin{aligned} \sigma ^\mathrm{SI}_p\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} \end{aligned}$$
(15)

in the s-channel region where \(m_Y \simeq 2 m_\chi \). The \(\sigma ^\mathrm{SI}_p\) constraint provides important limits on the vector mediator mass \(m_Y\) and the product \(g_\mathrm{DM}g_\mathrm{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 ^\mathrm{SD}_p\) from the PICO-60 experiment [79], and that on the spin-dependent DM-neutron scattering cross-section \(\sigma ^\mathrm{SD}_n\) from the XENON1T experiment [80]. We evaluate theoretical predictions for \(\sigma ^\mathrm{SD}_p\) and \(\sigma ^\mathrm{SD}_n\) using MicrOMEGAs 5.0.2, which uses estimates of the spin-dependent scattering matrix elements consistent with [81].

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/cm\(^3\), with negligible uncertainty. However, several recent analyses favour a larger local density, and here we follow [82] in assuming a central value for the local DM density of 0.55 GeV/cm\(^3\). 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 [83]. However, this modification of 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 [82]. For another estimate of the local density of DM based on Gaia data, see [84].

When implementing the experimental constraints on \(\sigma ^\mathrm{SI}_p\) and \(\sigma ^\mathrm{SD}_p, \sigma ^\mathrm{SD}_n\), we assume an overall uncertainty of 30% [82, 85], in the local density.

2.2 LHC experimental constraints

2.2.1 Monojet constraints

Spin-1 mediators decaying to missing energy are constrained by analyses targeting monojet final states by ATLAS [74] and CMS [75]. In this analysis we use the result obtained by CMS Collaboration [75] using 35.9/fb of data from collisions at 13 TeV. The ATLAS analysis of the axial coupling case has similar sensitivity. In particular, we use the 95 % CL upper bounds \(R_i^\mathrm{UL}(\mathbf{m})\), where i labels the vector and axial-vector mediators, imposed on \(\sigma (\mathbf{m}) / \sigma ^\mathrm{fix}(\mathbf{m})\) evaluated at each point in the \(\mathbf{m} = (m_\chi , m_Y)\) plane, where \(\sigma ^\mathrm{fix}\) is the cross-section evaluated with \((g_\mathrm{SM}, g_\mathrm{DM}) = (0.25, 1)\).

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

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

where \(\sigma _\mathrm{MG}\) is the monojet cross-section calculated using MG5_aMC 2.4.3 [65].

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 [75], which are indicated by dashed red lines.

Fig. 2
figure 2

The likelihoods for our implementations (16) (colours and solid red lines) of the CMS monojet constraints [75] (dashed green lines) for a spin-one mediator having (left panel) vector-like couplings and (right panels) axial couplings to quarks and DM particles

2.2.2 Dijet constraints

The ATLAS and CMS Collaborations have both published constraints from dijet invariant-mass distributions on \(Z'\) resonances [87,88,89,90,91,92] with up to 139/fb of data at 13 TeV [93]. We use a combination of the limits presented in the plane of the \(Z'\) mass and its \({{\bar{q}}} q\) coupling, \(g_q\). In recasting this constraint for the \(pp \rightarrow Y \rightarrow qq\) process in our DM models, we demand

$$\begin{aligned} \sigma _{Y \rightarrow qq} \; \le \; \sigma (g_q^*) \end{aligned}$$
(17)

at \(m_Y = m_{Z'}\), where \(g_q^*\) is the upper limit placed on the coupling of the \(Z'\) to quarks at a given \(m_{Z'}\).Footnote 2 We note that \(\sigma _{Y \rightarrow qq}\) and \(\sigma ({g_q^*})\) can be written as (using \(\Gamma _x \equiv \Gamma (Y \rightarrow xx)\))

$$\begin{aligned} \sigma _{Y \rightarrow qq} = c \, \frac{g_\mathrm{SM}^4}{\Gamma _{q} + \Gamma _{\chi }}, \quad \sigma (g_q^*) = c \, \frac{(g_q^*)^4}{\Gamma ^*_{q}} \end{aligned}$$
(18)

with the same constant c at \(m_Y = m_{Z'}\), and \(\Gamma _q = \big ( \frac{g_\mathrm{SM}}{g_q^*} \big )^2 \Gamma _q^*\). Using these relations, the equality in Eq. (17) occurs when \(g_\mathrm{SM}\) takes the value

$$\begin{aligned} g_\mathrm{SM}^* \; \equiv \; \sqrt{\frac{(g^*_q)^2}{2} \left( 1 + \sqrt{1 + \frac{4 \Gamma _{\chi }}{\Gamma ^*_q}} \right) }, \end{aligned}$$
(19)

which serves as the upper bound on \(g_\mathrm{SM}\) for given \(m_Y\), \(\Gamma _q\) and \(\Gamma _{\chi }\). We model the corresponding likelihood function using

$$\begin{aligned} \Delta \chi ^2 = 4 \times \left( \left[ \frac{g_\mathrm{SM}^4}{\Gamma _{q} + \Gamma _{\chi }} \right] / \left[ \frac{(g_\mathrm{SM}^*)^4}{\Gamma _{q}(g_\mathrm{SM}^*) + \Gamma _{\chi }} \right] \right) ^2 \end{aligned}$$
(20)

where \(\Gamma _{q}\) and \(\Gamma _{\chi }\) are obtained with MG5_aMC 2.4.3.

As a further step in our validation, we have cross-checked our recasting with the CMS and ATLAS limits [87,88,89,90,91,92,93], as shown in Fig. 3.Footnote 3 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_\mathrm{SM}\) and that quoted by the experiments is a valuable cross-check of our implementation of the dijet invariant-mass constraints.

Fig. 3
figure 3

Comparison of our implementation of the LHC dijet invariant-mass constraint, indicated by the colour coding, with a compilation of the limits published by the ATLAS and CMS Collaborations  [87,88,89,90,91,92,93], for a \(Z^{'}\) boson, in the plane of its mass and coupling to the quarks

In addition to these invariant-mass searches, the CMS Collaboration has also published a limit on \(g_\mathrm{SM}\) as a function of \(m_{Z'}\) from an analysis of dijet angular distributions [94]. 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. We also do not consider the \(Y \rightarrow t {\bar{t}}\) mode, since the sensitivity is weaker than in the dijet mode, due to the large \(t {\bar{t}}\) background.

Please note also that in evaluating the dijet and monojet constraints on model couplings we work only to leading order in QCD, in the absence of a complete set of the model-dependent NLO corrections. We expect this treatment to be conservative, in the sense that these corrections are generally positive, so that the upper limits on couplings that we obtain would be strengthened if NLO corrections were included.Footnote 4

3 Analysis framework

Global-fitting frameworks such as MasterCode [22,23,24,25,26,27,28,29,30,31,32] 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_aMC 2.4.3 [65] to calculate mediator properties and the collider constraints on models implemented using DMSIMP [66, 67], and we use MicrOMEGAs 5.0.2 [64] for DM density and scattering calculations. We use the MultiNest 2.18 [68,69,70] 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 [97] to automatize the deployment of MasterCode inside Linux containers. This is a middleware suite developed in the context of the INDIGO datacloud project [98] 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 dataFootnote 5 and indirect searches for astrophysical DM annihilations should be considered,Footnote 6 but include all the masses for which LHC searches are sensitive. The couplings \(g_\mathrm{SM}\) and \(g_\mathrm{DM}\) were restricted to perturbative ranges \(< \sqrt{4 \pi }\).

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

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

Figure 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:

  •  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 \ge m_Y\);

  •  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.Footnote 7

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 \rightarrow Y Y\), is proportional to \(g_\mathrm{DM}^4\) and independent of \(g_\mathrm{SM}\). Therefore, the relic abundance can be brought to the observed value by tuning \(g_\mathrm{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_\mathrm{SM}\) small enough.

Fig. 4
figure 4

Preferred regions in the \((m_Y, m_\chi )\) plane in the vector-like model. We delineate with red (blue) contours, respectively, 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, respectively. We also use colour coding to illustrate the dominant mechanisms bringing the DM density into the allowed range: green for annihilation via t-channel exchange into pairs of mediator particles Y that subsequently decay into SM particles, and yellow for rapid annihilation directly into SM particles via the s-channel Y resonance

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 > rsim 100 \,\, \mathrm {GeV}\) in the t-channel region. This is because the annihilation process \(\chi \chi \rightarrow 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_\mathrm{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 > rsim 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

$$\begin{aligned} \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} \end{aligned}$$
(21)

for intermediate values of \(g_\mathrm{SM}\), with larger values of \(\Omega _\chi h^2\) when \(g_\mathrm{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_\mathrm{DM}\lesssim \sqrt{4 \pi }\).

Figure 5 uses logarithmic scales to show the \((g_\mathrm{SM}, g_\mathrm{DM})\) plane in the vector-like model. Again, we distinguish immediately two distinct preferred regions, separated where \(g_\mathrm{SM}\lesssim 3 \times 10^{-4}\) and \(g_\mathrm{DM} > rsim 0.3\). As can be seen from the colour coding, the region at small \(g_\mathrm{SM}\) and large \(g_\mathrm{DM}\) corresponds to the triangular t-channel region in Fig. 4, whereas the larger values of \(g_\mathrm{SM}\)  and \(g_\mathrm{DM}\) outside this region appear in the s-channel funnel region.Footnote 8 This region is bounded at large values of the product \(g_\mathrm{SM}\) \(g_\mathrm{DM}\) by the upper limit on \(\sigma ^\mathrm{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_\mathrm{SM}\) \(g_\mathrm{DM}\) would be allowed for larger \(m_Y\). In the region \(0.07> g_\mathrm{DM}> 3 \times 10^{-4}\), \(\Omega h^2 \propto m_Y^2 / g_\mathrm{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

$$\begin{aligned} 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}} \end{aligned}$$
(22)

where \(\sigma _\mathrm{SI}^\mathrm{UL}(m_\chi )\) is the mass-dependent experimental upper bound on \(\sigma ^\mathrm{SI}_p\). We note also that the dijet constraint is the strongest at the largest value of \(g_\mathrm{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_\mathrm{DM}\lesssim 3 \times 10^{-4}\), rising at smaller \(g_\mathrm{SM}\), is also an artefact of the scanning limit \(m_Y > 100 \,\, \mathrm {GeV}\) combined with the relic density constraint, as indicated.

Fig. 5
figure 5

Preferred regions in the \((g_\mathrm{SM}, g_\mathrm{DM})\) plane (on logarithmic scales) in the vector-like model, where we have marginalized over the masses, \(m_\chi \) and \(m_Y\). We again use colour coding to illustrate the dominant mechanisms bringing the DM density into the allowed range: green for annihilation via t-channel exchange into pairs of mediator particles Y that subsequently decay into SM particles, and yellow for rapid annihilation directly into SM particles via the s-channel Y resonance. The diagonal dotted line indicates where \(g_\mathrm{DM}= g_\mathrm{SM}\), and the band where \(1/3< g_\mathrm{DM}/gsm < 3\) is bounded by dashed lines and shaded a darker yellow

We have indicated by a diagonal dotted line where \(g_\mathrm{DM}= g_\mathrm{SM}\). As discussed later, we might expect \(g_\mathrm{DM}\) and \(g_\mathrm{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_\mathrm{DM}/g_\mathrm{SM}< 3\), see Sect. 5.

We do not show the one-dimensional likelihood function for \(g_\mathrm{SM}\), which is quite featureless apart from a sharp rise for \(g_\mathrm{SM} > rsim 0.3\), nor that for \(g_\mathrm{DM}\), which is also featureless apart from a steep rise for \(g_\mathrm{DM}\lesssim 3 \times 10^{-4}\) that is an artefact of the limit on our scanning range for \(m_Y\).

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

Fig. 6
figure 6

The likelihood functions for \(g_\mathrm{SM}\) (left panel) and \(g_\mathrm{DM}\) (right panel) as functions of \(m_Y\) in the vector-like model. We again use colour coding to illustrate the dominant mechanisms bringing the DM density into the allowed range: green for annihilation via t-channel \(\chi \) exchange into pairs of mediator particles Y that subsequently decay into SM particles, and yellow for rapid annihilation directly into SM particles via the s-channel Y resonance

We see in the left panel of Fig. 6 that \(g_\mathrm{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 [55], where leptophobic models with universal quark couplings were analyzed,Footnote 9 and identifying \(g_\mathrm{SM}\) with the combination \(g Y^\prime _q\) of the parameters defined in those models, we see that the precision electroweak data do not constrain our model sample.Footnote 10 The two-branch structure of the t-channel region is a consequence of likelihood profiling. In fact, parameter space regions characterized by both s- and t-channel mechanisms overlap when projected in that specific region of the \(m_Y\)-\(g_{SM}\) plane and, since they have very similar \(\chi ^2\), sometimes one is selected over the other due to a statistical fluctuation.

In the right panel of Fig. 6, the lower bound on \(g_\mathrm{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_\mathrm{DM}\), which is given roughly by

$$\begin{aligned} \Big ( \frac{m_Y}{2 \,\mathrm TeV} \Big )^2 \Big ( \frac{g_\mathrm{DM}}{0.01} \Big )^{-2} \lesssim 1, \end{aligned}$$
(23)

as follows from Eq. (21) for \(2 m_{\chi } \sim m_Y\). The scattered yellow (green) islands in the green (yellow) region in the top-left corners of the plot is caused by the same mechanism that creates the two-branch structure in the left panel, i.e., overlapping regions of very similar \(\chi ^2\).

Fig. 7
figure 7

The likelihood functions for \(g_\mathrm{SM}\) (left panel) and \(g_\mathrm{DM}\) (right panel) as functions of \(m_\chi \) in the vector-like model. We again use colour coding to illustrate the dominant mechanisms bringing the DM density into the allowed range: green for annihilation via t-channel \(\chi \) exchange into pairs of mediator particles Y that subsequently decay into SM particles, and yellow for rapid annihilation directly into SM particles via the s-channel Y resonance

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

Figure 8 shows the likelihood function in the \((m_\chi , \sigma ^\mathrm{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 > rsim 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 ^\mathrm{SI}_p\) is allowed in both the t- and s-channel regions. The upper limits on \(\sigma ^\mathrm{SI}_p\) at various confidence levels are determined by the combined experimental likelihood for the LUX [76], PANDAX-II [77] and XENON1T [78] 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 ^\mathrm{SI}_p\) are allowed when \(m_\chi \sim m_Y/2\) and small values of \(g_\mathrm{DM}\) and/or \(g_\mathrm{SM}\) are favoured by the relic density constraint, and small values of \(\sigma ^\mathrm{SI}_p\) are allowed in the t-channel region because small values of \(g_\mathrm{SM}\) are allowed, as discussed previously. On the other hand, we see that \(\sigma ^\mathrm{SI}_p\) may well lie within the range to be probed by upcoming experiments such as LUX-ZEPLIN (LZ) [100] and XENONnT [101], though \(\sigma ^\mathrm{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 [82], updating [102, 103].

Fig. 8
figure 8

Contours of the likelihood in the \((m_\chi , \sigma ^\mathrm{SI}_p)\) plane for the vector-like model, showing the current upper limits from the LUX  [76], PANDAX-II [77] and XENON1T  [78] experiments (rescaled to account for the different local DM density that we assume), together with the neutrino ‘floor’ [102, 103] (shown as the dashed orange line), and the range of \(\sigma ^\mathrm{SI}_p\) that will be probed by the upcoming experiments LZ [100] and XENONnT [101]. We again use colour coding to illustrate the dominant mechanisms bringing the DM density into the allowed range: green for annihilation via t-channel \(\chi \) exchange into pairs of mediator particles Y that subsequently decay into SM particles, and yellow for rapid annihilation directly into SM particles via the s-channel Y resonance. We indicate the effective lower limit on \(m_\chi \) that is imposed by our sampling limit on \(m_Y\)

4.2 Axial-vector DM couplings

We now turn to the case of DM with axial-vector couplings. Figure 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/(2 m_\chi ) < 2\).

This broadening occurs because, whereas the direct detection constraint is very severe for the vector mediator, so that \(g_\mathrm{DM} \, g_\mathrm{SM}\) cannot be large and the parameters should be near the peak of the resonance where \(m_Y \sim 2 m_{\chi }\), in order for the DM particles to annihilate sufficiently (see Fig. 2 of [47] and the accompanying text), the strong \(\sigma ^\mathrm{SI}_p\) constraint is absent in the axial-vector case, so that off-resonance regions of parameter space with larger values of \(g_\mathrm{SM}\, g_\mathrm{DM}\) are allowed where annihilation is p-wave or \(m_q\) suppressed (see the discussion around Eq. (7) of [48]). This opens up more parameter space in the \((m_{Y}, m_{\chi })\) plane, with the deviation from \(m_Y \sim 2 m_{\chi }\) bounded only by the dijet constraint. We note that this constraint is weaker when \(m_Y > 2 m_\chi \) because the decay channel \(Y \rightarrow \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.

Fig. 9
figure 9

Preferred regions in the \((m_Y, m_\chi )\) plane in the model with axial-vector DM couplings. We delineate with red (blue) contours, respectively, 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, respectively. We use colour coding to illustrate the dominant mechanisms bringing the DM density into the allowed range: green for annihilation via t-channel exchange into pairs of mediator particles Y that subsequently decay into SM particles, and yellow for rapid annihilation directly into SM particles via the s-channel Y resonance

Also as in the vector-like case, we see again in the axial case in Fig. 9 the lower bound \(m_\chi > rsim 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_\mathrm{SM}, g_\mathrm{DM})\) plane in the axial-vector model, using logarithmic scales. The s-channel region extends to larger values of \(g_\mathrm{SM}\) and \(g_\mathrm{DM}\) than in the vector case, because of the suppression of axial-vector-mediated s-channel annihilation discussed above. Values of \(g_\mathrm{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_\mathrm{SM}\sim 10^{-2}\) and \(g_\mathrm{DM}\sim 10^{-1}\). The upper bound on \(g_\mathrm{SM}\) comes from the dijet constraint. As in the vector case, the dotted line is where \(g_\mathrm{SM}= g_\mathrm{DM}\), and the deeper shading indicates where \(1/3< g_\mathrm{SM}/g_\mathrm{DM}< 3\), favouring the large-\(g_\mathrm{SM}\) part of the s-channel annihilation region, see the discussion in Sect. 5. As in the vector case, there is a region at low \(g_\mathrm{DM}\), rising at small \(g_\mathrm{SM}\), whose exclusion by the relic density constraint is an artefact of the sampling restriction \(m_Y > 100 \,\, \mathrm {GeV}\).

Fig. 10
figure 10

Preferred regions in the \((g_\mathrm{SM}, g_\mathrm{DM})\) plane (on logarithmic scales) in the model with axial-vector DM couplings, again using colour coding to illustrate the dominant mechanisms bringing the DM density into the allowed range: green for annihilation via t-channel exchange into pairs of mediator particles Y that subsequently decay into SM particles, and yellow for rapid annihilation directly into SM particles via the s-channel Y resonance. The diagonal dotted line indicates where \(g_\mathrm{DM}= g_\mathrm{SM}\), and the band where \(1/3< g_\mathrm{DM}/g_\mathrm{SM}< 3\) is bounded by dashed lines and shaded a darker yellow

We display in Fig. 11 the \((m_Y, g_\mathrm{SM})\) and \((m_Y, g_\mathrm{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 > rsim 2.6 \,\, \mathrm {TeV}\) up to the sampling limit of \(5 \,\, \mathrm {TeV}\). The upper limit on \(g_\mathrm{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_\mathrm{SM}\). This is due to the p-wave nature of the s-channel annihilation via the axial mediator, which leads to the scaling behaviour , 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_\mathrm{DM}\) is given by the relic density constraint since. even exactly at \(m_Y = 2 m_\chi \), and DM is overproduced for sufficiently small \(g_\mathrm{DM}\).

Fig. 11
figure 11

The likelihood functions for \(g_\mathrm{SM}\) (left panel) and \(g_\mathrm{DM}\) (right panel) as functions of \(m_Y\) in the model with axial-vector couplings. We again use colour coding to illustrate the dominant mechanisms bringing the DM density into the allowed range: green for annihilation via t-channel \(\chi \) exchange into pairs of mediator particles Y that subsequently decay into SM particles, and yellow for rapid annihilation directly into SM particles via the s-channel Y resonance

Figure 12 shows the \((m_\chi , g_\mathrm{SM})\) and \((m_\chi , g_\mathrm{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_\mathrm{SM}\) where the t-channel mechanism dominates and at larger \(g_\mathrm{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_\mathrm{DM}\) comes from the relic density constraint as can be seen in Eqs. (12), (13) and (14).

Fig. 12
figure 12

The likelihood functions for \(g_\mathrm{SM}\) (left panel) and \(g_\mathrm{DM}\) (right panel) as functions of \(m_\chi \) in the model with axial-vector couplings. We again use colour coding to illustrate the dominant mechanisms bringing the DM density into the allowed range: green for annihilation via t-channel \(\chi \) exchange into pairs of mediator particles Y that subsequently decay into SM particles, and yellow for rapid annihilation directly into SM particles via the s-channel Y resonance

The one-dimensional \(\Delta \chi ^2\) function for \(g_\mathrm{DM}\) in the axial case rises sharply below \(\sim 4 \times 10^{-3}\), and that for \(g_\mathrm{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 ^\mathrm{SD}_p\), inferred from the PICO-60 search with a C\(_3\)F\(_8\) target) [79] (left panel) and to that on a neutron (\(\sigma ^\mathrm{SD}_n\), inferred from a search with the XENON1T detector) [80] (right panel), again accounting for the different local DM density that we assume. Since we consider here leptophobic mediators, the constraints on \(\sigma ^\mathrm{SD}_p\) provided by Super-Kamiokande [104] and IceCube [105] 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 ^\mathrm{SD}_p\). We also show in the left panel the estimated neutrino ‘floor’ applicable to experiments using a C\(_3\)F\(_8\) target (shaded blue), adapted from [106] using a similar factor as in [82]. We see that \(\sigma ^\mathrm{SD}_p\) approaches the PICO-60 limit most closely for a small range \(100~\mathrm{GeV} \lesssim m_\chi \lesssim 200\) GeV, and that \(\sigma ^\mathrm{SD}_p\) may be accessible to the LZ experiment [100] or the PICO-500 experiment [106] 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 ^\mathrm{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 [106] (shaded blue). Encouragingly, we note that there are significant regions of the axial-vector parameter space where both \(\sigma ^\mathrm{SD}_p\) and \(\sigma ^\mathrm{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 ^\mathrm{SD}_n\), since that currently provides the stronger constraint, the longer-term advantage may lie with searches for \(\sigma ^\mathrm{SD}_p\), since the expected ‘floor’ is lower in that case.

Fig. 13
figure 13

Predictions for \(\sigma ^\mathrm{SD}_p\) and \(\sigma ^\mathrm{SD}_n\) in the DMSM with axial couplings. Left panel: Contours of the likelihood function in the \((m_\chi , \sigma ^\mathrm{SD}_p)\) plane for the axial-vector model, showing the rescaled upper limit from the PICO-60 experiment [79] and the prospective sensitivity of the PICO-500 [106] and LZ [100] experiments to \(\sigma ^\mathrm{SD}_p\), as well as the neutrino ‘floor’ applicable to an experiment such as PICO-500 that uses C\(_3\)F\(_8\) (shaded blue). Right panel: Contours of the likelihood function in the \((m_\chi , \sigma ^\mathrm{SD}_n)\) plane for the axial-vector model, showing the rescaled upper limit from the XENON1T experiment [80] and the prospective sensitivity of the LZ [100] experiment to \(\sigma ^\mathrm{SD}_n\), as well as the neutrino ‘floor’ applicable to an experiment that uses Xenon (shaded blue)

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 [9,10,11,12,13], 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 \(\mathcal{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 [54], 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_\mathrm{DM}/g_\mathrm{SM}= \mathcal{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 \( > rsim \mathcal{O}(0.1)\), favouring parts of the funnel regions away from the regions where t-channel annihilation is dominant.

We have included in Figs. 5 and 10 diagonal dotted lines where \(g_\mathrm{DM}= g_\mathrm{SM}\) and shaded the strips where \(1/3< g_\mathrm{SM}/g_\mathrm{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_\mathrm{DM}\gg g_\mathrm{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_\mathrm{SM}/g_\mathrm{DM}< 3\), as we discuss later. We now discuss in more detail the implications of the constraint \(1/3< g_\mathrm{SM}/g_\mathrm{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_\mathrm{SM}/g_\mathrm{DM}< 3\). We see that both \(g_\mathrm{SM}\) and \(g_\mathrm{DM}\) are \(> 10^{-3}\) when this selection is imposed, and that values of \(g_\mathrm{SM}\) and \(g_\mathrm{DM}\lesssim 0.3\) lie within the favoured range, consistent with unification scenarios. The ranges \(m_Y > rsim 2 \,\, \mathrm {TeV}\) and \(m_\chi > rsim 1 \,\, \mathrm {TeV}\) are compatible with \(g_\mathrm{SM}, g_\mathrm{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_\mathrm{SM}\) and \(g_\mathrm{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_\mathrm{SM}\) for \(m_Y > rsim 1.8\) TeV (\(m_\chi > rsim 0.9\) TeV). These features are also visible in the region of Fig. 5 with the darker yellow band.

Fig. 14
figure 14

The likelihood functions for \(g_\mathrm{SM}\) (left panel) and \(g_\mathrm{DM}\) (right panel) as functions of \(m_\chi \) in the model with vector couplings after the selection \(1/3< g_\mathrm{SM}/g_\mathrm{DM}< 3\). We again use colour coding to illustrate the dominant mechanisms bringing the DM density into the allowed range: green for annihilation via t-channel \(\chi \) exchange into pairs of mediator particles Y that subsequently decay into SM particles, and yellow for rapid annihilation directly into SM particles via the s-channel Y resonance

Figure 15 shows the corresponding parameter planes for the axial-coupling case with the selection \(1/3< g_\mathrm{SM}/g_\mathrm{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_\mathrm{SM} > rsim 0.1\). We see from the right panel of Fig. 11 that \(g_\mathrm{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_\mathrm{SM}/g_\mathrm{DM}< 3\) condition then implies that \(g_\mathrm{SM} > rsim 0.1\) (0.3) for \(m_Y = 500\) GeV (2 TeV). However, such a large value of \(g_\mathrm{SM}\) is not compatible with the dijet constraint in this range of \(m_Y\). Therefore, no region with \(m_Y > rsim 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_\mathrm{SM} > rsim 10^{-3}\) and \(g_\mathrm{DM} > rsim 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_\mathrm{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_\mathrm{SM}> 0.1\).

In the lower panels of Fig. 15 we see the corresponding \((m_\chi , g_\mathrm{SM})\) and \((m_\chi , g_\mathrm{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_\mathrm{SM}\lesssim 0.3\), and in the right plot its range is \(0.05 \lesssim g_\mathrm{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 > rsim 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_\mathrm{DM}, g_\mathrm{SM}> 0.1\) are possible for all the range of \(m_\chi \).

Fig. 15
figure 15

The likelihood functions for \(g_\mathrm{SM}\) (left panel) and \(g_\mathrm{DM}\) (right panel) as functions of \(m_\chi \) in the model with axial-vector couplings after the selection \(1/3< g_\mathrm{SM}/g_\mathrm{DM}< 3\). We again use colour coding to illustrate the dominant mechanisms bringing the DM density into the allowed range: green for annihilation via t-channel \(\chi \) exchange into pairs of mediator particles Y that subsequently decay into SM particles, and yellow for rapid annihilation directly into SM particles via the s-channel Y resonance

In the vector-like case, the main changes in the one-dimensional \(\Delta \chi ^2\) functions for \(g_\mathrm{DM}\) and \(g_\mathrm{SM}\) due to the selection \(1/3< g_\mathrm{SM}/g_\mathrm{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_\mathrm{DM}\lesssim 0.5\) and \(g_\mathrm{SM} > rsim 10^{-3}\), while all values of \(m_\chi > rsim 50 \,\, \mathrm {GeV}\) and \(m_Y > rsim 100 \,\, \mathrm {GeV}\) still have very low \(\Delta \chi ^2\).

Finally, we show in Fig. 16 the preferred ranges of \((m_\chi , \sigma ^\mathrm{SI}_p)\) in the vector-like model (left panel) and \((m_\chi , \sigma ^\mathrm{SD}_n)\) in the axial model (right panel) after the selection \(1/3< g_\mathrm{SM}/g_\mathrm{DM}< 3\). We see that there are good prospects for detecting spin-independent DM scattering in the vector-like case, since \(\sigma ^\mathrm{SI}_p\) lies partially above the neutrino ‘floor’ within the range of \(m_\chi \) sampled, though smaller values of \(\sigma ^\mathrm{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_\mathrm{SM}/g_\mathrm{DM}< 3\), we see in the right panel of Fig. 16 that in the t-channel exchange region the predicted values of \(\sigma ^\mathrm{SD}_n\) are relatively close to the XENON1T limit [80], and hence may offer prospects for detection with the planned LZ experiment [100] within the sampled range of DM masses. However, values of \(\sigma ^\mathrm{SD}_n\) may be considerably lower in the s-channel exchange region.Footnote 11 The disfavored region for \(10^{-45}~\mathrm {cm}^2 \lesssim \) \(\sigma ^\mathrm{SD}_n\) \(\lesssim 10^{-42}~\mathrm {cm}^2\) and \(m_\chi > rsim 300\) GeV is caused by the constraint coming from dijet searches.

The reciprocal pattern of the two mechanisms 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_\mathrm{SM}\) and \(g_\mathrm{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_\mathrm{SM}\) and \(g_\mathrm{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.

Fig. 16
figure 16

Predictions for \(\sigma ^\mathrm{SI}_p\) and \(\sigma ^\mathrm{SD}_n\) after the selection \(1/3< g_\mathrm{SM}/g_\mathrm{DM}< 3\). Left panel: Contours of the likelihood function in the \((m_\chi , \sigma ^\mathrm{SI}_p)\) plane for the vector-like model, showing the current upper limits from the LUX [76], PANDAX-II [77] and XENON1T [78] experiments together with the neutrino ‘floor’ [102, 103] (shown as the dashed orange line), and the range of \(\sigma ^\mathrm{SI}_p\) that will be probed by the upcoming experiments LZ [100] and XENONnT [101]. Right panel: Contours of the likelihood function in the \((m_\chi , \sigma ^\mathrm{SD}_n)\) plane for the axial-like model, showing the upper limit from the XENON1T experiment [80] and the prospective sensitivity of the LZ experiment that also uses a Xenon target [100]

It is interesting that, with the selection \(1/3< g_\mathrm{SM}/g_\mathrm{DM}< 3\) motivated by the prospect of UV completion, searches for \(\sigma ^\mathrm{SI}_p\) and \(\sigma ^\mathrm{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_\mathrm{DM}\) of the mediator to the dark matter particle, and the coupling \(g_\mathrm{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 \le 5 \,\, \mathrm {TeV}\) and dark matter particle masses \(m_\chi \le 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_\mathrm{DM}/g_\mathrm{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_\mathrm{SM}\) and relatively large values of \(g_\mathrm{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_\mathrm{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_\mathrm{DM} > rsim 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 ^\mathrm{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 ^\mathrm{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_\mathrm{SM}\) is quite featureless, as is that for \(g_\mathrm{DM} > rsim 10^{-2}\). We find that the spin-dependent dark matter scattering cross sections \(\sigma ^\mathrm{SD}_p\) and \(\sigma ^\mathrm{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_\mathrm{DM}/g_\mathrm{SM}< 3\), as might be suggested in some ultraviolet completions of the dark matter simplified models considered here. In this case, values of \(g_\mathrm{SM}, g_\mathrm{DM}\sim 0.1\) are possible, 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_\mathrm{DM}/g_\mathrm{SM}\) disfavours the t-channel mechanism. However, in the axial case the t-channel mechanism is still possible, yields values of \(\sigma ^\mathrm{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 ^\mathrm{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 [9,10,11,12,13] using the MasterCode framework. Our analysis also provides a more detailed study of the vector and axial leptophobic models considered in [55]. 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.