Global analysis of dark matter simplified models with leptophobic spin-one mediators using MasterCode

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 ≲10-3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lesssim \, 10^{-3}$$\end{document} 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.


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 weaklyinteracting 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 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].
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 χ . The Lagrangian for the spin-one mediator Y interactions with the DM particle χ is and that for its interactions with quarks is For simplicity, following the DMWG [12] we assume that the interactions between the Y -boson and quarks are generationindependent, and we also assume that they are the same for charge-2/3 and -1/3 quarks (g V,A q := g V,A q i ∀ i). 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: and another in which they are purely axial: 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 SM , g 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 χ and m Y , and the mediator couplings to the DM and SM particles, g DM and g 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 χ h 2 requires values of g DM and g SM that depend on m χ 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 χ 50 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.

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]: We assume here that the dominant source of this cold DM is the WIMP χ , so that χ h 2 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, χ h 2 ∝ 1/ σ v . If the annihilation is dominated by the s-channel process χχ → Y ( * ) → Standard Model particles, the leading expression for the cross-section takes the form [72,73] σ v s c Y 32π where c Y is an 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: where X = χ (DM), q (SM), and where 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 SM 1, and one finds that in the limit g SM → 0 so the scaling of the relic density is expected to be near the Y peak. However, we recall that χχ annihilation takes place in a thermal bath, with collisions taking place at a range of centre-of-mass energies ≥ 2m χ . For this reason, annihilations at an effective centre-of-mass energy within O( Y ) of m Y do not dominate when Y is very small. We find numerically that (11) captures the g DM dependence of χ h 2 , and that it depends only weakly on g SM over an intermediate range of g SM , but that χ h 2 increases at both large and small g SM : for g SM g DM and for g SM g DM . The dependence on g DM in this case depends on the relative sizes of the width and mass difference between m Y and 2m χ [73].
However, if the width is very small (for example if g SM is very small and m Y < 2m χ ), then the denominator in (6) is dominated by the mass difference and χ h 2 is inversely proportional to g 2 DM , as seen in Eq. (9) for the axial-vector case.
Another important process is the t-channel process χχ → Y Y for m χ > m Y . The leading contribution to the t-channel cross-section is given by [72] and is kinematically forbidden when m χ < m Y . We use MicrOMEGAs 5.0.2 [64] to evaluate χ h 2 numerically in each model parameter set, including a ±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 χ h 2 (colours and solid blue line) with the χ 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 V SM = 0.25, g V DM = 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 χ h 2 > 0.12, and that the DM particle is underdense in the blue region at large m χ .

Spin-independent DM scattering
We evaluate the constraint on the spin-independent DM scattering cross-section σ 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 σ SI p using MicrOMEGAs 5.0.2. This constraint is important for the vector mediator case, where the spinindependent cross-section has the following scaling [11] σ SI p ∝ g 2 DM g 2 SM /m 4 Y , with the following approximate numerical value: in the s-channel region where m Y 2m χ . The σ SI p constraint provides important limits on the vector mediator mass m Y and the product g DM g SM .
We allow for an overall uncertainty of 10% in the spinindependent cross section.

Spin-dependent DM scattering
We also evaluate the constraint imposed by the upper limit on the spin-dependent DM-proton scattering cross-section σ SD p from the PICO-60 experiment [79], and that on the spindependent DM-neutron scattering cross-section σ SD n from the XENON1T experiment [80]. We evaluate theoretical predictions for σ SD p and σ 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.

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 χ 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 σ SI p and σ SD p , σ SD n , we assume an overall uncertainty of 30% [82,85], in the local density.

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 UL i (m), where i labels the vector and axial-vector We model the likelihood function from this search using the procedure outlined in Fastlim [86]: where σ 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.

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 itsqq coupling, g q . In recasting this constraint for the pp → Y → qq process in our DM models, we demand 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 . 2 We note that σ Y →qq and σ (g * q ) can be written as (using x ≡ (Y → x x)) 2 We note that σ (g * q ) is the cross-section limit assuming the finite width given by g * q . Due to the presence of the Y → χχ mode, Y and Z may be different. However, when this limit is important, the Y → qq branching ratio must be large, which justifies our approximation that with the same constant c at m Y = m Z , and q = g SM g * q 2 * q . Using these relations, the equality in Eq. (17) occurs when g SM takes the value which serves as the upper bound on g SM for given m Y , q and χ . We model the corresponding likelihood function using where q and χ 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. 3 The coloured lines correspond to the various experimental analyses, which can directly compared to the MasterCode χ 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 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 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 → tt mode, since the sensitivity is weaker than in the dijet mode, due to the large tt 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 modeldependent 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  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 data 5 and indirect searches for astrophysical DM annihilations should be considered, 6 but include all the Table 1 The ranges of the DMSM parameters sampled, together with the numbers of segments into which they were divided during the sampling Total # of segments 320 masses for which LHC searches are sensitive. The couplings g SM and g DM were restricted to perturbative ranges < √ 4π . 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 χ /m Y deviates from unity by < 10 −3 , so as to sample adequately annihilations near the Y peak when Y m Y . We used similar sampling procedure for both the vector and axial-vector DMSMs, generating ∼ 100 million parameter sets in each case. The parameter regions with χ 2 < 2.30 (5.99), which are favoured at the 68% (95%) CL and regarded as proxies for 1-(2-)σ 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 χ 2 in the displayed 2dimensional projection of the four-dimensional parameter space using colour coding:

Vector DM couplings
• Green: annihilation via t-channel χ exchange into pairs of mediator particles Y that subsequently decay into SM particles, in the region where m χ ≥ m Y ; • Yellow: rapid annihilation directly into SM particles via the s-channel Y resonance, in the region where 0.9 < m Y /(2m χ ) < 1.1.
We see clearly two separated favoured regions, one with m χ m Y /2 where rapid annihilation via the Y funnel dom-Footnote 6 continued neutrinos are not competitive with direct searches for spin-dependent DM scattering, as discussed below.

Fig. 4
Preferred regions in the (m Y , m χ ) plane in the vector-like model. We delineate with red (blue) contours, respectively, the parameter regions with χ 2 < 2.30(5.99), which are favoured at the 68% (95%) CL and regarded as proxies for 1-(2-)σ regions, respectively. We also use colour coding to illustrate the dominant mechanisms bringing the DM density into the allowed range: green for annihilation via tchannel 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 inates, and another with m χ > 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 phasespace limit in the t-channel case. 7 We see that the t-channel region is located at smaller m Y than the s-channel region for any value of m χ , subject to the kinematic constraint m Y < m χ . In this region the annihilation cross-section of this channel, χχ → Y Y , is proportional to g 4 DM and independent of g SM . Therefore, the relic abundance can be brought to the observed value by tuning g DM , unless m χ 100 GeV, whilst the experimental constraints from the LHC and the direct DM detection can be avoided by taking g SM small enough.
As can be seen in Fig. 4, the lower limit m Y > 100 GeV in our sample enforces a corresponding lower limit m χ 100 GeV in the t-channel region. This is because the annihilation process χχ → Y Y is kinematically blocked for m χ < m Y , so the cross-section for χχ annihilation becomes very small, resulting in DM overdensity. On the other hand, in the s-channel region the relic density constraint requires g DM to be very small for m Y ∼ 100 GeV, in which case We find that the χ 2 function is negligible for all values of m Y above the 100 GeV cut, and we also find that the χ 2 function is also very small and essentially featureless for m χ > 50 GeV. Moreover, we find no upper limits on m χ,Y within the ranges sampled. Our numerical results indicate that for intermediate values of g SM , with larger values of χ h 2 when g 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 ∼ O(700) TeV for g DM √ 4π . Figure 5 uses logarithmic scales to show the (g SM , g DM ) plane in the vector-like model. Again, we distinguish immediately two distinct preferred regions, separated where g SM 3 × 10 −4 and g DM 0.3. As can be seen from the colour coding, the region at small g SM and large g DM corresponds to the triangular t-channel region in Fig. 4, whereas the larger values of g SM and g DM outside this region appear in the s-channel funnel region. 8 where σ UL SI (m χ ) is the mass-dependent experimental upper bound on σ SI p . We note also that the dijet constraint is the strongest at the largest value of g SM ∼ 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 DM 3×10 −4 , rising at smaller g SM , is also an artefact of the scanning limit m Y > 100 GeV combined with the relic density constraint, as indicated.
We have indicated by a diagonal dotted line where g DM = g SM . As discussed later, we might expect g DM and g 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 DM /g SM < 3, see Sect. 5.
We do not show the one-dimensional likelihood function for g SM , which is quite featureless apart from a sharp rise for 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 DM = g SM , and the band where 1/3 < g DM /gsm < 3 is bounded by dashed lines and shaded a darker yellow g SM 0.3, nor that for g DM , which is also featureless apart from a steep rise for g DM 3 × 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 SM , g DM . We see in the left panel that within the t-channel region there is a strong upper limit on g SM 10 −3 , which is enforced by the combination of the upper limit on σ SI p , which constrains the product g DM g SM , and the lower limit on g DM visible in the right panel, which is another result of our limited scan to m Y > 100 GeV. Since σ SI p scales as g 2 DM g 2 SM /m 4 Y , the limit is stronger for smaller m Y , and g SM must be smaller than 10 −4 for m Y 100 GeV. The upper limit on m Y in the t-channel region visible for small g SM is an artefact of the limited scanning range m χ < 2.5 TeV. The upper limit on g SM in the s-channel region comes mainly from the LHC dijet constraint, with the σ SI p constraint also playing a role when m Y 500 GeV.
We see in the left panel of Fig. 6 that g SM 10 −2 for m Y ∼ 100 GeV and 0.1 for m Y > 1 TeV. Comparing with Fig. 11 of [55], where leptophobic models with universal quark couplings were analyzed, 9 and identifying g SM with the combination gY q of the parameters defined in those models, we see that the precision electroweak data do not constrain our model sample. 10 The two-branch structure of Fig. 6 The likelihood functions for g SM (left panel) and g 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 χ 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 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 χ 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 DM in the tchannel 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 DM , which is given roughly by as follows from Eq. (21) for 2m χ ∼ m Y . The scattered yellow (green) islands in the green (yellow) region in the topleft 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 χ 2 . Figure 7 displays the (m χ , g SM ) plane in the left panel and the (m χ , g DM ) plane in the right panel. The upper bounds on g SM and the lower bounds on g DM are the same as in Fig. 6, all increasing with m χ . Figure 8 shows the likelihood function in the (m χ , σ 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 χ 50 GeV are allowed, for the reason mentioned previously, namely the interplay of the cut m Y > 100 GeV (as indicated) and the relic density constraint. Above this value, a large range of values of σ SI p is allowed in both the t-and s-channel regions. The upper limits on σ SI p at various confidence levels are deter-Footnote 10 continued interference, would contribute to cross sections and rates only at the two-loop level, and hence can be neglected. mined 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 σ SI p are allowed when m χ ∼ m Y /2 and small values of g DM and/or g SM are favoured by the relic density constraint, and small values of σ SI p are allowed in the t-channel region because small values of g SM are allowed, as discussed previously. On the other hand, we see that σ SI p may well lie within the range to be probed by upcoming experiments such as LUX-ZEPLIN (LZ) [100] and XENONnT [101], though σ 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].

Axial-vector DM couplings
We now turn to the case of DM with axial-vector couplings. Figure 9 displays the (m Y , m χ ) 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-σ 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 χ ) < 2.
This broadening occurs because, whereas the direct detection constraint is very severe for the vector mediator, so that g DM g SM cannot be large and the parameters should be near the peak of the resonance where m Y ∼ 2m χ , in order for the DM particles to annihilate sufficiently (see Fig. 2 of [47] and the accompanying text), the strong σ SI p constraint is absent in  [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 σ 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 χ 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 χ that is imposed by our sampling limit on m Y the axial-vector case, so that off-resonance regions of parameter space with larger values of g SM g 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 χ ) plane, with the deviation from m Y ∼ 2m χ bounded only by the dijet constraint. We note that this con- Fig. 9 Preferred regions in the (m Y , m χ ) plane in the model with axial-vector DM couplings. We delineate with red (blue) contours, respectively, the parameter regions with χ 2 < 2.30(5.99), which are favoured at the 68% (95%) CL and regarded as proxies for 1-(2-)σ 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 straint is weaker when m Y > 2m χ because the decay channel Y → χχ 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 χ > 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 χ 50 GeV due to the interplay of the sampling limit m Y > 100 GeV and the relic density constraint. Above this value of m χ , the one-dimensional Fig. 10 Preferred regions in the (g SM , g 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 DM = g SM , and the band where 1/3 < g DM /g SM < 3 is bounded by dashed lines and shaded a darker yellow χ 2 likelihood function for m χ is featureless, as is that for m Y > 100 GeV.
In Fig. 10 we show the (g SM , g DM ) plane in the axialvector model, using logarithmic scales. The s-channel region extends to larger values of g SM and g DM than in the vector case, because of the suppression of axial-vector-mediated s-channel annihilation discussed above. Values of g DM as large as the sampling limit √ 4π are allowed. We see a separation between the regions of this plane where the s-and t-channel mechanisms are dominant when g SM ∼ 10 −2 and g DM ∼ 10 −1 . The upper bound on g SM comes from the dijet constraint. As in the vector case, the dotted line is where g SM = g DM , and the deeper shading indicates where 1/3 < g SM /g DM < 3, favouring the large-g 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 DM , rising at small g SM , whose exclusion by the relic density constraint is an artefact of the sampling restriction m Y > 100 GeV.
We display in Fig. 11 the (m Y , g SM ) and (m Y , g 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 2.6 TeV. The limit visible in Fig. 9, which is due to the sampling limit m χ < 2.5 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 2.6 TeV up to the sampling limit of 5 TeV. The upper limit on g SM for m Y < 5 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 SM . This is due to the p-wave nature of the s-channel annihilation via the axial mediator, which leads to the scaling behaviour h 2 ∝ ∼ m 2 Y /(g 2 SM g 2 DM ), as discussed around Eq. (9).
We see again in the right panel that the t-channel mechanism is relevant for m Y 3 TeV, whereas the s-channel mechanism is dominant for larger m Y . The lower limit on g DM is given by the relic density constraint since. even exactly at m Y = 2m χ , h 2 ∝ ∼ m 2 Y /(g 2 SM g 2 DM ) and DM is overproduced for sufficiently small g DM . Figure 12 shows the (m χ , g SM ) and (m χ , g DM ) planes (left and right panels, respectively) in the axial-vector model. We see again in the left panel that values of m χ up to the sampling limit of 2.5 TeV are allowed both at smaller g SM where the t-channel mechanism dominates and at larger g SM where the s-channel mechanism dominates. There is a region at low m χ where the s-channel mechanism dominates. The lower limit on g DM comes from the relic density constraint as can be seen in Eqs. (12), (13) and (14).
The one-dimensional χ 2 function for g DM in the axial case rises sharply below ∼ 4 × 10 −3 , and that for g SM rises above ∼ 0.3. The χ 2 functions for m χ 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 (σ SD p , inferred from the PICO-60 search with a C 3 F 8 target) [79] (left panel) and to that on a neutron (σ 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 σ SD p provided by Super-Kamiokande [104] and IceCube [105] limits on the annihilations into τ + τ − of DM particles trapped in the Sun are not relevant, and the limits on hadronic annihilations are not competitive with the direct constraints on σ 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 σ SD p approaches the PICO-60 limit most closely for a small range 100 GeV m χ 200 GeV, and that σ SD p may be accessible to the LZ experiment [100] or the PICO-500 experiment [106] for m χ 1 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 σ 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 σ SD p and σ SD n may be detectable above the corresponding Eur. Phys. J. C (2019) 79:895 Fig. 11 The likelihood functions for g SM (left panel) and g 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 χ 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 Fig. 12 The likelihood functions for g SM (left panel) and g DM (right panel) as functions of m χ 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 χ 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 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 σ SD n , since that currently provides the stronger constraint, the longer-term advantage may lie with searches for σ SD p , since the expected 'floor' is lower in that case.

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 comple-tion, 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 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)×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 axiallike couplings to quarks, and DMSMs with purely vectoror 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 DM /g SM = 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 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 DM = g SM and shaded the strips where 1/3 < g SM /g 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 DM g 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 SM /g DM < 3, as we discuss later. We now discuss in more detail the implications of the constraint 1/3 < g SM /g 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 SM /g DM < 3. We see that both g SM and g DM are > 10 −3 when this selection is imposed, and that values of g SM and g DM 0.3 lie within the favoured range, consistent with unification scenarios. The ranges m Y 2 TeV and m χ 1 TeV are compatible with g SM , g 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 SM and g 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 SM for m Y 1.8 TeV (m χ 0.9 TeV). These features are also visible in the region of Fig. 5 with the darker yellow band. Figure 15 shows the corresponding parameter planes for the axial-coupling case with the selection 1/3 < g SM /g DM < 3. In the upper panels, we see that the tchannel region is confined to a small region where m Y 500 GeV and g SM 0.1. We see from the right panel of Fig. 11 that g DM must be larger than ∼ 0.3 (0.9) at m Y = 500 GeV (2 TeV) to account for the DM density. The 1/3 < g SM /g DM < 3 condition then implies that g SM 0.1 (0.3) for m Y = 500 GeV (2 TeV). However, such a large value of g SM is not compatible with the dijet constraint in this range of m Y . Therefore, no region with m Y 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 SM 10 −3 and g DM 3 × 10 −3 (see the darker yellow band in Fig. 10), are needed to compensate this suppression. However, there is a narrow band of g 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 SM > 0.1. In the lower panels of Fig. 15 we see the corresponding (m χ , g SM ) and (m χ , g DM ) planes. We see in the left plot that in this projection the t-channel region appears in the coupling range 0.03 g SM 0.3, and in the right plot its range is 0.05 g DM 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 χ direction, and any value of m χ 300 GeV is in principle compatible with this mechanism. In the s-channel region all values of m χ between 50 GeV and 2.5 TeV are allowed. We note also that g DM , g SM > 0.1 are possible for all the range of m χ .
In the vector-like case, the main changes in the onedimensional χ 2 functions for g DM and g SM due to the selection 1/3 < g SM /g DM < 3 are reductions in their upper limits to ∼ 0.3, and there are negligible changes in the one-dimensional χ 2 functions for m Y and m χ . In the axial case the main changes when the selection is made are that g DM 0.5 and g SM 10 −3 , while all values of m χ 50 GeV and m Y 100 GeV still have very low χ 2 . Finally, we show in Fig. 16 the preferred ranges of (m χ , σ SI p ) in the vector-like model (left panel) and (m χ , σ SD n ) in the axial model (right panel) after the selection 1/3 < g SM /g DM < 3. We see that there are good prospects for detecting spin-independent DM scattering in the vector-like case, since σ SI p lies partially above the neutrino 'floor' within the range of m χ sampled, though smaller values of σ SI p could also occur for any value of m χ > 50 GeV. In the case of spin-dependent scattering after the selection 1/3 < g SM /g DM < 3, we see in the right panel of Fig. 16 that in the t-channel exchange region the predicted values of σ 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 σ SD n may be considerably lower in the s-channel Fig. 15 The likelihood functions for g SM (left panel) and g DM (right panel) as functions of m χ in the model with axial-vector couplings after the selection 1/3 < g SM /g 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 χ 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 exchange region. 11 The disfavored region for 10 −45 cm 2 σ SD n 10 −42 cm 2 and m χ 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 SM and g 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 SM and g DM than in the s-channel region for a large range of values of m χ . The intermediate 'hole' is enforced by the LHC dijet constraint, in particular.
It is interesting that, with the selection 1/3 < g SM /g DM < 3 motivated by the prospect of UV completion, searches for 11 Note that we have again not included in the right panel the upper limits and prospective detection 'floors' presented by the Super-Kamiokande [104] and IceCube [105] Collaborations, which are based on specific assumptions about the DM annihilation channels that are not applicable in the leptophobic DMSMs studied here. σ SI p and σ SD n are able to probe both the s-and t-channel regions of the DMSM parameter spaces

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 χ . Each of these models is characterized by four free parameters: m Y and m χ , the coupling g DM of the mediator to the dark matter particle, and the coupling g 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  [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 σ 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 χ , σ 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] upper limits on spin-independent and -dependent scattering on nuclei. We have scanned mediator masses m Y ≤ 5 TeV and dark matter particle masses m χ ≤ 2.5 TeV, delineating the regions of the model parameters with χ 2 < 2.30(5.99), which are favoured at the 68% (95%) CL and regarded as proxies for 1-(2-)σ 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 χ 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 DM /g 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 SM and relatively large values of g DM . The one-dimensional likelihood functions for both m χ and m Y are quite small and flat above thresholds ∼ 50 GeV and ∼ 100 GeV. Thus the LHC still has interesting prospects for discovering DM and mediator particles in these simplified models. Any value of g SM between ∼ 10 −6 and the dijet limit of ∼ 0.3 is possible without any χ 2 penalty, as is the case for g DM 10 −3 , where the lower limit is due to the upper limit on the sampling range for m Y . The spinindependent dark matter scattering cross section σ 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 σ 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 χ and m Y are again featureless above thresholds ∼ 50 GeV and 100 GeV, respectively, and that for g SM is quite featureless, as is that for g DM 10 −2 . We find that the spindependent dark matter scattering cross sections σ SD p and σ 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 DM /g SM < 3, as might be suggested in some ultraviolet completions of the dark matter simplified models considered here. In this case, values of g SM , g DM ∼ 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 DM /g SM disfavours the t-channel mechanism. However, in the axial case the t-channel mechanism is still possible, yields values of σ 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 σ SD n are possible if the t-channel mechanism dominates.
In this paper, we have performed a first global studyincluding 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.