Direct detection of dark matter: Precision predictions in a simplified model framework

We present a calculation of the next-to-leading order QCD corrections for the scattering of Dark Matter particles off nucleons in the framework of simplified models with s- and t-channel mediators. These results are matched to the Wilson coefficients and operators of an effective field theory that is generally used for the presentation of experimental results on spin-independent and spin-dependent direct detection rates. Detailed phenomenological studies illustrate the complementary reach of collider searches for Dark Matter and the direct detection experiments CRESST and XENON. In the case of cancellation effects in the tree-level contributions, one-loop corrections can have a particularly large impact on exclusion limits in the case of combined s+t-channel models.


Introduction
The However, apart from its gravitational properties very little is known about the nature of DM. In the context of particle physics, it is often assumed that DM consists of a new type of particle, not accounted for within the SM [3]. Prominent examples of such postulated particles are new types of neutrinos [4,5], axions -light particles introduced to resolve the strong CP problem of QCD [6,7] that could account for DM if their masses were in the meV range [8] -or weakly interacting massive particles (WIMPs) with a mass in the few-hundred GeV range (see, e.g., Ref. [9] for a recent review). Such WIMPs emerge, for instance, in supersymmetric models where the lightest stable supersymmetric particle is constituted by a neutralino. More recently there have also been attempts to account for WIMPs in a more generic way by the construction of so-called simplified models which aim at capturing the main features of these new particles and their interactions, in particular their mass, spin, and couplings to a specific mediator that in turn couples to SM particles (see, e.g., Ref. [10] for a recent discussion of several simplified models in the context of LHC searches). While not providing an ultimately UV-complete theory, such models have the advantage of being simple to use and featuring a relatively small number of parameters that are to be determined from experimental measurements. Once data have constrained the parameters of such a simplified model, theorists can use this input for the construction of more sophisticated extensions of the SM.
From the experimental side, this approach requires ways to constrain mass, spin, and coupling strengths of the DM candidate and the mediator particle. Ideally, such information is extracted from and cross-checked among conceptually entirely different types of searches. Indeed, one can distinguish three major types of experiments searching for Dark Matter: Indirect, astrophysical searches aiming to detect SM signatures resulting from DM annihilation processes; searches for DM production at high-energy colliders; and direct detection experiments that are designed to identify the recoil a DM particle causes in a nuclear target. Here, we will concentrate on the latter. Depending on the design of a specific direct detection experiment, its sensitivity is typically tailored to a particular mass range of the DM candidate: Detectors making use of dual-phase time projection chambers are particularly effective in the search for DM particles with a mass above about 10 GeV [11]. Cryogenic experiments are most sensitive to very light DM candidates, with a mass below a few GeV, because of their low energy threshold [11].
The recoil rate in direct detection experiments depends on several quantities of astrophysical origin, such as the local DM density and the velocity distribution of DM, as well as nuclear properties of the detector material and the nuclear scattering cross section.
For a multitude of simplified models, the leading expressions for the scattering cross sections can be obtained from the literature (see, e.g., Ref. [12]). In order to improve the theoretical predictions for the recoil rate, the uncertainties on all ingredients must be reduced. In particular, the calculation of the scattering cross section can be improved with conventional methods in particle physics, i.e. by including higher-order corrections in the perturbative expansion. For specific models, such higher-order terms as well as loop-induced contributions have been discussed, for instance in [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29].
In this work, we aim to provide accurate predictions for direct detection rates in the context of simplified models, including next-to-leading order (NLO) QCD corrections to the genuine tree-level processes where a DM particle scatters off a (anti-)quark from the nucleon. Including these perturbative corrections will allow us to obtain results with theoretical uncertainties related to unknown higher-order corrections reduced as compared to pure tree-level estimates. In order to assess remaining perturbative uncertainties in the considered observables we perform a detailed assessment of the dependence on unphysical scales. We show representative predictions within specific simplified models and discuss how exclusion limits of direct detection experiments can be used to derive constraints on model parameters such as the mass of the DM particle and mediator, or coupling strengths. We furthermore investigate the complementarity of results from direct detection experiments to those obtained at colliders, in particular recent measurements by the CRESST [30] and XENON [31,32] experiments, and at the LHC [33].
In Sec. 2 we describe the simplified model framework we are using as the basis of this study. We then explain how to connect these underlying models based on the properties of elementary particles to experimentally accessible direct detection rates in a nuclear medium. Numerical results are presented in Sec. 3, followed by our conclusions in Sec. 4.
Technical details of our calculation are summarized in the appendix.

Theory overview
In the following, we introduce the DM models that we are using in our analyses and describe the steps to obtain predictions for DM-nucleon cross sections, as relevant for direct detection experiments, up to NLO-QCD accuracy.

Simplified dark matter models
A particularly promising candidate for a DM particle is constituted by a -as of yet hypothetical -weakly interacting massive particle (WIMP) with a mass in the fewhundred GeV range that is subject only to the gravitational and weak interactions. Such a WIMP occurs in many extensions of the Standard Model, for instance models with additional Higgs doublets or in supersymmetric theories where the lightest neutralino provides a natural DM candidate.
Explicit calculations of experimentally accessible quantities (such as production rates and cross sections) can be performed within any such model. However, it can be advantageous to provide predictions of a more general nature that are not relying on the many parameters of a complex theory, such as supersymmetry, but only depend on those features that govern the interaction of the DM particle with the SM in the considered environment. To this end, two different strategies have been devised: The framework of effective field theories (EFTs), or so-called simplified models. Assuming that the interaction of the DM candidate with SM particles proceeds exclusively via a specific mediator particle, in the kinematic range where the mass of the mediator is much higher than all other energy scales of the considered reaction, in the EFT approach the heavy mediator particle is simply integrated out. Interactions between the SM and the DM sector can then be expressed in terms of effective operators. More control on details of a specific reaction is retained in simplified models. Such models are designed to capture the main features of reactions that are sensitive to the DM particle itself and the mediator between the SM and the DM sector. Contrary to the EFT approach, the dependence on the mediator particle is fully retained. Basic simplified models are constructed under the assumption that there is only a single type of DM particle, and a single type of mediator.
In this work, we consider simplified models with a single fermionic DM candidate.
The interaction between the SM and the DM fields proceeds either via the exchange of a neutral massive vector or axial-vector mediator, or via the exchange of a colorcharged scalar mediator. Because of the characteristic topologies of the dominant DM production modes associated with these models at colliders, the models with neutral or charged mediators are often referred to as "s-channel" or "t-channel" simplified models, respectively. We use the same notation for the terms in the Lagrangian as in our previous work [34].
For the simplified s-channel model we consider the interaction terms in the Lagrangian are of the form [35] where V denotes the vector mediator field with mass M V , χ the fermionic DM field with mass m χ , q a quark field with mass m q for each quark flavor q = u, d, s, c, b, t, and the coupling strength between the vector (axial-vector) mediator and fermionic For the t-channel models we assume a Lagrangian of the form [36][37][38]  are suppressed and we generically write "u" for all up-type quarks and "d" for downtype quarks. The λ Q L , λ u R , and λ d R are the Yukawa couplings of the mediator to the left-and right-handed quark fields, respectively. Due to the SU(2) L ×U(1) Y symmetry of the mediators, the left-handed quark coupling λ Q L is identical for up-and down-type flavors. Note that while in the s-channel model the mediator can be lighter than the DM particle, t-channel mediators are always heavier than the DM particles for them to remain stable.

Direct detection in the non-relativistic limit
The experimental signature of a DM particle of mass m χ scattering elastically off a nuclear target particle of type i and mass m i is given by the differential rate of nuclear recoil events dR per energy interval dE, where k i denotes the mass fraction of nuclear species i in the detector, σ i the DMnucleus scattering cross section, ρ 0 the local DM density, and µ i the reduced mass of the DM-nucleus system, The quantity η i depends on the distribution of the velocity v of the DM particles relative to the detector, where the integration range is determined by the minimum velocity v min,i of the DM particle required to cause a recoil energy E in the detector, and the galactic escape speed v esc beyond which DM particles are no longer gravitationally bound in the Milky Way.
The quantity f ( v) encodes the local velocity distribution in the detector rest frame.
DM particles move at non-relativistic velocities of about v ∼ 10 −3 [39]. Assuming the nuclei in the detector to be at rest, typical recoil energies, which are proportional to the momentum exchange q between the DM particles and the nuclei, are of the order of a few to a few hundred keV. We will therefore work in the limit of vanishing relative momentum between scattering DM particle and nucleus, where q 2 → 0, and refer to it as the non-relativistic limit.
In the differential detection rate of Eq. (3) all information on the microscopic DMnucleus interaction is contained in the elastic DM-nucleus cross section σ i , while the other quantities entering the equation are related to the detector composition and the DM relic density. In the non-relativistic limit, this elastic scattering cross section receives contributions only from so-called spin-independent (SI) and spin-dependent (SD) interactions. The cross section σ i therefore can be divided into an SI (σ SI i ) and an SD (σ SD i ) part. The strength of the contribution will differ for different types of nuclei, as the SI part is enhanced for heavy nuclei, while this enhancement is not present for the SD part (see, e.g., Ref. [40]). In order to allow for easier comparisons between different experimental setups using different detector materials, results are usually expressed in terms of the cross sections of the scattering between a single nucleon and a DM particle, keeping the division between the SI (σ SI N ) and SD (σ SD N ) parts. At low energies, elastic DM-nucleon scattering can be described by an effective fourfermion interaction which can be parameterized in terms of effective operators. We therefore start our discussion with an effective Lagrangian in relativistic notation with DM-quark interaction terms of the form with σ µν = i 2 [γ µ , γ ν ], and a priori unknown Wilson coefficients c j ≡ c j,q which in general depend on the quark flavor q. This Lagrangian contains the effective dimension-six operators O j accounting for scalar, vector, axial-vector, and tensor interactions of four fermions. These are the only operators that are not suppressed by small non-relativistic velocities. Mixed operators, such as e.g.χγ µ χqγ µ γ 5 q, are kinematically suppressed and will therefore not be discussed further. In the non-relativistic limit, the scalar and vector operators induce SI, the axial-vector and tensor operators SD interactions.
In order to derive from the EFT Lagrangian the scattering cross sections relevant for direct detection experiments, we can compute for each term in the Lagrangian the matrix element where for j = S, V, A, T the Γ j stands for 1, γ µ , γ µ γ 5 , i √ 2 σ µν , respectively. For computing quantities involving nucleons rather than elementary quarks, we adopt the conventional assumption that, in the non-relativistic limit, quark operators within nucleonic states are proportional to nucleonic operators. This effectively leads to a replacement of the quark fields by nucleon fields (see, e.g., Ref. [41]), with the nucleonic matrix elements f N j,q as proportionality factors, which encode the non-perturbative contributions from hadronic physics and are typically calculated in lattice-gauge theory or determined experimentally. The f N j,q depend on the type of the nucleon N = p, n, the interactions j, and the quark flavor q. In our calculations we use the numerical values for the f N j,q listed in Tab. 1. They have been obtained with the program micrOMEGAs [42,43], version 5.0.9. Note that values of the vector, axialvector, and tensor coefficients for neutrons can be obtained from the respective proton In order to capture all the quark content of the nucleon, a summation over all active quark flavors has to be performed. We therefore introduce the nucleonic couplings for each interaction type j. The total cross section for nucleon scattering in the nonrelativistic limit then simply reads: where the matrix element M j,N is obtained from M j,q by replacing the quark currents and Wilson coefficients with the corresponding nucleonic currents and couplings, and averaging (summing) over initial-state (final-state) polarizations and colors has been taken into account as indicated by the bar. In the non-relativistic limit, the center-ofmass energy squared of the DM-nucleon system can be written as s = (m χ + m N ) 2 , with m N being the mass of the considered nucleon.
Distinguishing between the contributions to SI and SD scattering and evaluating the fermionic traces within the squared amplitude, we obtain with the reduced mass µ N of the DM-nucleon particle system. The positive (negative) sign between the contributions from the different Wilson coefficients corresponds to DM (anti-DM) scattering.
In this work, we consider DM scattering in the framework of the simplified models introduced in Sec. 2.1. However, experimental results for direct detection experiments are often quoted in terms of EFT quantities. We therefore need to establish a strategy for a translation between the two approaches. We closely follow the procedure of Ref. [21] where a matching between neutralino-parton scattering in the framework of the MSSM and respective expressions in an EFT approach was presented. Instead of the SUSY Lagrangian of that reference we use the Lagrangian interaction terms of Eqs. (1)-(2) for the computation of scattering amplitudes in the framework of simplified models. To find the expressions of the coefficients c j,q in terms of the physical parameters of our models, we match, at the quark level, the simplified model amplitudes to the respective EFT expressions of Eq. (7), by imposing, at each order in perturbation theory, where M sim and M EFT denote the amplitude for DM-quark scattering in the simplified model and the EFT, respectively. We note that M EFT is a function of the Wilson coefficients c j,q , while M sim depends on the parameters of the simplified model. For specific settings of the simplified model parameters then a comparison with experimental limits on the Wilson coefficients can be performed, see, e.g., Ref. [44] for the CRESST-II experiment.

Radiative corrections in the non-relativistic limit
The leading-order (LO) Wilson coefficients contributing to SI and SD scattering are determined by matching the tree-level amplitudes of our simplified models, shown in Fig. 1 for the s-and t-channel models, respectively, to the lowest-order operator expression in the EFT, see purely on the model-specific couplings between the SM and DM sectors.
In order to improve the tree-level predictions for the Wilson coefficients, all the oneloop diagrams in the considered simplified models that contribute to the relevant operators at O(α s ) have to be calculated. In Fig. 3 representative one-loop diagrams in the s-and t-channel models are depicted. They can be classified as propagator, vertex, and box corrections, and have to be considered together with the corresponding counterterm diagrams. We denote these additional contributions as NLO-QCD or simply NLO corrections, as in this work we consider a perturbative expansion only in the QCD coupling α s .
Loop-induced processes with gluons in initial and final states contributing to DMnucleon scattering have been calculated in the framework of the minimally supersymmetric extension of the SM and similar models in [13][14][15]19,20]. We want to clarify that formally these contributions contribute only at O(α 2 s ) to the cross sections, which is one order higher than what we are discussing in our work, and we are thus not taking them into account.
Let us remark that while in general a full NLO-QCD calculation involves both virtual and real-emission corrections, in the non-relativistic limit real-emission contributions do not have to be considered. In this limit the in-and outgoing DM and quark momenta become equal, and no extra parton emission occurs. Consequently, the sum of all virtual corrections is infrared finite, which we have verified explicitly.
In order to convert the NLO corrections in the context of the simplified models we consider to one-loop expressions for the Wilson coefficients of the EFT, the matching condition of Eq. (12) has to be applied at NLO. Additionally, the one-loop operator corrections sketched in Fig. 2 (b) have to be considered. Spelled out in terms of effective operators, the NLO matching condition reads In the latter we do not compute terms like c 1-loop , as they contribute at O(α 2 s ), which is beyond the perturbative accuracy we are interested in. In the simplified s-and t-channel models that we consider, c tree S = c tree T = 0 (see App. A for details). At tree level we therefore obtain the following matching condition: At O(α s ), after inserting Eq. (14), we find:

Numerical analysis
We now present numerical results for the SI and SD DM-nucleon scattering cross sections, and discuss the effects of the O(α s ) corrections. Furthermore, we assess the potential of various direct detection experiments to explore regions in their respective parameter spaces when confronted with limits of LHC searches for DM production. As the numerical difference between the proton and neutron contribution to the SI cross section is marginal and the direct detection experiments typically publish limits irrespective of the nature of the nucleon, we limit ourselves to showing only the DM-proton cross sections for SI interactions. On the other hand, since the sensitivities of the experiments for the SD interaction differ between proton and neutron scattering, and the CRESST collaboration has only published limits for neutron scattering, we will focus on the DM-neutron cross sections in the case of SD interactions.
In the following, we do not discuss the pure s-channel model with a vector mediator further, since the O(α s ) corrections vanish for such a scenario, as discussed in App. A.1.
The impact of s-channel contributions is only considered in the context of interference effects in models that feature s-and t-channel topologies at the same time (see also App. A.3). For the t-channel model, we assume the DM particle to be a Dirac fermion.
We choose a scenario where all right-handed couplings are turned off, i.e. λ u R = λ d R = 0 for all quark flavors, and the mediators thus only couple to the left-handed SU(2) L quark doublets. This choice is motivated by the most recent mono-jet searches of the ATLAS experiment [33], where limits are presented for the same selection of couplings.
In the following, we use the short-hand notation λ ≡ λ Q L , and assume the mediators to couple to all left-handed quark doublets with equal strength. Unless explicitly specified otherwise, we choose λ = 1.
All quarks in our calculations are assumed to be massive and we use for the quark masses the current values compiled by the particle data group [45]. However, as in our with all dependence on α s explicitly factored out of the expansion coefficients σ (i) . Varying µ R thus affects both, the value of α s , and the size of the NLO correction term σ (1) which depends explicitly on renormalization scale logarithms. In the following discussion of our results, we will study the effects of varying the renormalization scale.
The above discussion does not yet include any renormalization-group running of the effective operators and, consequently, of the Wilson coefficients. If such running effects are taken into account, they will lead to an additional source of renormalization scale dependence, affecting both the σ (0) and σ (1) terms. We take into account the running effects of the operators appearing in our models following Ref. [46]. In practice, we find that the running is numerically only relevant for the axial-vector operator, though: While there are no running effects for the vector operator, the contributions of the scalar and tensor Wilson coefficients are too small for any running effects to be significant.
Using the specified setup, we first discuss the impact of NLO-QCD corrections and scale uncertainties on SI and SD scattering cross sections. After this assessment of the genuine features of experimentally accessible observables we turn to a systematic comparison of limits on DM models from colliders and direct detection experiments in the context of t-channel models. Subsequently, we investigate the impact of interference effects in combined s + t-channel models.

Theoretical aspects of the perturbative prediction
Before turning to a phenomenological discussion of the considered DM models, we explore the genuine features and theoretical uncertainties of our perturbative calculation.
In Fig. 4 we show the SI and SD DM-nucleon scattering cross sections for the t- in the spin-independent case, i.e. for σ SI p , as a function of the DM mass. The value of the coupling λ is now set to 1, and for the renormalization scale µ R we choose the central value of 2 GeV. The mediator mass M med is fixed to 2 TeV. In order to illustrate the    to one at small m χ and large M med to large values above two for large m χ and small M med , where, considering the different ranges that we explored, we approach the condition m χ = M med . This means that the ratio between the masses is close to one, and the growth of the K factor in this scenario is in agreement with what we observed in Fig. 4.
As expected from our discussion of Fig. 6, we observe a similar behavior for the SI and the SD case.

Comparison of limits from the LHC and direct detection experiments
We now move on to more phenomenological aspects and discuss the effects of the NLO Ref. [47,48]. In particular, for collider production of DM, the DM relic abundance is of no relevance. Therefore, in the M med -m χ plane showing exclusion limits for a specific In Fig. 8 we present exclusion limits for the DM-nucleon scattering cross sections σ SI p and σ SD n at LO and NLO accuracy as derived from ATLAS [33], together with exclusion limits provided by the CRESST III [30] and XENON1T [31,32]

Phenomenology of a combined s + t-channel model
where we have explicitly included the squared one-loop correction which is of order O(α 2 s ) and thus formally only contributes beyond NLO. Note that we did not include this term in any of the results discussed above. When c s,tree + c t,tree is small, the interference term between the tree-level and the one-loop contributions as well as the squared oneloop term can become the dominant contribution to the cross section. If the relative sign between c s,tree and c t,tree changes, the hierarchy of terms is reversed, and the one-loop correction ceases to be the dominant contribution to the cross section.
As no exclusion limits from LHC exist for the s + t-channel model, for illustrative purposes we assume fictional limits which exclude mediator masses up to M med = 500 GeV for all m χ < M med . Figure 9 illustrates the impact of interference effects for the example of the SI DM-proton cross section, σ SI p . For a positive sign of the g V q coupling, no cancellation effects occur at LO, and the NLO contributions constitute only a small correction to the LO results. However, when the sign of the g V q coupling is taken negative, the interference pattern of Eq. (18) has a strong impact on the fictional LHC exclusion limits. Because the LO contribution is artificially small, the interference of the tree-level contributions with the one-loop corrections as well as the squared one-loop term are the dominant contributions at NLO, and much larger than the pure LO result.
Consequently, the NLO corrections strongly modify the LO results.
Analogous effects are found for the SD DM-neutron cross section σ SD n , as shown in Fig. 10. We note, however, that due to an opposite sign of the tree-level Wilson coefficient
In our calculations, we are using conventional dimensional regularization with d = 4 − 2ε dimensions and anticommuting γ 5 . We use a hybrid on-shell/MS renormalization scheme: While we apply the on-shell scheme for the quark field-strength renormalization, the MS scheme is used for the renormalization of couplings and mediator masses.
We are working in the non-relativistic limit (q 2 → 0), where the center-of-mass energy squared of the DM-quark system is given asŝ = (m χ + m q ) 2 and the square of the momentum transfer between incoming DM particle and outgoing quark reduces tô In the following, we list the expressions for DM-quark scattering and keep the explicitû dependence. Corresponding expressions for anti-DM-quark scattering can easily be obtained thereof by applying the crossing relationû ↔ŝ.

A.1. Wilson coefficients for the s-channel model
In the s-channel model, the only tensor structures appearing in the tree-level matrix elements that contribute to DM-quark scattering are the ones corresponding to the vector and axial-vector currents. Imposing the matching condition of Eq. (14) leads to the following Wilson coefficients: with ∆ = 1 ε − γ E + ln(4π), where γ E is the Euler-Mascheroni constant γ E ≈ 0.57721. Note that here and in the following, the 1/ε term generally receives contributions both from UV and IR divergences. It turns out that the counterterm contributions are the same for both the vector and axial-vector currents, and furthermore they are exactly opposite to the corrections of the vector current: for all tensor structures j = S, V, A, T and with the color factor C F = 4/3. We note that, as mentioned before, while the sum of the vertex correction and counterterm contribution to the vector current cancels independently of the matching, the corresponding contribution to the axial-vector current is non-zero and vanishes only after including the operator correction.

A.2. Wilson coefficients for the t-channel model
Also in the t-channel model, only the vector and axial-vector operators contribute at tree level. We have used Fierz identities (see, e.g., Ref. [12]) to arrange the spinor fields in the same order as in the operators of Eq. (6). From the tree-level matching condition, we then obtain the following Wilson coefficients: Here we have introduced a common tree-level factor, c t,tree , which will also appear in the one-loop expressions for each tensor structure. It is quoted for our choice of model parameters, i.e. for t-channel mediators which only couple to left-handed quarks. A transition to the general case is straightforward, i.e. by replacing λ 2 Q L by a sum of the allowed couplings.
At O(α s ), we now have several contributions coming from propagator, vertex, and box corrections. We write these contributions in terms of scalar and tensor integrals in the notation of Ref. [55] and use the Mathematica extension Package-X to obtain analytical expressions for the loop integrals. All one-loop integrals emerging in our calculation can be expressed in terms of a small set of master integrals that below we abbreviate as For the three-point integrals, we encounter C i with i = 0, 1, 2, while in the case of four-point integrals, those of type D i and D ij with i, j = 0, 1, 2, 3 emerge.
The contributions from the propagator correction can then be written as: The contributions from the vertex corrections are: The contributions from the box correction are: We also list the counterterm contributions which we have calculated from the renormalization constants of our model: for each tensor structure j = S, V, A, T .

A.3. Wilson coefficients for the s + t-channel model
For the s + t-channel model, corresponding to the sum of the Lagrangians of Eq. (1) and Eq. (2), the tree-level and one-loop Wilson coefficients can be readily obtained by summing the contributions of the s-and t-channel models quoted above. In particular, this means that at the one-loop level, we must add the vertex and counterterm contribu- which corresponds to Eq. (28) without the operator correction.

A.4. One-loop operator corrections
As shown in the Feynman diagram of Fig. 2 (b), there is only one vertex correction plus a corresponding counterterm contributing to each one-loop operator correction in the EFT.
In the non-relativistic limit, these one-loop terms are proportional to the corresponding tree-level operators, for each tensor structure j = S, V, A, T , parameterizing the full one-loop correction and the color factor C F explicitly factored out.
In the following, we do not consider the one-loop corrections to the scalar and tensor and for the counterterms:  as defined in Eqs. (20) and (21). Summing the contributions to obtain the full one-loop operator correction as in Eq. (32), we see that the corrections vanish for the vector operator, and lead to a finite contribution for the axial-vector operator: