Simplified dark matter models with loop effects in direct detection and the constraints from indirect detection and collider search

Although the existence of Dark Matter (DM) has been established by substantial cosmological and astronomical observations, the microscopic nature of DM particles is still unknown. An appealing candidate for DM is the Weakly Interacting Massive Particle (WIMP) arising from various extensions of the Standard Model (SM). The experimental searches for WIMP consist of four main categories, i.e. the direct detection (DD) of possible scattering between DM and SM target materials, the indirect detection (ID) looking for signals of DM annihilation/decay products from the sky, the collider searches for signals from DM production at high energy accelerators, and the gravitational and/or cosmological effects originating from the DM in the early and/or the current Universe.

Among the aforementioned four categories, the DD experiments have achieved significantly improved sensitivity in the past two decades, but yielded up to now null results and very stringent bounds on the WIMP-nucleon scattering cross section. A natural explanation of the absence of a confirmed DM signal is that the scattering rate is highly suppressed by the typically small value of transfer momentum in the process, and/or the relative velocity between DM and the nucleon. A simple but compelling scenario resulting in the suppressed rate at tree-level is that fermionic DM particles $ \chi $ scatter off the target nucleon N through a pseudo-scalar mediator in the t-channel scattering process [1, 2]. The corresponding tree-level DM-nucleon contact interaction can reduce to a non-relativistic contact operator

in the non-relativistic limit. Here, $ {{s}}_\chi $ ($ {{s}}_N $) is the DM (target nucleon) spin, and the scattering exchange momentum $ {{q}} $ is only of the order of 10 MeV. As a result, this scenario leads to a momentum suppressed spin-dependent (SD) scattering cross section and thus an undetectable signal rate.

The suppression of tree-level scattering rate makes it appealing to further scrutinize the high-order effects from one-loop induced processes (see early discussions in e.g. [3, 4]), and the possibly detectable signals at the upgraded DD experiments. Integrating out the one-loop diagrams can induce distinct scalar-type operators giving non-momentum suppressed spin-independent (SI) scattering cross section. The enhancement of loop-level SI cross section by the squared total nucleon number in a nucleus competes with the loop suppression and may dominate the WIMP-nucleus cross section over the suppressed tree-level scattering. This one-loop effect in direct DM detection has been investigated in both simplified frameworks and UV complete models [2, 5-13]. Recent progress in the pseudo-scalar mediator scenario goes beyond the one-loop processes for scalar-type DM-quark operator $ m_q\bar{\chi}\chi \bar{q}q $ , and includes the dedicated contributions from two-loop scattering diagrams for scalar-type DM-gluon operator $ {\alpha_s\over \pi}\bar{\chi}\chi GG $ , after integrating out both the heavy quarks and the mediator. It has been shown in Ref. [14], in the framework of the two Higgs doublet model, that the full two-loop calculations deviate considerably from the result obtained by the conventional relation between the scalar-type current for heavy quarks and the gluon. This discrepancy is caused by the failure of the quark momentum expansion for heavy quarks, and by the fact that the two-loop diagram for gluon emission is ignored when one uses the relation for DM-gluon scattering.

In this work we revisit the loop effect in DD in simplified DM models, including either a pseudo-scalar mediator or maximal CP violation. We consider the hypotheses with only one spin-0 mediator coupled to the SM quarks and fermionic DM particles [15, 16], which leads to momentum-suppressed WIMP-nucleus scattering cross section at tree-level. The latest approach for dedicated loop calculations is used for high-order contributions to the cross section, together with the estimate of the running effects of renormalization group equations (RGE). We find that the differential event rates with respect to recoil energy can provide particular spectral shapes at different energies induced by loop corrections. Besides the loop effects on direct DM detection in pseudo-scalar mediated models discussed in the recent works [14, 17], we take into account the constraints from the other DM detection categories, e.g. the DM relic abundance, the ID constraint in terms of gamma-ray emission, as well as the current status of collider search. These synergistic studies are regarded as completion and improvements of Refs. [14, 17].

The paper is organized as follows. In Sec. 2 we describe the simplified dark matter models. We then give the effective DM-nucleon interactions at tree-level and the corresponding DM-nucleus scattering cross sections. In Sec. 3, we present the effective operators and the Wilson coefficients at loop-level for the DM-nucleon cross section. The numerical results are given in Sec. 5. Our conclusions are drawn in Sec. 6. Some technical details of loop calculations are collected in the Appendix.

In this work, we focus on the simplified DM framework which consists of Majorana fermion DM $ \chi $ and a spin-0 mediator a coupled to $ \chi $ , and the SM quarks with strength $ g_\chi $ and $ g_q $ , respectively. We consider each of the following three scenarios at a time

Here, the $ a\bar{q}q $ coupling is also scaled by the SM-like Yukawa coupling, where $ v_0 = 246 $ GeV is the SM Higgs vacuum expectation value. The model D4 with a a pure pseudo-scalar, is designated as the pseudo-scalar mediator DM model in most literature. Models D2 and D3 are induced by more specific UV complete models with CP violation [17-22], and correspond to the case with maximal CP violation.

Based on the DM interactions with quarks and gluons at tree-level, the DM-nucleon contact interactions are described by the effective Lagrangians as follows

where the tree-level coefficients are defined as

The coefficient $ C_N^{\rm tree}({\rm D_3}) $ is equal to $ C_N^{\rm tree}({\rm D_4}) $ , as models D3 and D4 share the same quark bilinear form $ \bar{q}\gamma_5 q $. Here, $ \Delta_q^{N} $, $ f_q^{N} $ and $ f_G^{N} $ are quark/gluon-nucleon form factors, as numerically used in micrOMEGAs [23]. The quark level constant is defined as $ C_q = g_\chi g_q{m_q\over v_0} $. Consequently, the differential DM-nucleus scattering cross sections read as

where $ m_T $ is the nucleus mass, v is the DM speed in the Earth frame, $ E_R $ is the nuclear recoil energy and $ F_{\Sigma''}^{(N,N')}(q^2), F_{M}^{(N,N')}(q^2) $ are the form factors defined in Ref. [24]. The tree-level WIMP-nucleus scattering cross sections in the above simplified models are all dependent on the transfer momentum $ q = \sqrt{2E_R m_T} $. As seen above, they are suppressed by $ m_N^2q^2/ m_a^4 $, $ m_\chi^2 q^2/ m_a^4 $ and $ q^4/m_a^4 $ for models D2, D3 and D4, respectively.

In this section, we derive the loop-level effect in direct DM detection in the above simplified models, followed by the estimate of the scale effects in terms of renormalization group evolution.

The general Lagrangian for the non-momentum suppressed DM-nucleon SI cross section is given by

where $ {\cal{O}}^q_{\mu\nu} = {i\over 2}\bar{q}\left(\partial_\mu\gamma_\nu+\partial_\nu\gamma_\mu-{1\over 2}g_{\mu\nu}{\not\!\! \partial }\right)q $ is the twist-2 operator. For the models we consider, the Wilson coefficients in Eq. (11) are all zero at tree-level but can be generated at loop-level, denoted by $ C'_q = C_q^{\rm box} $, $ C_G = C_G^{\rm box} $, $ C_q^{(1)} = C_q^{\rm (1)box} $, $ C_q^{(2)} = C_q^{\rm (2)box} $. The coefficients for scalar-type DM-quark operator and the twist-2 operator, i.e. $ C_q^{\rm box} $, $ C_q^{\rm (1)box} $ and $ C_q^{\rm (2)box} $, are generated by the box diagrams in the top panels of Fig. 1. The two-loop diagrams in Fig. 1 , with only the heavy quark Q in the quark loop, contribute to the scalar-type DM-gluon operator and the coefficient $ C_G^{\rm box} $.

Figure 1.  Loop diagrams for the DM-quark currents (top) and DM-gluon currents (bottom).

Following the non-relativistic limit used in Ref. [14], we expand the small momentum of valence quarks in the amplitude of the DM-quark scattering box diagrams. The coefficients $ C_q^{\rm box} $, $ C_q^{\rm (1)box} $ and $ C_q^{\rm (2)box} $ are then obtained by reading the DM-quark effective operators. For the DM-gluon coefficient $ C_G^{\rm box} $, one needs to calculate the amplitude of two-loop diagrams and find the effective operator $ \bar{\chi}\chi G^a_{\mu\nu}G^{a\mu\nu} $. The complete two-loop calculations ensure the validity of the obtained $ C_G^{\rm box} $ for any value of mediator mass $ m_a $. For model D4, the above Wilson coefficients are equivalent to those in the pseudo-scalar mediator model as derived in Ref. [14]

where the loop functions $ X_{001}, X_{111}, F $ are given in Ref. [14] and the references therein. Following the same procedure, we obtain the corresponding Wilson coefficients for models D2 and D3, which are related to those in model D4

The new loop functions are collected in Appendix. In our numerical calculations, we use Package-X [25] to compute the above loop functions.

Based on the above effective operators for SI DM-nucleon scattering and the corresponding Wilson coefficients, we define the DM-nucleon constant at loop-level

where the second moments of the parton distribution functions for quarks $ q^N(2) $ and anti-quarks $ \bar{q}^N(2) $ are taken from CTEQ PDFs [26]. The SI cross section for the DM interaction with a nucleon is thus given by

In terms of the form factor function $ F_{M}^{(N,N')}(q^2) $, the differential SI cross section of the DM interaction with a nucleus with mass $ m_T $ is

Another manifestation of the loop effect is the mixing of operators according to RGE. The RGE effect can be important if one considers DM phenomenology at vastly different energy scales. For instance, the DM annihilation typically occurs at the electroweak scale, while the DM particles are possibly produced near the TeV scale at colliders. The energy scale for the DM-nucleon scattering in DD experiments is of the order of the hadron scale $ \mu_{\rm had} $. One usually starts with a gauge-invariant renormalizable DM model defined near or above the electroweak (EW) scale $ \mu_{\rm EW}\simeq m_Z $, but studies the non-relativistic DM-nucleon scattering rate at a characteristic scale $ \mu_{\rm had}\simeq 1 $ GeV. A series of effective field theories (EFT) should be properly constructed by integrating out particles heavier than the current EFT scale $ \mu_{\rm EFT} $ and reasonably matched when passing the thresholds of particles lighter than $ \mu_{\rm EFT} $ , where they are integrated out in a similar way. Between the thresholds, the evolution and mixings of EFT operators should be performed according to RGE. The scale of the first EFT constructed in the whole analysis determines the procedure of RGE and threshold matching. The above procedures have been well elaborated in e.g. [27-31], and implemented in packages such as DirectDM [32], runDM [33], Wilson [34] for specific or generic models.

An important difference should be emphasized between the scalar-type and twist-2 operators in Eq. (11). The scalar-type form factors $ f_q^{N} $ for the light quarks $ q = u,d,s $ are attributed to the non-perturbative QCD effects with energy scale around or below 1 GeV, and are obtained from the lattice QCD simulations. Thus, the scalar-type Wilson coefficients for the light quarks $ u,d,s $ and gluon must take values around 1 GeV. Depending on the scale of the first EFT constructed in the whole analysis, e.g. at $ \mu = m_Z $, this implies the procedures of RGE and threshold matching when calculating the scalar-type operator contributions. At the scale of about 1 GeV, the heavy quarks $ Q = c,b,t $ have been integrated out into the scalar-type gluon operator using the full two-loop calculations, as emphasized in Refs. [14, 17]. In turn, the twist-2 form factors $ q^N(2), \bar{q}^N(2) $ in Eq. (11) can be calculated perturbatively using parton PDFs at various scales [3, 27-31], e.g. 1 GeV or $ m_Z $. One can choose a convenient scale to calculate the twist-2 contributions, with the proper active field content (e.g. 5 flavor quarks $ u,d,s,c,b $ at $ \mu = m_Z $), and the Wilson coefficients and form factors evaluated at that scale. Note that we ignored the negligible contributions from the twist-2 gluon operator, since its Wilson coefficient is suppressed by an additional $ \alpha_s/\pi $ due to the operator definition [28, 29].

We note that in Refs. [14, 17], the values of $ q^N(2), \bar{q}^N(2) $ are evaluated at $ \mu = m_Z $ in the calculations. Since the same (similar) box diagrams and model parameters are used to obtain the scalar-type and twist-2 Wilson coefficients for quarks (gluons), a more consistent implementation should involve the RGE running effects for the scale-type operators from $ m_Z $ to 1 GeV as discussed above. Since the coupling between the mediator and the SM quarks is chosen to mimic the SM Yukawa structure, the DM-nucleon constant from scalar-type DM-gluon interaction $ C_G f_{G}^N $ in Eq. (22) is dominated by the top quark loop. The constant $ C_G f_{G}^N $ also dominates over the scalar-type and twist-2 DM-quark interactions. To have a conservative estimate of the scale effects on the scalar type DM-gluon operator $ -{9\alpha_s\over 8\pi}\bar{\chi}\chi G^a_{\mu\nu}G^{a\mu\nu} $, we used the package DirectDM to perform its RGE running from $ m_Z $ to 1 GeV. We found that the scale effects give a negative correction of 1%~2%, and thus do not affect our main conclusions in this work.

In this section, we consider the DM constraints in the above simplified models from the other categories mentioned in the Introduction, including the relic abundance, indirect detection and collider search.

Assuming that the DM particles have frozen out in the early Universe as standard thermal relics, they acquire their present abundance through annihilation processes. The pair of DM particles $ \chi $ in the simplified models can either annihilate into SM quark or gluon pairs via the s-channel processes $ \chi \chi\to a\to q\bar{q}, gg $ , or annihilate into two mediators $ \chi \chi\to a a $ when kinematically allowed [15, 16, 35-37]. The amplitudes of the two annihilation channels are governed by $ g_\chi g_q $ and $ g_\chi^2 $, respectively. We assume all kinematically accessible final states of DM annihilation and use micrOMEGAs 5.0 [38] to calculate the relic abundance. Note that the WIMP candidate may account for only a fraction of the total DM of the Universe, referred to as multi-component DM scenario [39, 40]. In this scenario, the DM energy density measured by PLANCK [41] is imposed as an upper limit on the WIMP relic abundance.

Dwarf galaxies are the search targets for DM annihilation into gamma rays. The Fermi Large Area Telescope (LAT) has not detected any excess of gamma ray emission from the dwarf spheroidal satellite galaxies (dSphs) of the Milky Way. Thus, an upper limit on the DM annihilation cross section can be placed from a combined analysis of multiple Milky Way dSphs [42, 43]. For individual dwarf galaxy target, Fermi-LAT provided tabulated values of delta-log likelihood as a function of the energy flux bin-by-bin. The gamma ray energy flux from DM annihilation for the jth energy bin and the kth dwarf is given by

where $ J_k $ is the J factor for the kth dwarf. The energy flux only depends on $ m_{\chi} $, $ \langle \sigma v\rangle $ and $ J_k $, and is thus calculable for the DM annihilation process from the above simplified models. We use the PPPC4DMID package [44] to obtain the spectrum of photons $ {\rm d}N_\gamma/{\rm d}E $. The likelihood for kth dwarf is

where $ {\cal{L}}_{j,k} $ is the tabulated likelihood provided by Fermi-LAT for each dwarf and energy flux. The uncertainty of the J factors is taken into account by profiling over $ J_k $ in the likelihood below [42]

with the measured J factor $ \bar{J}_k $ and error $ \sigma_k $. A joint likelihood for all dwarfs can then be performed as

where $ \mathbb{J} $ is the set of J factors $ J_k $. In our numerical implementation, we adopt the corresponding values of $ {\cal{L}}_{j,k} $ and $ \bar{J}_k, \sigma_k $ for 19 dwarf galaxies considered in Ref. [43].

According to the maximum likelihood analysis adopted by Fermi-LAT, the delta-log likelihood is given by

where $ \widehat{\langle \sigma v\rangle} $ and $ \widehat{\mathbb{J}} $ maximize the likelihood at any given $ m_{\chi} $, and $ \widehat{\widehat{\mathbb{J}}} $ maximizes the likelihood for given $ m_{\chi} $ and $ \langle \sigma v\rangle $. The 95% C.L. upper limit of the annihilation cross section for a given $ m_{\chi} $ is determined by demanding $ -2\Delta \ln{\cal{L}}(m_{\chi},\langle \sigma v\rangle)\leqslant 2.71 $. We perform the likelihood analysis and obtain the upper limit using Minuit [45]. Once the annihilation cross section calculated by a certain set of model parameters is larger than the limit, we claim that the corresponding parameter values are excluded by Fermi-LAT dSphs.

The Large Hadron Collider (LHC) performed the search for DM in association with energetic jets [46, 47], or the third generation quarks [48-50] for simplified DM models with spin-0 mediator at $ \sqrt{s} = 13 $ TeV collisions. The most severe limits are from the final states with $ t\bar{t} $ and missing transverse momentum [51]. For model D4 with pseudo-scalar mediator, assuming unitary couplings $ g_\chi = g_q = 1 $, the range of mediator mass between 15 and 25 GeV is excluded [49]. This limit is valid for all DM masses as long as $ m_a>2m_\chi $ . Closely related models D2 and D3 should have a very similar collider constraint.

By combining the theoretical calculations of DM-nucleus scattering cross sections with a certain velocity distribution for DM particles, we can calculate WIMP signal rates for DD experiments. The differential event rate with respect to the recoil energy is given by

where $ \rho_\chi $ is the local DM density, which is fixed to $ 0.3 \ {\rm GeV}/{\rm cm}^3 $. We take the cut-off Maxwell-Boltzmann velocity distribution for $ f(\vec{v}) $ with the escape velocity $ v_{\rm esc} = 544 \ {\rm km}/{\rm s} $. Together with the differential scattering cross sections obtained above, the differential event rate can be evaluated. We employ DMFormFactor [24, 52] for the numerical calculations on the Xenon nucleus $ ^{129}{\rm Xe} $. The recoil energy spectra of models D2, D3 and D4 are shown in Figs. 2, 3 and 4, respectively. Model D2 leads to SI scattering cross section with a strong enhancement for large nuclei. Thus, although there is momentum suppression at tree-level, the spectrum at tree-level for model D2 still dominates over the loop-level contribution. As it is suppressed by $ q^4 $ at tree-level, the loop-level spectrum is much greater than at tree-level for model D4. Although models D3 and D4 both lead to SD cross section at tree-level, this discrepancy is smaller in model D3 as its tree-level scattering cross section is suppressed by $ q^2 $ only. The loop-level spectrum of model D3 dominates over the tree-level only in the range of low recoil energy.

Figure 2.  (color online) Recoil energy spectra for model D2 with $ m_a = 10 $ GeV, $ g_\chi g_q = 1 $ and $ m_\chi = 40 $ GeV (left) or $ m_\chi = 400 $ GeV (right). The tree-level and loop-level spectra are denoted by red and black curves, respectively.

Figure 3.  (color online) Recoil energy spectra for model D3, as labeled in Fig. 2.

Figure 4.  (color online) Recoil energy spectra for model D4, as labeled in Fig. 2.

As the recoil spectrum of the SI scattering induced by loop diagrams is dominant in model D4, the prediction of SI DM-nucleon cross section in Eq. (23) can be compared directly to the limits set by DD experiments to yield a bound on $ g_\chi g_q $. As shown in Fig. 5, for $ g_\chi g_q \leqslant 1 $ in model D4, the SI scattering cross sections in the red region are below the XENON1T exclusion limit [53, 54] but above the neutrino floor. The green region gives cross sections below the neutrino floor. Future DD experiments could thus be more sensitive to the region of $ g_\chi g_q > 0.4 $ and $ m_\chi < 200 $ GeV in the case of $ m_a = 10 $ GeV.

Figure 5.  (color online) The region $ g_\chi g_q $ vs. $ m_\chi $ below the Xenon1T exclusion limit and above the neutrino floor (red), and the region below the neutrino floor (green) for model D4, with $ m_a = 10 $ GeV (left) and $ m_a = 30 $ GeV (right). Assuming $ g_\chi = 1 $, the light blue values are excluded by Fermi-LAT dSphs, and the DM relic abundance satisfies $ \Omega h^2\leqslant 0.12 $ above the black curves.

The cross section for DM annihilation into SM quarks in model D4 is proportional to $ m_q^2/v_0^2 $, and thus the $ t\bar{t} $ channel dominates if kinematically allowed. The annihilation into mediator pairs is governed by $ m_a/m_\chi $ and plays a crucial role in the small $ m_\chi $ range. Fermi-LAT dSphs excludes a major part of the parameter region $ m_{\chi}\lesssim 100 $ GeV and $ m_{\chi}\gtrsim m_t $, as shown by the light blue region in Fig. 5. The region around $ m_{\chi}\simeq 100 $ GeV is not subject to the ID constraint due to the fact that the channel $ \chi\chi\to a\to t\bar{t} $ is not kinematically allowed and the annihilation into mediator pairs is suppressed [36, 37]. Above the black contours in Fig. 5, the DM relic abundance satisfies $ \Omega h^2\leqslant 0.12 $. In the multi-component DM scenario, other DM particles are required to explain the measured value of the thermal DM energy density in the region above the black curves. In that case the dwarf spheroidal galaxies contain not only $ \chi $ but also other particles. The number density of $ \chi $ is thus smaller than the value estimated from the J-factors. We expect the ID constraint on $ \chi $ in the multi-component scenario to be less stringent than the exclusion shown in Fig. 5. The detailed study of the multi-component DM scenario is beyond the scope of this work.

We reexamined the loop-level correction of the WIMP-nucleon scattering cross section in the framework of simplified DM models. A spin-0 mediator is assumed to couple with fermionic DM particle and the SM quarks in each model. The simplified models include either a pseudo-scalar mediator or maximal CP violation. The cross sections in these models are low-momentum suppressed at tree-level, but receive sizable loop-induced SI contribution. Following the recent progress in loop-level correction, we perform one-loop calculations for scalar-type and twist-2 DM-quark operators, and complete two-loop calculations for scalar-type DM-gluon operator. The differential event rates with respect to recoil energy could provide particular spectral shapes induced by loop correction in models with pseudo-scalar mediator or CP violation. By including the loop-level SI cross section, we found that future DD experiments could be sensitive to a fraction of parameter space with $ g_\chi g_q > 0.4 \ (0.8) $ and $ m_\chi < 200 \ (80) $ GeV for the pseudo-scalar mass of 10 (30) GeV, which gives no detectable signal with only tree-level contribution. The sensitivity of DD experiments to these models is also complementary to the constraints from ID and collider search.

T. L. would like to thank Tomohiro Abe, Motoko Fujiwara and Junji Hisano for helpful discussion. P. W. would like to thank Chengcheng Han for the beneficial argument.

The loop functions that are used for models D2 and ${\rm{D3}} $ are

and

with

The $ Y_i (i = 1,2,3) $ functions are defined in Ref. [14] and the references therein.