Constraints on Z' models from LHC dijet searches and implications for dark matter

We analyse a combination of ATLAS and CMS searches for dijet resonances at run I and run II, presenting the results in a way that can be easily applied to a generic Z' model. As an illustrative example, we consider a simple model of a Z' coupling to quarks and dark matter. We first study a benchmark case with fixed couplings and then focus on the assumption that the Z' is responsible for setting the dark matter relic abundance. Dijet constraints place significant bounds on this scenario, allowing us to narrow down the allowed range of dark matter masses for given Z' mass and width.


Introduction
Resonant structures in the invariant mass distribution of dijet events are amongst the most generic signatures for "exotic" new physics at the LHC, since any new heavy particle produced in the s-channel at hadron colliders can decay back into a pair of jets. Searches for dijet resonances are therefore a high priority at both ATLAS and CMS and have been among the first searches carried out at a centre-of-mass energy of 13 TeV [1,2]. These searches are complemented by earlier searches at 8 TeV [3,4], as well as a dedicated search for dijet resonances with an invariant mass below 1 TeV at CMS based on a novel data scouting technique [5]. Among the many models probed by such searches are Randall-Sundrum (RS) gravitons [6], excited quarks [7,8] and models with a leptophobic Z [9][10][11] (see [12] for the implications of these searches on more generic Z models).
The latter model of a massive spin-1 boson that couples predominantly to quarks has received a significant amount of interest in the context of dark matter (DM) production at hadron colliders. The reason is that a leptophobic Z can have large couplings to the DM particle and thereby mediate the interactions that keep DM in thermal equilibrium in the early Universe. We can then hope to experimentally probe these interactions with a range of different DM search experiments. Two sensitive probes of such scenarios are direct detection experiments searching for evidence of DM-nucleus scattering and searches for missing energy at the LHC [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]. In fact, DM models with a leptophobic Z have served to inspire a class of so-called simplified DM models, which are now commonly used to optimise LHC searches for DM [29,30].
However, as emphasised in a number of previous works [13,25,31,32], DM models with a massive spin-1 mediator cannot only be probed by conventional DM searches, but also by direct searches for the mediator, in particular searches for dijet resonances. Due to the presence of invisible decays of the mediator, the width of resonance may be broadened, making it harder to distinguish the signal from the smoothly-falling QCD backgrounds. The purpose of the present work is to derive combined limits on such models by re-analysing all available LHC searches for dijet resonances.
Previous combinations of dijet constraints have either focussed on narrow resonances [13], resonances that decay exclusively into quarks [10] or on individual DM models with specific choices of couplings [25,31]. The present work aims to take a largely model-independent approach, so that our results can be applied to a range of different Z models. For this purpose, we take the width of the Z as a free parameter, which we allow to be as large as Γ/m Z = 0.3. The resulting bound on the Z -quark coupling can then be applied to models where the Z decays exclusively to quarks (for example if couplings to DM are negligible or absent, or if decays of the Z into DM available are kinematically forbidden), to models where the width of the Z is given exclusively in terms of its couplings to quarks and DM (as is the case in simplified DM models) and to models where additional unobserved decay channels may be present that further broaden the Z width. To illustrate our approach, we show how the resulting dijet constraints can be applied to simplified DM models. Moreover, we develop a new method to combine dijet constraints with information on the relic abundance of DM to determine those regions of parameter space where thermal freeze-out via a Z is incompatible with constraints from the LHC.
The outline of our work is as follows: We first perform a combined analysis of searches for dijet resonances from the CMS and ATLAS experiments at both 8 and 13 TeV in section 2. We present the results of this analysis in terms of an upper limit on the Z -quark coupling as a function of the Z mass and width, which applies to a wide range of Z models irrespective of what other couplings and particles are present in the model. In section 3 we then show how these results may be applied to a specific model, namely a model with a Z mediator coupling to quarks and DM. We first consider a simplified model similar to the one presently employed by the LHC collaborations and then combine such a model with relic density constraints to place limits on the coupling and the mass of the DM particle. Additional technical details are provided in the appendices.
2 Limits on generic Z models from dijet resonance searches at the LHC In this section we describe our technique for the re-analysis of searches for dijet resonances from the LHC run I and II and present the resulting combined limits. Our basic approach can be divided into three separate steps. We first implement a fully general Z model in a Monte Carlo generator to produce dijet events at the relevant centre-of-mass energies and apply the selection cuts corresponding to the various analyses. We then compare the predicted distributions of the dijet invariant mass to the experimental data, employing the same strategy as the experimental collaborations to model the background contribution. Finally, we combine the experimental tension from all data sets in a statistically consistent way in order to determine the largest signal strength that can still be compatible with all experimental data. The details of each of these steps are discussed in the following subsections.

Dijet event generation
The first step in our analysis is to generate dijet events resulting from a new Z mediator with mass 500 GeV ≤ m Z ≤ 4 TeV, which interacts with the standard model quarks q via the vectorial coupling g q . To generate this signal, we add the following terms to the SM Lagrangian: Although we do not specify any other couplings of the Z , such couplings may in principle be present. 1 We therefore do not calculate the total decay width Γ of the Z in terms of its quark coupling but instead take Γ to be a free parameter of the model. An important advantage of taking Γ and g q as independent parameters is that the shape of the dijet invariant mass distribution depends only on the two parameters m Z and Γ, whereas the total magnitude of the signal is proportional to g 4 q irrespective of whether the Z is produced on-shell or off-shell. 2 We can therefore generate events for different values of m Z and Γ and a fixed value of g q and apply the simple rescaling σ ∝ g 4 q to obtain signal predictions for the full three-dimensional parameter space.
Our simulation of dijet events is carried out using a pipeline of the publicly available software packages FeynRules_v1.6.11 [33], MadGraph_v3.2.2 [34], Pythia_v8.186 [35,36] and FastJet_v3.0.5 [37]. First, we implement the model Lagrangian in FeynRules to calculate the Feynman rules and generate a UFO model file [38]. In MadGraph we then generate matrix elements for all processes involving a virtual Z and a pair of u, d, s, c or b quarks in the final state.
The output from MadGraph is interfaced with Pythia, which we use both as a Monte Carlo event generator and to simulate showering and hadronisation. For our simulations, we use the CTEQ5L parton distribution function [39]. We neglect next-to-leading order corrections, which are expected to lead to somewhat larger cross sections [40], so we give a 1 In particular, we assume at this point that the Z has vectorial couplings to quarks. We will discuss below how our results can be applied to models with axial couplings to quarks and to a combination of vectorial and axial couplings. 2 Note in particular that we allow the unphysical situation that the branching ratio into quarks, given by Γ qq /Γ, can become larger than unity. This does not pose a problem as long as the resulting constraints are only evaluated for physical combinations of gq and Γ.  Table 1. Jet parameters chosen for the anti-k T algorithm for the five experimental searches. The radius parameter is defined as R = (∆η) 2 + (∆φ) 2 with φ the azimuthal angle. For the CMS analyses, jets are first reconstructed with a radius parameter of 0.5 (0.4) at 8 TeV (13 TeV) and are then combined into two fat jets with radius parameter 1.1.
conservative bound. The resulting final states are clustered with FastJet using the antik T algorithm [41]. We choose the jet parameters (cone-size R, maximum pseudorapidity η and minimum transverse momentum p T min of the jet) to match those adopted by each experiment under consideration (see table 1).
Once the jets have been reconstructed, we apply the experimental selection cuts outlined in table 2. For each event that passes these cuts we calculate the invariant mass of the dijet system. Rather than performing a full detector simulation, we approximate the uncertainties in reconstructing the energy and momentum of the jet events arising from detector performance by applying a Gaussian smearing to the invariant dijet mass. For this purpose, we take the detector resolution in both ATLAS and CMS to be σ(m jj ) = 1.8GeV m jj /GeV , (2.3) which was determined by fitting the smeared signals to shapes given by the CMS experiment [4] for a RS graviton benchmark model. The smeared invariant masses are then binned according to the bin sizes given by the different experiments and the resulting histograms are converted into differential cross sections dσ Z /dm jj .

Compatibility of a dijet signal with LHC data
Once the dijet invariant mass distributions have been generated, the second step is to determine the compatibility of such a signal with observations at the LHC. For this purpose, we follow the approach of the experimental collaborations and assume that the SM background can be described by a smooth function of the form:  Table 2. Experimental cuts adopted by the five experimental searches. p T,j1 refers to the transverse momentum of the leading jet, while p T,j2 refers to the subleading jet. H T is the scalar sum of all jet p T for jets with p T > 40 GeV and |η| < 3, and ∆η jj refers to the rapidity separation of the leading and subleading jets.
the first term depends on the unknown parameters P i , while the second term depends on the assumed values for m Z , Γ and g q .
To compare the model prediction to experimental data, we calculate the usual χ 2 test statistic where the index i denotes the bin number in a given experiment, d i is the observed differential cross section with corresponding (statistical) error σ i and s i is the predicted signal containing both the SM contribution and the new-physics signal. We now fix the unknown parameters P i by finding the minimum of the χ 2 distribution with respect to these parameters (calledχ 2 ). We now want to place an upper bound on the magnitude of the new-physics signal, beyond which the sum of signal and background are incompatible with the data. For this purpose, we employ a ∆χ 2 method. We first calculateχ 2 in the absence of a contribution from the Z mediator (calledχ 2 0 ) and then define ∆χ 2 (m Z , Γ, g q ) =χ 2 (m Z , Γ, g q ) −χ 2 0 . For ∆χ 2 < 0, the data actually prefers a non-zero contribution from the Z mediator. Positive values of ∆χ 2 , on the other hand, are disfavoured by the data. In such a case, we can calculate the p-value, i.e. the probability to observe at least as large a value of ∆χ 2 from random fluctuations in the data as where CDF 1, χ 2 is the cumulative distribution function for the χ 2 distribution with one degree of freedom. 4 4 We note that the ∆χ 2 test statistic as we define it does not exactly follow a χ 2 -distribution. The reason is that we takeχ 2 0 to be the value ofχ 2 for gq = 0 rather than finding the value of gq that actually minimiseŝ χ 2 (called gq,0) in order to avoid the problem that the data may prefer a negative signal contribution. Sincê χ 2 (0) ≥χ 2 (gq,0), our definition yields slightly smaller values for ∆χ 2 than the one obtained from minimisinĝ χ 2 with respect to gq. Using a χ 2 -distribution to calculate the p-value therefore means that we slightly overestimate the p-value and consequently place more conservative bounds. We have verified that the error made by this approximation is small by determining the actual distribution of ∆χ 2 from a Monte Carlo simulation for specific parameter points. σ ⨯ BR ⨯ A [pb] Our limit CMS observed limit CMS expected limit Figure 1. Comparison of our method for setting limits from dijet data (blue, dashed) with those from the CMS experiment [4] (green, solid) presented in terms of the cross section times acceptance, for the benchmark model of an RS graviton [6].
As discussed above, the new-physics signal is proportional to g 4 q . As we increase g q , keeping m Z and Γ fixed, we will reach the point where ∆χ 2 becomes so large that P becomes unacceptably small. For P < 5%, we can exclude the corresponding value of g q at the 95% confidence level. If P > 5%, the value of g q cannot be excluded by the experiment under consideration, but it may still be excluded by the combination of results from several experiments. Such a combination is necessarily model-dependent in the sense that it requires an assumption on the ratio of the production cross section of the resonance at 8 TeV and 13 TeV. Since we use a Z -model to generate dijet events, our combination is valid for any resonance produced dominantly from light quarks (with equal couplings to all flavours).
For a given signal hypothesis, we can follow the procedure described above to obtain a value of ∆χ 2 for each experiment under consideration. Crucially, the parameters P i are fitted independently for each experiment. Since the two CMS searches at 8 TeV are not statistically independent, we use ref [4] for m Z ≥ 1 TeV and ref [5] for smaller Z masses. We can then simply add up the individual contributions to ∆χ 2 to obtain ∆χ 2 total . This test statistic is again expected to approximately follow a χ 2 -distribution with one degree of freedom, so we can calculate the combined p-value with eq. (2.6).

Validation based on the RS-graviton model
To validate our limit-setting procedure, we have applied our analysis chain to the RS graviton model [6], which is used as a fiducial benchmark by CMS [4]. The model contains two free parameters, namely the mass of the RS graviton m RS and the curvature of the five-dimensional bulk k/M Pl whereM Pl is the reduced Planck mass. The latter is taken to be k = 0.1M Pl , which fully determines the width and the couplings relevant for the generation of dijet events for a given value of m RS . We adopt these values to calculate the total production cross section and to generate dijet invariant mass distributions. We then multiply the resulting distributions with a rescaling factor µ in order to determine the largest signal strength that is still compatible with data. Applying the resulting rescaling factor to the total cross section then yields an upper bound on the production cross section as a function of the RS graviton mass. Figure 1 shows the comparison of the bounds we obtain with the results from the CMS analysis. We conclude that our approach yields good agreement with the results from the CMS collaboration over a wide range of masses.

Results
Having described the strategy for deriving bounds on the Z -quark coupling g q , we now present the results of our analysis in a way that can be applied to a wide range of Z models. We consider m Z masses between 500 GeV and 4 TeV in steps of 50 GeV. For each mediator mass, we consider six different widths: 1%, 2%, 5%, 10%, 20% and 30% (in units of m Z ). These values were chosen to best sample the variation in the constraints for different widths. In particular, we note that for a Z width smaller than 1% of m Z , the shape of the dijet invariant mass distribution is dominated by the detector resolution and therefore becomes independent of Γ. For each combination of m Z and Γ we then determine the largest value of g q that is compatible with the experimental data at 95% confidence level (called g q,95% ).
While the resulting values of g q,95% typically depend on m Z in a non-monotonic way (due to different random fluctuations from bin to bin), the dependence on Γ is typically very smooth, with larger values of Γ corresponding to larger values of g q,95% (i.e. weaker bounds). We can make use of this observation to interpolate between the values of Γ considered in our simulation. Specifically, we find that it is possible to fit g q,95% for fixed m Z using a function of the form where the values of a, b and c are listed in table 3 in appendix A as a function of m Z . We show g q,95% as a function of m Z and Γ as obtained from the interpolation functions in the left panel of figure 2. For m Z 1.5 TeV dijet constraints are able to exclude values of g q between 0.1 (for a narrow width) and 0.3 (for a broad width). For larger masses, these bounds become somewhat weaker and reach up to g q,95% ≈ 0.6 for m Z ∼ 4 TeV and Γ/m Z > 0.2. We observe that rather weak bounds are obtained for m Z ≈ 1.6-1.7 TeV. The reason is that in this mass range all four experiments see an upward fluctuation in the data, so that the observed bound is weaker than the expected one (see also [42,43]). 5 Since we have consistently treated Γ and g q as independent parameters, our results can be applied to any Z model (with universal vector-like couplings to quarks) by applying the following procedure: 1. For given Z mass and given couplings of the Z to all other particles in the theory, calculate the total decay width Γ.
2. Look up g q,95% for this value of Γ and the assumed Z mass.
3. If g q,95% is larger than the assumed Z -quark coupling, the parameter point is allowed. Otherwise, it is excluded at 95% confidence level.
While the procedure detailed above applies to Z models with universal vector couplings to all quarks, it is also possible for us to constrain more complicated models. For this purpose, we make use of the narrow-width approximation (NWA), which is valid as long as the width of the Z is small compared to its mass (typically Γ/m Z < 0.3). The NWA states that the cross section for the production of dijet events via a resonance factorises into the production cross section of the resonance and the probability for this resonance to decay into a pair of jets: In the model we consider the Z production cross section is proportional to g 2 q , with a constant of proportionality that depends on the Z mass and the centre-of-mass energy. Consequently, the dijet signal in each experiment is proportional to g 2 q times the relevant branching ratio: Indeed, this relation is also correct for the differential cross section, i.e. the shape of the dijet invariant mass distribution is independent of the coupling g q for fixed mediator mass and width. This observation motivates a different way of presenting our results, namely to place an upper bound on the combination [13] j ≡ g 2 q × BR(Z → jj) . (2.10) As discussed above, j is proportional to g 4 q for fixed Γ. We can then calculate the upper bound on j at 95% confidence level, called j 95% , by evaluating j for g j,95% . We emphasise that it is perfectly acceptable for this calculation to yield a branching ratio larger than unity. In this case the conclusion would simply be that the experimental bounds cannot exclude any value of g q compatible with the chosen value of Γ.
Our results for j 95% are shown in the right panel of figure 2. The advantage of this approach is that j 95% can also be used to constrain models beyond the one considered here. In particular, our analysis can be applied to the following cases: • For a Z with both vector (g V q ) and axial (g A q ) couplings to quarks, the production cross section for a Z is proportional to (g V q ) 2 + (g A q ) 2 . In such a model, one should therefore calculate (g V q ) 2 + (g A q ) 2 × BR(Z → jj) and compare the result to j 95% as shown in the right panel of figure 2. • The Z production is typically dominated by up and down quarks in the initial state.
Consequently, for a Z with different couplings to the three generations, one can obtain an approximate bound by calculating g 2 1 × BR(Z → jj), where g 1 is the coupling to the first generation, and comparing the result to the bound on g 2 For the convenience of the reader, we provide a plain text version of j 95% as a function of m Z and Γ in the supplementary material accompanying this paper. The advantage of these bounds is that they are independent of other interactions that may be present in a given Z model.
This discussion concludes the general analysis of dijet constraints. The remainder of the paper is dedicated to a specific application of the results shown above, which serves as an illustration for the procedure described above and uses this procedure in order to constrain a model of particular interest.

Constraints on a leptophobic Z coupling to DM
We now show how the results from the previous section can be applied to a specific model. For this purpose, we consider a simple model of a leptophobic Z coupling to DM (see [44][45][46][47]), which is similar in spirit to the spin-1 s-channel simplified DM model discussed in refs. [18,20,23,26,27,29,30]. Assuming DM to be a Majorana fermion ψ, the interactions between the SM and the dark sector are defined by the following Lagrangian: The Majorana nature of the DM particle ensures that it can only have an axial coupling to the Z , which significantly reduces constraints on the model from direct detection experiments. On the SM side, the couplings of the Z are assumed to be purely vectorial, which is consistent with the assumption that the Z couples neither to leptons nor to the SM Higgs [48]. We do not specify the additional dark Higgs necessary to generate the Z mass and the DM mass [48], assuming that this particle is sufficiently heavy and sufficiently weakly mixed with the SM Higgs to be irrelevant for LHC phenomenology. While additional heavy fermions are needed to cancel anomalies, these can be colour-neutral and vector-like with respect to the SM gauge group, making them very difficult to observe at the LHC (see [12] for a discussion of anomaly-free models).
The resulting decay widths are then given by 7 where we have assumed m Z > 2m t , 2m DM . The equations above enable us to calculate the total decay width Γ, which is required in order to apply the dijet bounds derived above.

Bounds for fixed couplings
The model introduced above has four free parameters (the two masses m Z and m DM and the two couplings g q and g DM ). Since kinematic distributions at the LHC depend more sensitively on the masses than on the couplings, it is interesting as a first step to study LHC constraints for fixed couplings and varying masses. This approach is consistent with the most common way of presenting LHC constraints on DM simplified models.
Note however that our model is not identical to any of the simplified models presently used by the ATLAS and CMS collaborations, because we consider a Majorana DM particle and a different coupling structure (in our model the Z has vector couplings to quarks and axial couplings to DM). Using a Majorana fermion instead of a Dirac fermion leads to an invisible width smaller by a factor of two giving slightly stronger bounds from dijet searches. The different coupling structure is not expected to significantly alter dijet bounds which typically depend on (g V q ) 2 + (g A q ) 2 (see above). It does, however, change the relic density constraint compared to the one obtained for the simplified models used by the experimental collaborations.
Following the recommendations from [30,49], we consider the case q q = 0.25, g DM = 1. For these couplings the width of the Z varies between 2.5% (for m Z < 2m DM , 2m t ) and 4.3% (for m Z 2m DM , 2m t ) of its mass. For each combination of m Z and m DM , we calculate the width Γ and then read off the largest allowed value for g q from figure 2. Whenever g q,95% < 0.25, the parameter point is excluded by our combined dijet bounds. The results of this analysis are shown in figure 3 with the dijet excluded regions shown 7 The pre-factor 1/(24π) for the decay into DM results from the fact that there are two identical particles in the final state.  Figure 3. Excluded regions of parameter space in the mass-mass plane for fixed couplings, following the recommendations of the DM LHC working group [49]. The region in red is excluded by our combined dijet analysis at the 95% confidence level while the green lines represent parameter points which reproduce the observed relic density of DM in the Universe. In the grey region perturbative unitarity is violated.
in red. We find that Z masses between 500 GeV and 1600 GeV are excluded irrespective of the value of m DM . For m Z between 1600 GeV and 3 TeV, the model is excluded for heavy DM particles, such that the invisible branching ratio of the Z is small and decays into dijets dominate. These bounds are somewhat stronger than the ones found by the individual experiments due to our combined analysis. For comparison, we also show the parameter region (in grey) where the coupling of the DM particle to the longitudinal component of the Z violates perturbative unitarity [25,48] and the masses for which the assumed couplings reproduce the observed DM relic abundance, Ωh 2 = 0.12 [50], calculated using micrOMEGAs_v4.1.8 [51]. Details on the relic density calculation can be found in the following subsection.
We find that there are only two regions where the relic density is compatible with dijet constraints: A low-mass region with m Z < 500 GeV and m DM < 150 GeV and a high-mass region with m Z > 3 TeV and m DM > 1200 GeV. 8 This conclusion is, however, obviously dependent on our choice of couplings. For example, smaller values of g q would reduce the production cross section of the resonance, while larger values of g DM would increase the invisible branching ratio (provided m DM < m Z /2), so that dijet constraints could be significantly relaxed. 8 We emphasise that in the low-mass region there may be additional dijet constraints from previous hadron colliders, as well as constraints from dijet resonances produced in association with SM gauge bosons [11,25]. Moreover, this region of parameter space is tightly constrained by mono-X searches, in particular searches for jets in association with missing transverse energy [52][53][54][55]. These searches are very sensitive to Z masses below about 1 TeV, but lose sensitivity very quickly towards larger masses, where dijet constraints can still be sensitive.
To study whether the intermediate mass range can be made compatible with dijet constraints for different choices of couplings, one could in principle repeat the analysis from above for many different combinations of g q and g DM (or simply scan the entire parameter space). Instead, we will take a more systematic approach and develop a new method that can be used to establish the compatibility of relic density constraints and dijet searches across the entire parameter space of our model (a similar comparison between LHC searches for dilepton resonances and the DM relic abundance in the context of gauged B − L has been performed in [56]).

Combining di-jet bounds and relic density
Out of the four-dimensional parameter space of our model, we are particularly interested in those combinations of masses and couplings for which the thermal freeze-out of the DM particle can reproduce the relic abundance which is the result from Planck CMB observations combined with Baryon Acoustic Oscillations, supernova data and H 0 measurements [50]. We will approximate the relic density as Ωh 2 = 0.12 in the rest of this work. We emphasise that the relic density requirement can be relaxed if the dark sector consists of multiple components or if the thermal history of the Universe is non-standard. Nevertheless, it is certainly of interest to consider those parameters for which the simplest assumptions are already sufficient to match observations. If these parameters can be excluded experimentally, the model would require additional ingredients in the dark sector (such as additional annihilation channels, additional stable states or a mechanism to produce additional entropy after DM freeze-out), which by itself would be an important conclusion.
The remainder of this section focusses on how to reduce the parameter space of our model by imposing the relic density constraint. We first discuss some general aspects of the relic density calculation and then introduce a convenient set of free parameters that can be used to combine the relic density requirement with dijet constraints. Finally, we apply the dijet constraints from above to place bounds on the simple thermal freeze-out scenario.
To first approximation, we can obtain the relic density by calculating the cross section for DM annihilation into a pair of quarks, σ ψψ→qq , and expanding the result in terms of the relative velocity v of the two DM particles: The relic abundance is then approximately given by where x fo ∼ 20-30 is the ratio of the DM mass and the freeze-out temperature and g * ∼ 80-90 is the number of relativistic degrees of freedom during freeze-out. For our model, we For m DM ≈ m Z /2, the denominator in eq. (3.8) becomes very small and DM annihilation receives a resonant enhancement. In this case, an expansion in terms of the velocity of the two DM particles is insufficient for an accurate calculation of the relic density and numerical methods are needed. We therefore calculate the relic density using micrOMEGAs_4.1.8 [51], including two modifications under the instruction of the authors (see appendix B). 9 For given m Z , m DM and g DM we can then numerically determine the value of g q that is required to reproduce the observed relic density. 10 As long as m DM is well below the resonance region, i.e. m DM m Z /2, eq. (3.8) implies that the annihilation cross section is proportional to g 2 q g 2 DM m 2 DM /m 4 Z . Therefore it is always possible to fix g q in such a way that the observed value of Ω h 2 is matched and the solution is always unique. In the resonance region, the annihilation cross section is proportional to g 2 q g 2 DM /Γ, which is still a monotonic function of g q so that any solution is unique. However, since the expression g 2 q g 2 DM /Γ remains finite for g q → ∞, it is possible that no solution exists. In short, as long as g DM is large enough, there will always be a unique value of g q that reproduces the relic abundance.
We therefore obtain a function g q (m Z , m DM , g DM ), which is illustrated in figure 4 as a function of g DM for various fixed values of m Z and m DM . The resulting curves have the following features: 1. Since the annihilation cross section grows monotonically with the DM mass (for fixed couplings and mediator mass), the lines for different DM masses never cross, i.e. smaller values of m DM always require larger couplings.
2. For sufficiently small DM masses, the curves are hyperbolas (g q ∝ 1/g DM ), whereas for larger values of m DM , the curves are steeper at small g DM and flatter at large g DM due to the resonance effects discussed above.
We will make use of these properties below to choose a particularly convenient set of free parameters for the analysis of our model.

Relic density constraints for a fixed width
Having constructed the function g q (m Z , m DM , g DM ) from the relic density requirement, one could now simply proceed to scan the remaining three-dimensional parameter space. One subtlety arises, however, from the fact that the width Γ -and therefore the bound from dijet constraints -depends on all three parameters in a non-trivial way. For example, for fixed m Z and m DM one would naively expect stronger dijet constraints for smaller g DM corresponding to larger g q (implying both a larger production cross section of the resonance and a larger branching fraction into dijets). However, if at the same time Γ increases, it is conceivable that dijet constraints are weakened sufficiently to evade experimental bounds and that in fact larger values of g q are less constrained than smaller couplings.
To avoid this complication, we take both m Z and Γ as free parameters. As shown in figure 2, for fixed values of these two parameters we can always place an unambiguous upper bound on g q . A second important advantage of this approach is that m Z and Γ are the two parameters that are most directly observable at the LHC. While the DM mass is very difficult to measure at the LHC and coupling constants can only be inferred in the context of a specific model, an observation of a new resonance in the dijet channel would immediately enable us to determine the mass and the width of the mediator from the invariant mass distribution. 11 To be able to directly interpret such an observation in the context of the present model, it therefore makes sense to construct all bounds in terms of these two most apparent observables.
In order to treat the width Γ as a free parameter, we need to determine those combinations of m DM , g DM and g q that reproduce the observed relic density while at the same time matching the required width. For this purpose we first of all observe from eqs. (3.3)-(3.4) that for fixed m Z and m DM < m Z /2 the total width Γ is an ellipse in the couplings. 12 We can now consider ellipses of constant width Γ in the same g q -g DM -plane used in figure 4 to study the relic density constraints. Since the relic density curve is convex while the constant-width curve is concave, the two curves will have either exactly two intersects or zero intersects (neglecting those special cases where the two curves just touch at exactly one point). In other words, for fixed values of m Z , Γ and m DM , there is either no combination of g q and g DM that reproduces the relic density constraint or there are two separate solutions corresponding to the desired value of Γ. Whenever there are two solutions, we 11 This argument assumes that the width of the resonance is large compared to the detector resolution.
Nevertheless, for a narrow resonance it is still possible to determine the mass and place an upper bound on the width. 12 For mDM ≥ m Z /2, the total width is independent of gDM and hence a straight line gq = const.  Figure 5. Examples of how we find pairs of couplings that satisfy the relic density constraint (blue) for a given fixed width (red).
define Solution I to be the one with larger g DM (and therefore smaller g q ) and Solution II to be the one with smaller g DM (larger g q ). Two examples are shown in figure 5.
As noted above, increasing the value of m DM will shift the relic density curve towards smaller couplings. Conversely, the constant-width curve will be shifted to larger couplings (due to the larger phase-space suppression of Z → ψψ). This means that for each value of m Z and Γ there is a minimum value of m DM , called m DM,min (m Z , Γ), such that there is no solution for m DM < m DM,min and two solutions for m DM > m DM,min . Increasing m DM beyond m DM,min will shift Solution I to larger values of g DM and smaller values of g q and vice versa for Solution II. As m DM approaches m Z /2, Solution I will yield arbitrarily large values of g DM and thus ultimately violate the perturbativity bound g DM < √ 4π. Solution II, on the other hand, will approach a small but non-zero minimum value of g DM , called g DM,min . 13 Figure 6 shows both m DM,min and g DM,min as a function of m Z and Γ.
With these considerations in mind, we can now eliminate either m DM or g DM in favour of Γ and proceed with either (m Z , Γ, g DM ), where g DM,min < g DM < √ 4π, or with (m Z , Γ, m DM ), where m DM,min < m DM < m Z /2. While in the first case we find a single value of g q for each set (m Z , Γ g DM ) compatible with the relic density constraint and consistent with the required width, the second case yields two separate solutions as discussed above. We discuss both possibilities below, as they each offer different physical insights.

Dijet bounds on the DM coupling
We have shown above that for fixed m Z and Γ smaller values of g DM correspond to larger values of g q . Using figure 2 to read off the upper bound on g q from dijet searches thus allows us to place a lower bound on g DM . This lower bound on g DM is shown in figure 7. Wherever no bounds from dijet searches can be placed, we simply show the smallest value of g DM for which the relic density curve and the constant-width curve intersect (called g DM,min above). If on the other hand the lower bound from dijet searches is so strong that it requires g DM to be larger than √ 4π, we conclude that it is impossible to find perturbative values of g q and g DM such that the width Γ and the relic density can be reproduced without violating dijet constraints. The corresponding regions are shaded in orange in figure 7.
We observe that rather large values of g DM are required in order to avoid dijet constraints. While the consistency of the relic density requirement and the assumed width only required g DM,min ∼ 0.1-0.3 (see figure 6), dijet constraints require g DM > 1 in almost the entire parameter space that we consider. For large Z width, even larger values of g DM are required in order to reduce the branching ratio of the Z into dijets. For m Z 1.5 TeV and Γ/m Z 0.2 as well as for 1.7 TeV m Z 3.3 TeV and Γ/m Z 0.25, all perturbative values of g DM that reproduce the relic abundance are excluded by dijet searches. For larger Z masses, LHC dijet searches lose sensitivity, but significant improvements in this mass range can be expected from upcoming runs of the LHC at 13 TeV.

Dijet bounds on the DM mass
Let us finally present our results from a complementary perspective by taking m DM as a free parameter and determining both g DM and g q from the relic density constraint and the requirement of a constant width. As discussed above, for each value of m DM we obtain two separate solutions, with Solution I (II) corresponding to larger (smaller) g DM . For each of the two solutions, we can directly read off from figure 2 whether the parameter point is excluded by the combined dijet constraints that we have derived above. These exclusion limits in turn allow us to determine the allowed range of DM masses as a function of m Z and Γ. We now discuss the two different solutions in turn. As noted above, for Solution I (i.e. larger values of g DM ), the DM coupling increases with the DM mass. The requirement to have a perturbative coupling, g DM < √ 4π, therefore gives an upper bound on m DM , called m DM,max . For some values of m Z and Γ we find that Solution I yields a non-perturbative value of g DM for all values of the DM mass, so for these combinations of Z mass and width only Solution II is of interest.
Conversely, for Solution I smaller DM masses correspond to larger values of g q . Since large values of g q are excluded by LHC dijet searches, we can use the LHC bounds to place a lower limit on m DM , called m DM,min . 14 The combination of the perturbativity requirement and LHC dijet searches therefore yield a range of permitted dark matter masses [m DM,min , m DM,max ], which satisfy all of our constraints. In other words, for m DM in this range, it is possible to find values of g q and g DM that yield the observed relic abundance and are consistent with all other constraints that we consider. Figure 8 shows one example, where we have fixed the Z mass to 1 TeV and show m DM,min (orange) and m DM,max (blue, dashed) as a function of Γ. To illustrate the impact of dijet searches, we also show the value of m DM,min that one obtains solely from the consistency of relic density and constant width (orange, dotted). For large values of Γ we find that m DM,max < m DM,min , i.e. all perturbative solutions are excluded by LHC dijet searches. In the specific case under consideration, this occurs for Γ/m Z 0.19. Figure 9 shows the largest allowed DM mass (left) and the smallest allowed DM mass (right) as a function of m Z and Γ. The plots can be read by picking a point on the plane (i.e. fixing m Z and Γ) and then reading of m DM,max and m DM,min from the two panels to find the range of permitted dark matter masses [m DM,min , m DM,max ] that satisfy all constraints. Those combinations of m Z and Γ for which m DM,max < m DM,min are shaded in orange. The  grey region indicates those combinations of Z mass and width for which no perturbative solutions can be found. As expected, the orange shaded region is identical to the one found in figure 7.
Turning now to Solution II, we note that for this solution perturbativity constraints will typically be less important (because we consider smaller values of g DM ), while dijet constraints will be more important (because the corresponding values of g q are larger). 15 Compared to the previous solution, the situation is now reversed: The requirement of perturbativity may raise m DM,min , while dijet constraints will lower m DM,max . We show the maximum and minimum allowed DM masses for Solution II in figure 10.
As expected, we find dijet constraints (shown in orange) to be significantly stronger than for Solution I. For large width, the entire range 500GeV ≤ m Z ≤ 3500GeV is excluded  by dijet constraints. For m Z ∼ 1200 GeV the dijet bounds extend down to very narrow resonances. As discussed above, dijet bounds are particularly weak around 1600 GeV, due to an intriguing upward fluctuation in the data. Finally, we note that we can always find a value of the DM mass such that Solution II corresponds to perturbative couplings (so the grey shaded region from figure 9 is absent).

Conclusions
We have presented a combination of all available searches for dijet resonances at the LHC in the context of a generic Z model. Taking the width of the resonance and its coupling to quarks as independent parameters allows us to obtain constraints that apply irrespective of whether the Z decays exclusively into quarks or dominantly into other states. The results of this analysis, summarised in figure 2 and table 3, are provided in such a way that they can be easily used to constrain a range of different models.
As a specific illustration of our approach, we have applied our constraints to a Z that couples to quarks and dark matter (DM), similar in spirit to a DM simplified model with a spin-1 s-channel mediator. It is straight-forward to map our constraints onto the parameter plane showing DM mass versus mediator mass for fixed couplings, which is conventionally used to present LHC results from missing energy searches. We show that for the typical choice of couplings (g DM = 1, g q = 0.25), dijet searches can exclude the range 500 GeV < m Z < 3 TeV for almost all values of the DM mass (see figure 3). These findings suggest that future searches for simplified DM models should focus on smaller values of g q and larger g DM , which would relax constraints from searches for dijet resonances while still allowing for sizeable interactions between DM and quarks.
Finally, we have focussed on the special case that the Z mediates the interactions of DM and quarks responsible for thermal freeze-out, so that one of the parameters of the model can be eliminated by the requirement to reproduce the observed relic abundance. We have constructed a novel way of studying this set-up by making explicit the parameters that can be directly probed by searches for dijet resonances, i.e. the mass and the width of the Z . The remaining free parameter can then be taken to be either the DM coupling (figure 7) or the DM mass (figures 9 and 10). We find that for very broad widths (Γ/m Z 0.25) and Z masses below about 3 TeV, LHC searches already exclude the possibility that the DM-quark interactions mediated by the Z are responsible for setting the DM relic abundance.
Furthermore, these figures provide a useful tool for interpreting future searches for dijet resonances at the LHC. Should an excess be seen in such a search, the mass and the width of the resonance can be determined from the data in a model-independent way. One can use these figures to look up whether the new state could conceivably act as the mediator into the dark sector. If a solution to the relic density requirement exists, the plots then provide the allowed ranges of the DM mass and coupling. Presently there is still ample room for such an interpretation, so there is much to be learned from the upcoming LHC data at 13 TeV.
where R is an arbitrary number that is by default fixed to 2.7. For |q 2 − m 2 | > √ R 2 + 1 m Γ, on the other hand, the width in eq. (B.1) is set to zero. In the intermediate region the width Γ is replaced by a function of q 2 that interpolates between the two cases. For fairly large widths, as considered in the present work, this interpolation procedure can lead to kinks in the relic density calculation. To remove such kinks one can simply increase the value of R from its default value by changing the value of the variable BWrange [58]. We have found that R = 100 is sufficient to remove the kinks in our plots. 16 Furthermore, as pointed out in ref. [31], the width Γ of the resonance depends in general on the momentum transfer q 2 , i.e. Γ = Γ(q 2 ). For the case of narrow widths (Γ/m 1), eq. (B.1) gives a good approximation, because Γ is only relevant for q 2 ≈ m 2 and one can therefore approximate Γ ≈ Γ(q 2 = m 2 ). Since we consider widths as large as 30% in this work, it is however appropriate to modify the BW formula in order to take the momentum dependence of the width into account.
Following appendix A of [31], we can approximately capture the momentum dependence by setting Γ(q 2 ) = √ q 2 m Γ(q 2 = m 2 ). This substitution yields for the shape of the BW resonance. This modification can be implemented by editing the function prepDen used in the CalcHEP code.