Higher-order QCD predictions for dark matter production at the LHC in simplified models with s-channel mediators.

Weakly interacting dark matter particles can be pair-produced at colliders and detected through signatures featuring missing energy in association with either QCD/EW radiation or heavy quarks. In order to constrain the mass and the couplings to standard model particles, accurate and precise predictions for production cross sections and distributions are of prime importance. In this work, we consider various simplified models with s-channel mediators. We implement such models in the FeynRules/MadGraph5_aMC@NLO framework, which allows to include higher-order QCD corrections in realistic simulations and to study their effect systematically. As a first phenomenological application, we present predictions for dark matter production in association with jets and with a top-quark pair at the LHC, at next-to-leading order accuracy in QCD, including matching/merging to parton showers. Our study shows that higher-order QCD corrections to dark matter production via s-channel mediators have a significant impact not only on total production rates, but also on shapes of distributions. We also show that the inclusion of next-to-leading order effects results in a sizeable reduction of the theoretical uncertainties.

Introduction

Various astrophysical and cosmological observations provide strong hints for the existence of dark matter (DM). Yet, very little is known about the nature of DM, and constraints on models from various direct/indirect detection experiments and cosmology still allow for a wide range of DM masses and couplings to the Standard Model (SM) particles. Various types of DM searches are sensitive to different regions of the DM model parameter space [1, 2]. In order to maximise the chances for discovering – or at least excluding – DM models, it is hence imperative to perform both astrophysical and collider searches for DM. The most promising way to detect signals of weakly interacting DM particles at the LHC is through their associated production with jets, EW bosons and heavy quarks, leading to signatures with missing transverse energy (MET). Searches for DM have been performed at LHC Run I (see e.g. Refs. [3, 4]) and are one of the central goals of LHC Run II [5].

While the complementarity of different DM searches is a powerful tool, it is intrinsically limited in that the comparison of collider and other results introduces some degree of model dependence. As the nature of DM is still unknown, there exist a myriad of DM models and mechanisms to be explored, spanning a wide range of complexity and ambition. In this context, it is of utmost importance to follow an approach where the model dependence is limited, while the salient features that can provide a useful characterisation of possible signals are kept intact, i.e., the “simplified models” approach [6].

Many simplified models for DM have been proposed in the past (see Refs. [5, 7, 8] and references therein). In their simplest form, these models assume DM to be a single massive particle which interacts weakly with the SM particles. The interaction of DM with the SM can be mediated by a new field, which we dub the mediator. When the mediator is heavier than the energy scales the experiment can probe, the interaction becomes point-like and the Lagrangian can be organised in terms of a tower of higher-dimensional operators in the framework of an effective field theory (EFT). However, when the experimentally accessible energy becomes comparable or higher than the mediator mass, on-shell effects become important and a properly defined quantum field theory featuring the mediator state in the spectrum is needed [9-14]. The LHC can explore a large range of DM and mediator masses, as well as coupling strengths and possible combinations of DM and mediator spins. Collider results within the framework of simplified models can then be combined with direct searches for the mediators, e.g. in Drell–Yan processes, as well as with cosmological and astrophysical constraints on DM.

LHC searches for DM will rely on precision calculations to impose the most stringent bounds on DM models, and hopefully characterise possible signals. Higher-order corrections in QCD to DM production processes at the LHC are hence vital for extracting precise information about DM from the LHC results.

Previous analyses [15-19] studied next-to-leading order (NLO) QCD corrections for DM production in MET in the case of EFT, i.e., in the limit of heavy mediators.1 In this article, we consider simplified models with s-channel mediators. We analyse the impact of the higher-order corrections to mono-jet signals in various benchmark scenarios with spin-1 (vector and axial-vector) mediators, calculate DM production cross sections (both total and differential) at NLO accuracy with up to two jets, and also merge the corresponding samples via the FxFx procedure [20]. To our knowledge, such accuracy has so far not been achieved in the context of DM simulations/predictions for the LHC Run II.

In addition, we consider associated production and compute NLO cross sections and distributions for spin-0 (scalar and pseudo-scalar) mediators in representative cases, including those with a light mediator. Predictions for this class of processes at NLO in QCD also represent a novelty in the context of DM computations for LHC Run II.

The first goal of this work is to illustrate the feasibility of having a fully general implementation of DM simplified models in the FeynRules [21]/MadGraph5_aMC@NLO [22] (MG5aMC henceforth) framework, accurate up to NLO in QCD. To this aim we start with the simplest (yet non-trivial) case of s-channel mediators (colour singlet, spin 0 and 1 bosons) coupling to DM and quarks. We assume that DM is a Dirac fermion for concreteness, yet our implementation is not limited by the choice of DM spin or nature (Dirac or Majorana). We show how predictions and event generation for this class of models can be achieved at NLO QCD accuracy, in a fully automatic way, for a wide set of observable/final states, while also employing the most recent matching/merging multi-jet techniques [20].

The second goal of this work is of phenomenological nature, i.e. to investigate the impact of the NLO corrections on the current and future searches for DM at the LHC. We consider two examples, among several promising final state signatures:for a spin-1 mediator model, andfor a spin-0 mediator model, where X is a DM particle. We do not only calculate NLO QCD corrections to the overall production rates, but also study in detail the effects of higher-order corrections on the differential distributions of relevant observables.

Our simulation set-up is based on the FeynRules/MG5aMC framework. The FeynRules package provides the relevant Feynman rules starting from any local Lagrangian [21, 23, 24], as well as the UV/ counterterms [25] necessary for the NLO computations via Nloct [26]/FeynArts [27]. Our simplified DM model files are publicly available at the FeynRules repository [28]. With these ingredients, which are only based on the model and are not process dependent, MG5aMC computes tree-level amplitudes, loop-amplitudes [29-31] and subtraction terms for a desired process, as well as their integration over phase space [32]. Event generation is obtained by matching short-distance events to the parton shower employing the MC@NLO method [33], which is implemented for Pythia6 [34], Pythia8 [35], Herwig6 [36] and Herwig++ [37]. We note that our DM model files can be exported not only to event generators but also to tools for DM relic abundance as well as DM direct and indirect detections such as MicrOMEGAs [38, 39] and MadDM [40, 41], allowing for more comprehensive DM studies.

The paper is organised as follows. In Sect. 2 we introduce simplified models for DM and specify the relevant interactions and model parameters of the implementation. We discuss the impact of the NLO QCD corrections on DM pair production with jets through spin-1 mediators in Sect. 3, which includes a discussion of inclusive cross sections, differential distributions as well as a discussion of the impact of parton showers and the NLO merging of events with different jet multiplicities. NLO QCD predictions for DM production in association with a top-quark pair through spin-0 mediators are discussed in Sect. 4. We provide our conclusions and an outlook in Sect. 5.

Simplified dark matter models: the s-channel mediator case

We start by defining the particle content and the interactions of the simplified model, which we dub DMsimp. We assume that DM is described by a single, massive and weakly interacting particle, that communicates with the SM through the exchange of a mediator. For simplicity, we assume that the mediator is not part of the SM.2 The first very general classification stems from the class of vertices that characterise the model: Lagrangians featuring DM–DM-mediator and SM–SM-mediator type interactions identify models with an s-channel mediator, while Lagrangians characterised by DM–SM-mediator interactions define a t-channel mediator. The former interactions, for example, arise in models featuring an extra scalar or which couples to a pair of DM particles, while the latter is common in supersymmetric models. From the point of view of QCD corrections the two classes are very different, as an s-channel mediator is necessarily a colour singlet, while a t-channel mediator can be either neutral or coloured. In this work we focus on the s-channel models, leaving the implementation and validation of t-channel models to forthcoming work. The s- and t-channel classes can be further organised by the quantum numbers of the DM particle and the mediator. To start with, we focus on the case of Dirac DM with spin-1 or spin-0 mediators coupling to the matter fields of the SM. Changing the spin or the nature of the fermion (Dirac or Majorana) of the DM particle or including a coupling of the mediator to the SM bosons is straightforward [42]. On the other hand, extending our analysis to spin-2 mediators, while feasible in principle, entails dedicated validation work, as such models are, in general, not renormalisable. We defer such an extension to the future.

Spin-1 mediator model

In the framework of our simplified model, the interaction Lagrangian of a spin-1 mediator () with a Dirac fermion DM () is given byand with quarks bywhere d and u denote down- and up-type quarks, respectively, () are flavour indices, and are the vector/axial-vector couplings of DM and quarks. Note that we adopt this notation according to the actual implementation in FeynRules. The model file, including an alternative choice for the spin of DM particle (complex scalar ), can be downloaded at the FeynRules repository [28].

The pure vector and pure axial-vector mediator scenarios are given by setting the parameters in the Lagrangians (3) and (4) toandrespectively, where we assume quark couplings to the mediator to be flavour universal and set all flavour off-diagonal couplings to zero. With this simplification of a single universal coupling for the SM- interactions, the model has only four independent parameters, i.e. two couplings and two masses:We note that the mediator width is calculated from the above parameters.

Finding a signal of DM in this parameter space (or to constrain these parameters) is the primary goal of the DM searches at the LHC Run II [5], and the most important signature in this model is mono-jet plus MET. The di-jet final state via the Drell–Yan process can be an important complementary channel.

Spin-0 mediator model

Similarly, in the case of a spin-0 mediator () interacting with the Dirac fermion DM and the SM particles, we define the interaction part of the Lagrangians asandwhere are the scalar/pseudo-scalar couplings of DM and quarks. Assuming a UV complete description of the scalar theory with the couplings of the mediator to the SM particles proportional to the particle masses, we normalise these couplings to the SM Yukawa couplings, , and set all flavour off-diagonal couplings to zero. This implies that, in a five-flavour scheme with massless bottom quarks, only top quarks are relevant for DM production in this model. Extension to a four-flavour scheme with massive bottom quarks is possible. The model file for the spin-0 mediator case, including other choices for the spin of the DM particle (real scalar and complex scalar ), is also available at the FeynRules repository [28].

The pure scalar and pure pseudo-scalar mediator scenarios are given by setting the parameters in the Lagrangians (10) and (11) toandrespectively. All the other and parameters are irrelevant. Similar to the spin-1 case, the model has only four independent parameters as in (9).

In the spin-0 mediator model with Yukawa-type couplings, the most relevant tree-level process at the LHC is DM pair production associated with a top-quark pair. On the other hand similarly to Higgs production, at one loop, gluon fusion can give rise to MET + jets signatures which are in general phenomenologically relevant. For the heavy mediator case, the four-top final state can be also relevant.

At this stage, we do not see compelling reasons to introduce couplings of the mediator to leptons, even though it is straightforward to do. We do not include effective gluon–gluon– interactions either, for several reasons. The first is that this operator is higher dimensional () and therefore might lead to unitarity-violating effects that need to be studied on a model and benchmark basis. The second is that a simplified model assumes no other new physics particle beyond the dark matter particle X and mediator Y at the weak scale. If such new particle effects decouple with their mass, the main contribution to the gluon–gluon– coupling would then be due to loop of SM quarks. Depending on the masses involved (that of the DM, the mediator and the quarks contributing in the loop) and the momentum transfer of the process, such interactions might be considered point-like. This is, however, very much model and process dependent. To be accurate one should first calculate the loop-induced processes exactly and study the range of applicability of such an effective interaction by explicit comparison. This is possible in MG5aMC [43] and has been considered with the same DM implementation as presented here in Ref. [44]. Other studies of the loop-induced process for mono-jet + MET can be found in Refs. [45-48]. We also note that couplings of the mediator to the SM gauge bosons can be introduced easily [42].

Dark matter production with jets

In this section, we consider a spin-1 mediator scenario and discuss the impact of the NLO-QCD corrections on DM pair production with jets, i.e.,In MG5aMC the code and events for the above process can be automatically generated by issuing the following commands:

We have checked that our model can reproduce the SM predictions for by adjusting the corresponding coupling and mass parameters.

To illustrate the effect of the higher-order corrections, we consider pure vector, Eqs. (5) and (6), or pure axial-vector, Eqs. (7) and (8), couplings with a simplified flavour structure. We takeas our benchmark for the spin-1 mediator scenario.

We assume that the mediator can only decay into the DM particle and the SM quarks (if kinematically allowed) through the interactions specified in Eqs. (3) and (4), so that the mediator width, , is determined by the particle masses and the couplings and . In our framework, the width is automatically computed by using the MadWidth module [49] for each parameter point. The above benchmark coupling strength in (17) leads to for and for , both for the vector and axial-vector cases. Note that, if we take with , the width becomes very large as for .

We provide LO and NLO QCD predictions for at the center-of-mass energy  TeV. The central value for the renormalisation () and factorisation () scales is set to , where is the sum of the transverse momenta of all jets in the event and the missing transverse energy. The scale uncertainty is estimated by varying the scales and , independently, by a factor two around . We adopt the five-flavour scheme and the LO and NLO NNPDF2.3 set [50] through the LHAPDF interface [51], with the corresponding values of and , for the LO and NLO predictions, respectively. The PDF uncertainties are computed automatically [52], following the prescription summarised in [53].

Where relevant, we apply a parton shower to the events using Pythia8 [35]. We then define jets using the anti- algorithm [54] as implemented in FastJet [55] with the jet cone radius , where we require  GeV and for all jets in the event.

Total cross sections

In Table 1 we present LO and NLO cross sections (in pb) for DM pair production in association with a jet at fixed order (FO) in perturbation theory. We show results for a pure vector mediator by fixing the parameters as in Eqs. (5), (6) and (17). We cover various benchmark points suggested by the ATLAS/CMS DM forum [5] in the – plane, representing four different cases: on-shell () and off-shell () production of the mediator, in the threshold regime () and in the EFT limit (). We also present scale and PDF uncertainties in % as well as K factors which we define as the ratio of the central values of the NLO and LO cross sections. We compute the table entries with different MET cuts: 150, 300, and 500 GeV. For convenience, we also show a graphical summary of our results in Fig. 1. As a reference, the cross sections for are also shown, where the vector mediator is produced on-shell and does not decay. For , the mono-jet rate is given by in the narrow width approximation.

Table 1

LO and NLO cross sections and corresponding K factors for DM pair production in association with a jet for the vector mediator scenario at the 13-TeV LHC, where different MET cuts are imposed. The uncertainties represent the scale and PDF uncertainties in per cent, respectively. We show several benchmark model points for the mediator and DM masses with the coupling parameters and

(m_Y,m_X) [GeV]			Vector
			MET >150 GeV	MET >300 GeV	MET >500 GeV
		σ_LO [pb]	2.923×10^2-8.9+10.7 ±1.6%	1.734×10^1-11.9+14.2 ±1.1%	1.695×10^0-14.0+17.4 ±1.8%
10	Undecayed	σ_NLO [pb]	5.093×10^2-8.2+10.3 ±0.5%	2.689×10^1-9.1+10.4 ±0.6%	2.433×10^0-10.0+11.1 ±1.1%
		K factor	1.74	1.55	1.44
		σ_LO [pb]	1.605×10^2-8.9+10.7 ±1.6%	0.978×10^1-12.0+14.3 ±1.1%	0.970×10^0-14.1+17.4 ±2.0%
(10, 1)	m_Y>2m_X	σ_NLO [pb]	2.818×10^2-8.1+10.1 ±0.5%	1.517×10^1-8.9+10.0 ±0.6%	1.345×10^0-9.6+10.5 ±1.1%
		K factor	1.76	1.55	1.39
		σ_LO [pb]	2.434×10^0-10.1+11.8 ±1.5%	2.843×10^-1-12.5+15.0 ±1.2%	3.786×10^-2-14.5+18.0 ±2.4%
(10, 50)	m_Y<2m_X	σ_NLO [pb]	3.198×10^0-5.4+5.6 ±0.5%	3.485×10^-1-6.3+5.9 ±0.7%	4.325×10^-2-7.8+7.3 ±1.3%
		K factor	1.31	1.23	1.14
		σ_LO [pb]	6.968×10^-3-14.0+17.4 ±4.3%	2.314×10^-3-15.0+18.9 ±4.6%	7.317×10^-4-16.1+20.6 ±5.6%
(10, 500)	m_Y<2m_X	σ_NLO [pb]	7.698×10^-3-6.4+5.4 ±2.2%	2.385×10^-3-6.9+5.7 ±2.3%	6.800×10^-4-7.1+5.5 ±2.6%
		K factor	1.10	1.03	0.93
		σ_LO [pb]	2.148×10^2-9.3+10.6 ±1.5%	1.616×10^1-12.0+14.4 ±1.0%	1.644×10^0-14.1+17.4 ±1.9%
100	Undecayed	σ_NLO [pb]	3.011×10^2-5.9+6.6 ±0.5%	2.121×10^1-7.1+7.3 ±0.6%	1.955×10^0-8.2+8.1 ±1.2%
		K factor	1.40	1.31	1.19
		σ_LO [pb]	1.100×10^2-9.3+10.6 ±1.5%	0.822×10^1-12.0+14.4 ±1.1%	0.862×10^0-14.1+17.4 ±1.9%
(100, 1)	m_Y>2m_X	σ_NLO [pb]	1.530×10^2-5.7+6.5 ±0.5%	1.100×10^1-7.2+7.4 ±0.6%	1.059×10^0-8.1+8.0 ±1.2%
		K factor	1.39	1.34	1.23
		σ_LO [pb]	1.117×10^1-9.6+11.0 ±1.5%	0.988×10^0-12.2+14.7 ±1.1%	1.140×10^-1-14.2+17.6 ±2.0%
(95, 50)	m_Y≲2m_X	σ_NLO [pb]	1.512×10^1-5.5+6.0 ±0.5%	1.281×10^0-6.8+6.8 ±0.6%	1.325×10^-1-7.6+7.2 ±1.2%
		K factor	1.35	1.30	1.16
		σ_LO [pb]	7.043×10^-3-14.0+17.4 ±4.3%	2.329×10^-3-15.0+18.9 ±4.6%	7.395×10^-4-16.1+20.6 ±5.6%
(100, 500)	m_Y<2m_X	σ_NLO [pb]	7.804×10^-3-6.4+5.3 ±2.2%	2.411×10^-3-6.8+5.5 ±2.3%	6.908×10^-4-7.1+5.5 ±2.6%
		K factor	1.11	1.04	0.93
		σ_LO [pb]	2.248×10^0-13.2+16.1 ±3.2%	6.865×10^-1-14.3+17.7 ±3.3%	1.979×10^-1-15.5+19.6 ±4.1%
1000	Undecayed	σ_NLO [pb]	2.601×10^0-6.0+5.1 ±1.7%	7.393×10^-1-6.4+5.2 ±1.8%	1.909×10^-1-6.8+5.3 ±2.1%
		K factor	1.16	1.08	0.96
		σ_LO [pb]	1.093×10^0-13.3+16.4 ±3.1%	3.278×10^-1-14.4+18.0 ±3.3%	9.182×10^-2-15.6+19.7 ±4.1%
(1000, 1)	m_Y>2m_X	σ_NLO [pb]	1.215×10^0-5.5+4.2 ±1.7%	3.399×10^-1-6.0+4.5 ±1.7%	8.743×10^-2-6.5+4.8 ±2.0%
		K factor	1.11	1.04	0.95
		σ_LO [pb]	1.094×10^0-13.3+16.4 ±3.1%	3.268×10^-1-14.4+18.0 ±3.3%	9.137×10^-2-15.6+19.7 ±4.1%
(1000, 50)	m_Y>2m_X	σ_NLO [pb]	1.221×10^0-5.6+4.3 ±1.7%	3.416×10^-1-6.0+4.6 ±1.7%	8.807×10^-2-6.6+4.9 ±2.0%
		K factor	1.12	1.05	0.96
		σ_LO [pb]	2.169×10^-1-13.3+16.4 ±3.4%	6.777×10^-2-14.4+18.0 ±3.6%	1.981×10^-2-15.6+19.7 ±4.4%
(995, 500)	m_Y≲2m_X	σ_NLO [pb]	2.497×10^-1-6.2+5.3 ±1.8%	7.223×10^-2-6.6+5.5 ±1.9%	1.914×10^-2-6.8+5.3 ±2.1%
		K factor	1.15	1.07	0.97
		σ_LO [pb]	8.487×10^-6-14.3+18.0 ±4.3%	2.666×10^-6-15.7+20.0 ±5.5%	8.238×10^-7-17.0+22.0 ±7.3%
(10000, 1)	m_Y≫s^	σ_NLO [pb]	8.835×10^-6-5.1+3.1 ±2.5%	2.579×10^-6-5.5+3.1 ±3%	7.148×10^-7-7.0+5.0 ±4.4%
		K factor	1.04	0.97	0.87

Fig. 1

Summary plot of NLO cross sections and corresponding K factors in Table 1

The production rate strongly depends on the both masses as well as on the kinematic cuts, and varies by orders of magnitude in the parameter scan. On the other hand, the K factors, i.e. higher-order effects, are not so sensitive to the mass spectra; for the heavy-mediator and/or heavy-DM cases, while –1.4 for the mediator with light DM, assuming the  GeV cut. We find that in the case of a relatively light mediator with a very light DM,  GeV, the K factor can reach a value as large as 1.8.

Different benchmark points probe different Bjorken-x regions of the parton distribution functions. As heavy mediators/DM are produced from very high-x partons, the dominant contribution comes from the initial state, as the gluon PDF is sub-dominant in the high-x region. For light mediators with light DM, on the other hand, a large contribution arises from the qg initial state. For instance, we find that the ratio of production cross sections via and ug initial states, , is 1.4 in the case of while 0.2 for , at LO.

As expected, most of the results at NLO accuracy display significantly smaller scale uncertainties compared to the LO calculations. An exception is provided by the  GeV case, which we discuss in detail at the end of this subsection. The PDF uncertainties are sub-leading in both the LO and NLO results and reduced by going from LO to NLO. Furthermore, the scale and PDF uncertainties increase when the mass scale of the mediator and/or DM increases.

The higher cut leads to smaller K factors and to larger scale and PDF uncertainties, which one can clearly see in the MET distributions in Fig. 2.

Fig. 2

MET distributions at FO (N)LO accuracy for at the 13-TeV LHC for and (1000, 50) GeV, where we assume a pure vector mediator and Dirac DM. The middle and bottom panels show the differential scale uncertainties and K factors, respectively

In Table 2, we present the pure axial-vector mediator case by fixing the parameters as in Eqs. (7), (8) and (17). The resulting cross sections are very similar compared to the pure vector case for , while in the off-shell regime, we find that the cross sections are suppressed compared to the production via pure vector mediators. In the off-shell situation the DM pair is produced at threshold. A pair of (Dirac) DM originating from a decay of a spin-1 mediator will be in a state with . If the coupling is vector-like the DM pair can be in a state, while if it is axial-like it will be in a state, i.e. suppressed at threshold. The similar argument holds in case of case, as we show in Sect. 4.

Table 2

Same as Table 1, but for the axial-vector mediator scenario

(m_Y,m_X) [GeV]			Axial-vector
			MET >150 GeV	MET >300 GeV	MET >500 GeV
		σ_LO [pb]	2.130×10^2-9.3+10.6 ±1.6%	1.573×10^1-12.0+14.4 ±1.1%	1.633×10^0-14.0+17.3 ±1.9%
100	Undecayed	σ_NLO [pb]	3.063×10^2-6.1+6.9 ±0.5%	2.153×10^1-7.4+7.7 ±0.6%	2.055×10^0-8.3+8.4 ±1.6%
		K factor	1.44	1.37	1.26
		σ_LO [pb]	1.101×10^2-9.3+10.6 ±1.6%	0.825×10^1-12.1+14.4 ±1.1%	0.854×10^0-14.1+17.4 ±2%
(100, 1)	m_Y>2m_X	σ_NLO [pb]	1.549×10^2-6.0+6.8 ±0.5%	1.127×10^1-7.2+7.4 ±0.6%	1.063×10^0-8.2+8.2 ±1.2%
		K factor	1.41	1.37	1.24
		σ_LO [pb]	3.070×10^0-10.0+11.6 ±1.5%	3.359×10^-1-12.4+14.9 ±1.2%	4.457×10^-2-14.3+17.7 ±1.8%
(95, 50)	m_Y≲2m_X	σ_NLO [pb]	4.093×10^0-5.7+6.0 ±0.5%	4.302×10^-1-6.9+6.7 ±0.7%	5.079×10^-2-7.4+6.9 ±1.3%
		K factor	1.33	1.28	1.14
		σ_LO [pb]	2.298×10^-3-14.5+18.1 ±5%	7.839×10^-4-15.4+19.5 ±5.3%	2.558×10^-4-16.5+21.2 ±6.3%
(100, 500)	m_Y<2m_X	σ_NLO [pb]	2.502×10^-3-6.8+5.9 ±2.5%	7.972×10^-4-7.3+6.2 ±2.6%	2.383×10^-4-7.5+6.1 ±3.0%
		K factor	1.09	1.02	0.93

The NLO effects are very similar to the vector mediator scenario for all the mass combinations as well as the cuts. Although we do not show the mixed scenario of vector and axial-vector, one can easily compute such scenarios by changing the coupling parameters in our simplified model.

The parameter point  GeV warrants special attention, as it illustrates a case of large NLO corrections (so-called “giant K factors” [56]), which might arise in the limit where . The giant K factors in the process occur due to the opening of the channel at NLO. This process can lead to a di-jet event topology with a soft, collinear emission of . In the regime of , the emission behaves similar to an emission of a massless gauge boson, where the diagrams with dijet topologies contribute factors of , and hence NLO K factors which scale as . Similar features commonly appear in calculations of SM W / Z+jets processes at high jet  [56].

Topologies leading to giant K factors are naturally suppressed in the case of DM production by the cut on MET. This restricts the calculation to regions of phase space which are insensitive to the soft and collinear double-logs of di-jet event topologies with a soft emission. Figure 3 (left) illustrates the effect in case of  GeV and  GeV. The region of low missing energy displays a two to three orders of magnitude difference in rate between the LO and NLO calculations, whereas we see that above , the K factor is drastically reduced. On the other hand, we see that already for  GeV in Fig. 3 (right), such logarithmic enhancements are only very weak.

Fig. 3

MET distribution at (N)LO+PS accuracy for at the 13-TeV LHC for (left) and (right). The lower panels provide information on the differential scale uncertainty and K factor

Differential distributions

We proceed with the discussion of the features of the differential distributions relevant for DM studies. We begin with Fig. 4 which shows the distributions at LO and NLO for four benchmark points of the simplified model, assuming a pure vector mediator and Dirac fermion DM. As seen in the total rates, the NLO effects in the distributions do not depend on the mass relation between the mediator and the DM, i.e. on-shell or off-shell, but do depend on the energy scale of the final state. In the top panels, the energy scale is  GeV for or . We find that the two benchmark points display striking similarities in the shape of the distributions, while the rate of the latter is suppressed due to off-shell production. The largest effects of NLO corrections are in the low regions, where NLO corrections reach K factors of about 1.4 for , with a steady decrease with increasing . We observe similar features also in the high-scale benchmark points of  TeV for or (bottom panels of Fig. 4), where the largest K factors are about 1.2 for . Comparing with the FO distributions in Fig. 2, we observe that the parton shower does not affect the MET distribution. Note that the NLO corrections are different for different regions, with the largest NLO corrections occurring in the lower MET regions where the rate is the highest. Hence the careful estimation of NLO effects is very important for accurate LHC studies of DM in each signal region.

Fig. 4

MET distributions for at the 13-TeV LHC for four benchmark points specified by (), where we assume a pure vector mediator and Dirac DM. The middle and bottom panels show the differential scale uncertainties and K factors, respectively

Fig. 5

distributions of the hardest jet for at the 13-TeV LHC for four benchmark points specified by (), where we assume a pure vector mediator and Dirac DM and the cut is imposed. The middle and bottom panels show the differential scale uncertainties and K factors, respectively

Fig. 7

distributions of the hardest jet for at the 13-TeV LHC for four benchmark points specified by (), where we assume a pure vector mediator and Dirac DM and the cut is imposed. The middle and bottom panels show the differential scale uncertainties and K factors, respectively

Next, we study the features of jet kinematic distributions produced in association with DM. Figures 5, 6, 7 and 8 show example and distributions of the hardest and second hardest jets for the four benchmark points as in Fig. 4, assumingDistributions of the hardest jet transverse momentum show very interesting features. In Fig. 5 we find that, in all benchmark points, the LO distributions match the NLO predictions at the peak, i.e.  GeV, to a very good degree. The agreement can be attributed to the imposed cut in (18), which forces the events into a back-to-back configuration of the leading jet and the mediator (on average). We also note that the NLO scale uncertainty in the peak region becomes very small compared to the LO estimates.

Fig. 6

Same as Fig. 5, but for the second hardest jet

Fig. 8

Same as Fig. 7, but for the second hardest jet

The NLO corrections to distributions affect not only the overall rate, but the shape of the distribution as well. In the low- region, K factors are about 1.2–1.5. In the high region, we find significant NLO effects again for the and (95, 50) GeV cases (top panels), but not for the and (1000, 50) GeV cases (bottom panels). We note that the scale uncertainty does not significantly reduce at NLO in the regions away from the peak, especially for the light mediator and DM case (top panels). Significant differences in NLO contributions and theoretical uncertainties in different regions of the spectrum suggest that the proper modelling of the hardest jet differential distributions has to go beyond the simple scaling by a constant K factor.

Apart from the highest jet which is modelled by the hard matrix element, all other jets in the LO simulation are generated by the parton shower. By contrast, the NLO corrections include real emission diagrams which can contain two hard and well-separated partons in the final state as well as virtual corrections to one parton emission. One could expect significant differences between LO and NLO in the kinematic distributions of the second highest jet. For the and (95, 50) GeV cases (top panels), we observe giant K factors in the high- tails of the distributions. The large difference between LO and NLO computations is a consequence of the inadequacy of the parton shower to accurately model high- emissions. In Fig. 6 (bottom panels), on the other hand, we find no significant differences between LO and NLO for the overall rate and shape of the second jet emission in case of very heavy mediators (i.e. ) or heavy DM (i.e. ), suggesting that the second hardest jet is described very well by the parton shower. This is because the scale of the shower is very high and therefore extra parton emission from the parton shower can be sufficiently hard.

Features similar to those observed in also occur in distributions of the hardest/second-hardest jet pseudo-rapidity (), shown in Figs. 7 and 8. For the light mediator/DM (top panels), we observe that the rate at which the hardest jet is emitted at NLO in the low rapidity region is enhanced by a factor about 1.5, with the corrections falling off with the increase in rapidity. However, even though the overall rate for the second-hardest jet increases by a factor of roughly the shape of the distribution is affected only mildly. In the case of heavy mediator/DM, the hardest jet is emitted at a lower rapidity (on average) at a significantly higher rate compared to light mediators as illustrated by the width of the distributions in Fig. 7. As the hardest jet typically recoils against , this explains why the spectrum falls off more quickly for lighter mediators than for the heavier ones.

Merging samples at NLO accuracy

In addition to total and differential production cross sections for the process, we study NLO effects for different jet multiplicities in the final state. For this purpose we utilise the FxFx merging procedure [20] within the framework of MG5aMC, and considerWe take the merging scales at 45 and 30 GeV for 2- and 1-jet merged samples, respectively.

Figure 9 shows the number of jets in the final state for NLO merged samples in case of  TeV and  GeV. The red and blue curves show the results of merging up to 1 and 2 jets, respectively, while the orange curve shows the process at NLO+PS without merging for comparison. An inspection of the three samples in the lowest panel of Fig. 9 shows that the effects of NLO merging are mild. The non-merged NLO sample over-estimates the production rate in the 0j and bins by 20 and 10 % respectively, and underestimates the rate in the 1j–2j bins by . The differences are even milder between the samples merged to 1j and to 2j. As the 0j bin is phenomenologically irrelevant, we can conclude that the effects of jet merging at NLO are within .

Fig. 9

Jet multiplicity for jets at the 13-TeV LHC with FxFx merging in case of  TeV with  GeV, where NLO+PS samples are merged up to one (red) and two (blue) jets. The 1-jet NLO+PS sample without merging (orange) is also shown for comparison. The middle panel shows the relative scale uncertainties, while the bottom panel presents the ratio of the 2j merged sample to the 1j merged one and to the non-merged one. We assume and jets with and for the purpose of illustration

We show effects of jet merging on the MET distribution in Fig. 10. Except in the low MET region, we find that the effects of NLO merging are again mild and within 10 %.

Fig. 10

Same as Fig. 9, but for the MET distribution

Comparison of signal distributions to the Standard Model

In discussions of NLO corrections to DM production, it is important to consider how the possible signal events at NLO look in the midst of large SM backgrounds. For the purpose of illustration, we consider only the largest background in mono-jet searches for DM, i.e. Z+jets and simulate it to NLO merged up to 2 extra jets via the FxFx method.

Figure 11 shows an example comparison of the Z+jets channel to several benchmark points of the simplified model discussed in previous sections. The shape of signal jet multiplicity distributions (upper left panel) resembles the Z+jets distribution to a good degree, while the overall rate varies wildly depending on the model point. Note that events containing one jet are produced at almost an identical rate to 2j events and comparable to 3j events, both in Z+jets and all of the benchmark model points we considered. The production rate for different jet multiplicities implies that it could be beneficial to consider DM searches beyond mono-jet, either inclusive, or at fixed jet multiplicities.

Fig. 11

Distributions for various signal benchmark points and jets, where and jets with and are considered. We assume a pure vector mediator and Dirac DM with the coupling parameters and

Next, the MET (as well as the hardest jet ) distributions for mediators of mass with light DM naturally resemble the Z+jets distribution in shape, as the kinematics of the MET and hardest jet are determined by the mass scale of the heavy object (i.e. the Z boson or the mediator). For heavy mediators, the MET and distributions fall off with a milder slope, suggesting that the signal to background ratio (S / B) in searches for DM could be improved by requiring a higher for the hardest jet.

The distribution of the second hardest jet, on the other hand, seems to display a similar shape in all model points, as well as the Z+jets background channel. Requiring a high on the second jet would hence not improve neither S / B nor the signal significance, suggesting that in searches which exploit the presence of a second jet, only minimum cuts on should be applied. Note that the modelling of the second jet from LO calculations grossly underestimates the overall rate in case of lighter mediators for , as illustrated in Fig. 6.

Dark matter production with a top-quark pair

In the spin-0 mediator model, due to the normalisation of the Yukawa couplings in the Lagrangian (11), the most relevant tree-level process at the LHC isSuch models have in the past been studied in the context of EFT [57-59] and simplified models [46, 48], and searched for at the LHC Run I [60, 61]. Past work on DM interactions with the top quarks has mainly focused on LO estimates, with only a few analyses including NLO corrections. Here we present a comprehensive study of NLO effects of DM interaction with top quarks in the framework of the simplified model. We note that a wide class of so called “top-philic DM” models exist where the LO production is via top loops. Reference [62] studied such a scenario in a minimally model-dependent framework, while more recently, Ref. [44] presented concrete predictions for loop-induced DM production for the current LHC13 run, using the same simulation framework as in this work.

The code and events for the above process can be automatically generated by issuing the following commands in MG5aMC:

We have checked that our model can reproduce the SM predictions for by adjusting the coupling and mass parameters. Note that we use the on-shell renormalisation for the NLO model construction. The top-quark decays can be subsequently performed by MadSpin [63], which keeps production and decay spin correlations.

To illustrate the NLO effects, we consider pure scalar, Eqs. (12) and (13), or pure pseudo-scalar, Eqs. (14) and (15), couplings, and takeas the default couplings. With these values, the scalar and pseudo-scalar mediator width is –0.1 for , while for and . The width for scalar is slightly smaller due to the additional factor where .

We start by showing total production rates of , where the top quark is considered stable. Table 3 shows the LO and NLO cross sections (in pb) for the scalar and pseudo-scalar mediator scenarios, where we use  GeV. The central renormalisation and factorisation scales are set to half the sum of the transverse mass of the top quarks and the missing transverse energy. We also present scale and PDF uncertainties in % as well as K factors.

Table 3

LO and NLO cross sections and corresponding K factors for DM pair production in association with a top-quark pair for the scalar and pseudo-scalar mediator scenario at the 13-TeV LHC. The uncertainties represent the scale and PDF uncertainties in per cent, respectively. We show several benchmark model points for the mediator and DM masses with the coupling parameters and

(m_Y,m_X) [GeV]			Scalar	Pseudo-scalar
		σ_LO [pb]	2.278×10^1-20.4+28.0 ±4.2%	5.202×10^-1-22.0+30.8 ±6.0%
10	Undecayed	σ_NLO [pb]	2.435×10^1-8.5+5.4 ±1.8%	5.431×10^-1-10.2+7.4 ±2.6%
		K factor	1.07	1.04
		σ_LO [pb]	2.294×10^1-20.5+28.0 ±4.2%	5.500×10^-1-22.1+30.8 ±6.0%
(10, 1)	m_Y>2m_X	σ_NLO [pb]	2.460×10^1-8.5+5.4 ±1.8%	5.739×10^-1-10.2+7.4 ±2.6%
		K factor	1.07	1.04
		σ_LO [pb]	2.415×10^-3-21.8+30.5 ±5.8%	3.329×10^-3-23.8+33.9 ±8.7%
(10, 50)	m_Y<2m_X	σ_NLO [pb]	2.340×10^-3-9.1+5.8 ±2.8%	3.133×10^-3-11.0+7.5 ±3.9%
		K factor	0.97	0.94
		σ_LO [pb]	8.226×10^-1-20.9+28.7 ±4.4%	2.442×10^-1-22.9+32.2 ±7.2%
100	Undecayed	σ_NLO [pb]	8.391×10^-1-8.6+5.3 ±2.1%	2.431×10^-1-10.7+7.6 ±3.2%
		K factor	1.02	1.00
		σ_LO [pb]	8.135×10^-1-20.9+28.8 ±4.4%	2.464×10^-1-23.0+32.4 ±7.2%
(100, 1)	m_Y>2m_X	σ_NLO [pb]	8.207×10^-1-8.3+4.8 ±2.1%	2.427×10^-1-10.4+7.0 ±3.2%
		K factor	1.01	0.98
		σ_LO [pb]	7.986×10^-3-21.3+29.5 ±5.0%	1.404×10^-2-23.3+32.9 ±7.8%
(95, 50)	m_Y≲2m_X	σ_NLO [pb]	7.897×10^-3-8.8+5.5 ±2.4%	1.362×10^-2-10.8+7.4 ±3.5%
		K factor	0.99	0.97
		σ_LO [pb]	1.571×10^-3-27.0+40.2 ±17.1%	1.827×10^-3-27.1+40.4 ±17.4%
1000	Undecayed	σ_NLO [pb]	1.127×10^-3-13.6+10.9 ±8.7%	1.297×10^-3-13.7+10.8 ±8.8%
		K factor	0.72	0.71
		σ_LO [pb]	7.499×10^-4-27.2+40.8 ±16.8%	8.174×10^-4-27.3+40.9 ±17.0%
(1000, 1)	m_Y>2m_X	σ_NLO [pb]	5.201×10^-4-12.7+8.4 ±8.4%	5.675×10^-4-12.9+8.6 ±8.5%
		K factor	0.69	0.69
		σ_LO [pb]	7.354×10^-4-27.2+40.7 ±16.8%	8.137×10^-4-27.3+41.0 ±17.0%
(1000, 50)	m_Y>2m_X	σ_NLO [pb]	5.125×10^-4-12.8+8.6 ±8.5%	5.595×10^-4-12.7+8.3 ±8.5%
		K factor	0.69	0.69

Fig. 12

Distribution of the invariant mass of the top-quark pair for at the 13-TeV LHC for different mediator masses with the DM mass fixed at 50 GeV, where we assume a pure scalar mediator and Dirac DM. The middle and bottom panels show the differential scale uncertainties and K factors, respectively

For the total production rates, the NLO effects are very mild for the light mediator case, while they are significant for the heavy case. The inclusion of NLO corrections results in a drastic improvement of the scale uncertainties, from up to 40 % at LO to typically only about 10 % at NLO. Also, the PDF uncertainties are reduced by approximately a factor of two when going from LO to NLO.

Table 3 also shows clear differences between the overall production rates in the cases of scalar and pseudo-scalar mediators. For mediator mass of  GeV we find that DM production cross section via scalar mediators is an order of magnitude larger compared to the production rate via the pseudo-scalar mediator with the same mass. The large difference occurs due to the fact that in case of , the production cross section is dominated by the fragmentation. In case of the scalar mediator, the fragmentation function contains terms with soft singularities of the form – where x is the momentum fraction carried by the mediator – causing enhancements in the production rate [64]. The soft-enhanced term is absent in the case of a pseudo-scalar mediator [65], explaining the order of magnitude difference between the total rates of the scalar and pseudo-scalar mediators.

In cases where either DM or the mediator is produced close to threshold, we observe that the production cross section in the pseudo-scalar mediator case is larger. The effect can be attributed to the production rate originating mainly from top fusion diagrams. The production of a DM (Dirac) pair via scalar mediators at threshold can proceed only via a P-wave () and is hence suppressed by extra two powers of [66]. Conversely, production of DM pair via pseudo-scalar mediators can proceed via an S-wave () and hence does not suffer any kinematic suppression.

For the study of differential distributions, we consider the invariant mass of the top-quark pair (), without inclusion of a parton shower. Figures 12 and 13 illustrate the scalar and pseudo-scalar results, respectively, for different mediator masses (off-shell, threshold, and on-shell) with the DM mass fixed at 50 GeV. In all cases the shape of is well modelled by the LO calculation, and including a constant K factor reproduces the NLO results to an excellent degree, except in the threshold region. However, the LO calculation suffers from significant scale uncertainties which tend to increase with , whereas the scale uncertainties are under much better control at NLO.

Fig. 13

Same as Fig. 12, but for the pseudo-scalar mediator scenario

Whether DM is produced via scalar or pseudo-scalar mediators can have a dramatic effect on the shape of the distribution. In Fig. 14 we compare the NLO distributions in Figs. 12 and 13, where we normalise the histograms to unit area to point out the shape differences. We observe that the shape of the distribution is particularly enhanced for in the case of the pseudo-scalar mediator, while the scalar mediator distribution displays a much more prominent peak at lower . The effect is severely damped in case of heavy mediators, where we find no clear differences between the shapes of the scalar and pseudo-scalar mediator distributions. Figure 14 suggests that scalar and pseudo-scalar mediators could be distinguished based on the shape of the distribution, as long as the mediator is sufficiently light and/or does not decay highly off-shell. An analogous observation has been made already in the case of the study of the CP properties of the Higgs boson [67].

Fig. 14

Comparison of invariant mass distributions between the scalar (blue) and pseudo-scalar (red) mediator models for different mediator masses

Summary

Searches for DM are one of the main endeavours at the LHC Run II. Accurate and precise predictions for production rates and distributions are necessary to obtain robust constraints on DM models and characterise possible DM signals. In this article we have provided a general implementation of the simplified DM model approach into a calculation/simulation framework that allows to systematically evaluate and include NLO QCD corrections to the production of DM at the LHC. We have considered a class of simplified models where DM is a Dirac fermion and couples to the SM via either spin-1 or spin-0 s-channel mediators, making no restrictions on chiral couplings. For the purpose of illustration, we analysed the NLO effects on the DM production via vector and axial-vector mediators in the context of mono-jet signals. In addition, we have presented detailed predictions of DM production in association with a top-quark pair via scalar and pseudo-scalar mediators. We presented our results for various DM and mediator masses to cover benchmark points suggested by the ATLAS/CMS DM forum [5].

For MET+jets in the spin-1 mediator model, our results show that higher-order corrections have a significant effect both on the overall production rate as well as on the shape of relevant differential distributions, with a sizeable reduction of the scale and PDF uncertainties. The NLO corrections to the LO production rates can be large, with K factors of up to , and typically occur in parts of the model parameter space where the mass scale of DM and mediator is  GeV. For such scenarios, we also find large NLO effects on the shape of differential distributions in and the transverse momentum of the associated jets. Simplified models with heavy (e.g.  TeV) mediators/DM do not receive large NLO corrections, and we find that LO predictions describe both total production rates and shapes of differential distributions quite accurately. Distributions of the second hardest jet in the event are well modelled by the parton shower for heavy mediator/DM cases. On the other hand, for mediators/DM with masses of  GeV, the inclusion of NLO effects is essential for a proper description of and distributions, especially in the high- tails, where the NLO effects can be an order of magnitude.

So-called “giant K factors” can occur in NLO computations of DM production rates in the regions where . Such effects can be extremely large when considering mono-jet production rates, especially in phase-space regions with low . Imposing a sufficiently large cut and hence avoiding the soft/collinear singularities associated with the mediator emissions from high- jets efficiently mitigates the effect of giant K factors.

In our analysis we have gone beyond FO in perturbation theory and studied the effects of jet multiplicity merging at NLO accuracy. We found that FO calculations model the jet multiplicity and other differential distributions adequately well, with no significant effects on the shapes or overall rates coming from jet sample merging.

Comparisons with the NLO predictions for the leading SM background channel reveal that considerations of either inclusive or exclusive jet samples beyond one jet could be beneficial for increasing the prospects for DM detection. The leading jet distributions in case of heavy mediators display a milder decrease with the increase in , suggesting that a significant improvement in S / B could be obtained by focusing on high- regions. The second hardest jet distribution, on the other hand, was characterised by a slope similar to the leading backgrounds, implying that more inclusive cuts on the second jet should be used.

For MET+ in the spin-0 mediator model, our results show that the NLO corrections are very mild for the light mediator case, while they are significant for the heavy case. We observed a drastic improvement of the theory uncertainties when going from LO to NLO. We have noted that the shape of the distribution can reveal the chiral structure of the DM–SM interactions, as long as the mediator is relatively light (i.e.  GeV).

The DM model we studied in this paper is publicly available in the FeynRules repository [28]. We emphasise that all results presented here have been obtained in the FeynRules/MG5aMC framework, and thus they can be easily reproduced and used in DM searches at the LHC Run II.