Recasting direct detection limits within micrOMEGAs and implication for non-standard dark matter scenarios

Abstract

Direct detection experiments obtain 90% upper limits on the elastic scattering cross sections of dark matter with nucleons assuming point-like interactions and standard astrophysical and cosmological parameters. In this paper we provide a recasting of the limits from XENON1T, PICO-60, CRESST-III and DarkSide-50 and include them in micrOMEGAs. The code can then be used to directly impose constraints from these experiments on generic dark matter models under different assumptions about the DM velocity distribution or on the nucleus form factors. Moreover, new limits on the elastic scattering cross sections can be obtained in the presence of a light t-channel mediator or of millicharged particles.

1 Introduction

Searches for dark matter (DM) through direct detection (DD) experiments have been pursued actively for decades [1,2,3,4,5,6,7,8]. None of the experiments with a good signal/background discrimination have found evidence for DM, thus could only set upper limits on the DM elastic scattering cross section on nucleons. For DM masses above roughly 3 GeV, the best limits for spin-independent (SI) interactions are currently obtained by XENON1T [1, 9]. For lower masses, searches are more challenging and require a very low threshold for nuclear recoil energy, thus the limits are typically much weaker. Currently the best limits are obtained from DarkSide [3], and CRESST [10] and a series of projects are concentrating their efforts in improving the reach at or even below the GeV [11, 12] in particular by using DM scattering on electrons [13,14,15,16,17,18]. For spin-dependent interactions on neutrons and protons, currently the best limits are obtained by XENON1T [2] and PICO-60 [4, 5] respectively. Currently, limits are generally interpreted in terms of DM elastic scattering on nucleons through a mediator with a mass much larger than the typical momentum exchange. Moreover they are obtained assuming equal proton and neutron spin-independent cross sections and for a specific choice of astrophysical parameters, notably that the DM velocity distribution is Maxwellian.

Although traditional WIMP models feature mediators at or above the electroweak scale (e.g., a Higgs, Z, a new boson or a new coloured particle), new classes of DM models have relinquished the link with the electroweak scale thus considerably extending the range of masses for both DM and mediators. In particular models with a very light mediator have been considered [19,20,21]. The motivation for a light mediator include the possibility to provide strong dark matter self-interactions and explain anomalies in galaxy clusters [22,23,24] as well as the possibility to enhance the direct detection signal in models with feebly coupled particles [25].

While it is straightforward for the experimental collaborations to obtain limits within a framework different than the default one chosen, the corresponding code is not publicly available. For example only PandaX [7, 26] and more recently XENON1T [9] have published limits obtained for both heavy and light mediators. Our goal is precisely to provide a tool that allows to reinterpret the 90% limits obtained by the experimental collaborations within their specific framework and apply them to a wider set of DM models and DM velocity distributions. The code is developed as a module of micrOMEGAs  [27, 28]. In this first version, recasts of the limits from XENON1T [1], DarkSide-50 [3], PICO-60 [5] and CRESST-III [10] are provided. These thus provide the best limits for the cases of spin independent and spin dependent interactions in neutrons and protons for DM masses above 1 GeV. Based on this recast, we give typical examples on how the code can be used to set limits on new models. The models considered include the case of a light mediator, in particular a $Z'$, as well as millicharged particles. Moreover the impact of alternate velocity distributions is analysed. Recasting of these limits as well as other recent direct detection experiments are also included in DDCalc [29,30,31] and in SuperIso [32]. Note that both these recasting reproduce well the XENON1T exclusions for DM masses at the weak scale or above, however they feature significant differences for masses near the sensitivity threshold when events are expected in the region at low nuclear recoil energy which is particularly challenging for experiments. Our implementation provides a better match to XENON1T in the case of light dark matter as will be described in the next section. Moreover since a heavy DM with a light mediator features a recoil energy distribution that resembles that of a light dark matter, in the sense that it peaks at smaller energies than the corresponding one for a heavy mediator, we expect a more reliable recast for the light mediator case. Considering the lack of complete information on the experimental data, for XENON1T we adopt a strategy which consists in tuning the efficiency for nuclear recoils in order to reproduce the SI experimental limit for DM interactions at all masses. We refer to this approach as ’inverse recasting’. Note that our approach can only be applied to the case where the DM signal is dominant at low recoil energy as will be discussed in Sect. 3.1. For exotic signals with interactions at large recoil energy, for example the ones studied in [33,34,35], our approach cannot be applied as it would lead to limits on the exclusion cross- section that are not severe enough. For other experiments we simply use the information provided in the publications to describe the detector efficiency and the background to reproduce the experimental limit.

The paper is organised as follows. After describing the formalism for the event rates in direct detection in Sect. 2, we describe our reconstruction of the XENON1T, DarkSide-50, PICO-60 and CRESST-III experimental limits on SI interactions in Sect. 3 and SD ones in Sect. 4. In Sect. 5 we show how these recasts allow to obtain limits in specific models involving a light mediator, a millicharge DM as well as generic DM velocity distribution. Section 6 contains our conclusions. All results obtained in our paper can be reproduced using the new micrOMEGAs functions described in the Appendix A.

Fig. 1

Left: predictions for the recoil energy distribution of Xenon nuclei for an exposure ${\mathcal {L}}=279\times 900$ kg$\cdot$days and $\sigma ^{90}=248.6,5.39,0.566,0.448,0.912,1.71\times 10^{-46}\,\mathrm{cm}^2$ for $M_\chi =6,10,20,40,100,200$ GeV respectively. Right: recoil energy distributions convoluted with the XENON1T acceptance $p_{Xe}(E)$ for $M_\chi =35(200)$ GeV and $\sigma ^{90}=4.71 (17.1)\times 10^{-47}\, \mathrm{cm}^2$

2 Dark matter scattering on nuclei

We first review the standard formalism for obtaining the nuclear recoil energy distribution for DM scattering on nuclei, relevant for direct detection experiments. Since the velocity of DM particles is about $v_0\approx 0.001c$, the maximum velocity of the nucleus that recoils cannot exceed 2$v_0$. Thus, the maximum transferred momentum in DM-nucleus collision is $q_{max}=2 v_0 M_A\approx 200$ MeV for a nucleus mass $M_A \approx 100 \mathrm{GeV}$. At such low momentum transfer, DM-nucleon interactions can be described by an effective Lagrangian leading to constant matrix elements. Moreover the amplitudes can be divided into spin-dependent (SD) and spin-independent (SI) interactions which do not interfere. The DM-nuclei interactions are simply related to the DM-nucleon interactions after introducing a nucleus form factor which depends on the momentum transfer $q=\sqrt{2 M_A E}$ where E is the nucleus recoil energy. The energy distribution of a recoil nuclei A produced by SI interaction with DM in a detector with total mass $M_{det}$ and exposure time T reads [27, 36]

$$\begin{aligned} \frac{dN^{SI}_A}{dE}= & {} \frac{2}{\pi }M_{det} T \frac{\rho _{\chi }}{M_\chi } I(E) (\lambda _p Z + \lambda _n(A-Z))^2 F_A^2(q)\nonumber \\ \end{aligned}$$

(1)

where Z and A are the atomic number and mass number of the detector material, $M_\chi$ is the DM mass, $\rho _{\chi }$ the DM local density, and $\lambda _N$ are DM-nucleon scattering amplitudes. SI interactions are typically generated from effective scalar or vector interactions of DM with nucleons. For example for an effective scalar interaction of Majorana fermions with nucleons N, ${\mathcal {L}}=\lambda _N \bar{\chi }\chi \bar{\psi }_N{\psi }_N$, the SI DM-nucleon cross section is given by

$$\begin{aligned} \sigma _{\chi N}^{SI} = \frac{4}{\pi } \mu _{\chi N}^2\lambda _N^2,\quad N=n,p \end{aligned}$$

(2)

where $\mu _{\chi N} = M_\chi M_N/(M_\chi +M_N)$ is the DM-nucleon reduced mass. The event rate also depends on the nucleus form factor, $F_A(q)$ and on the velocity distribution through,

$$\begin{aligned} I(E) = \int \limits _{\sqrt{E M_A/(2\mu _{\chi A}^2)} }^\infty \frac{f(v)}{v} dv, \end{aligned}$$

(3)

where f(v) is the DM velocity distribution in the detector rest frame normalized such that

$$\begin{aligned} \int \limits _0^\infty f(v)dv =1 \end{aligned}$$

(4)

The recoil energy distributions for various DM masses are displayed in Fig. 1-left.

In direct detection experiments after analysing the number of registered events and estimating the background, limits are set on $\sigma ^{SI}_{\chi p}$ assuming $\sigma ^{SI}_{\chi p}=\sigma ^{SI}_{\chi n}$. All experiments also assume a value for the DM local density near the Sun, $\rho _{\chi }=0.3~\mathrm{GeV}/\mathrm{cm}^3$, and a Maxwellian DM velocity distribution defined with the parameters

$$\begin{aligned} v_{Rot}=220~\mathrm{km/s},\quad v_{esc}=544~\mathrm{km/s},\quad v_{Earth}=232~\mathrm{km/s} \end{aligned}$$

(5)

where $v_{Rot}$, the rotation velocity of the Galaxy and $v_{esc}$, the escape velocity in the galaxy, characterize the DM velocity distribution in the Milky Way [36]. $v_{Earth}$ is the velocity of the Earth in the galactic frame.

The energy distribution of recoil events resulting from SD interactions of DM with nuclei in a detector with mass $M_{det}$ and exposure time T reads [27, 36,37,38]

$$\begin{aligned} \frac{dN^{SD}_A}{dE}= & {} {\mathcal M}_{det} T \frac{\rho _{\chi }}{M_\chi } I(E) \frac{8}{2 J_A+1}\Bigg (S_{00}(q)(\xi _p+\xi _n)^2 \nonumber \\&+S_{01}(q)(\xi _p^2-\xi _n^2)+ S_{11}(q)(\xi _p-\xi _n)^2\Bigg ) \end{aligned}$$

(6)

where $J_A$ represents the spin of the detector material, $\xi _{p,n}$ are the DM-nucleon amplitudes normalized such that

$$\begin{aligned} \sigma _{\chi {N}}^{SD} = \frac{12}{\pi }\mu _{\chi {N}}^2\xi _{N}^2. \end{aligned}$$

(7)

For example, an effective axial-vector interaction of Majorana fermions with nucleons, ${\mathcal {L}}=\xi _N \bar{\chi }\gamma _\mu \gamma _ 5\chi \bar{\psi }_N \gamma ^\mu \gamma _ 5{\psi }_N$ will lead to the above cross section while for Dirac fermions, the same cross section is obtained for a Lagrangian defined with $\xi _N \rightarrow 2 \xi _N$. $S_{ij}(q)$ are the nucleus SD form factors. Calculations or these form factors within nuclear models are reviewed in  [37] and more recent calculations are available in [38]. Another set of form factors is currently used by experimental collaborations, these form factors, $F_{44}^{ab}$, are defined in the effective field theory approach in Ref. [39], they are expressed as

$$\begin{aligned} F_{44}^{ab}= \frac{J_A(J_A+1)}{12}\left( F_{\Sigma ^{'}}^{ab}+F_{\Sigma ^{''}}^{ab}\right) \end{aligned}$$

(8)

and analytical expressions for $F_{\Sigma ^{'}}^{ab}, F_{\Sigma ^{''}}^{ab}$ can be found in the Appendix of Ref. [39]. Simple expressions allow to relate these form factors with those in Eq. 6.

$$\begin{aligned} F_{44}^{pp}(q)= & {} \frac{\pi }{4(2J_A+1)} \left( S_{00}(q)+S_{11}(q)+ S_{01}(q)\right) , \nonumber \\ F_{44}^{nn}(q)= & {} \frac{\pi }{4(2J_A+1)} \left( S_{00}(q)+S_{11}(q) - S_{01}(q)\right) \nonumber \\ F_{44}^{pn}(q)= & {} F_{44}^{np}(q)= \frac{\pi }{4(2J_A+1)} \left( S_{00}(q) - S_{11}(q)\right) . \end{aligned}$$

(9)

3 Spin-independent interactions: recasting experimental exclusions

3.1 XENON1T

To repeat exactly the XENON1T analysis would require detailed information on events distribution, background estimation, and the use of nuisance parameters for all points of event space characterized by scintillation signals cS1, cS2b and interaction positions Z and R [1]. Lacking this detailed information we propose instead to reconstruct an effective efficiency by using the 90% exclusion cross section obtained by XENON1T from their complete analysis, this will then be validated by comparing with the XENON1T upper limit as will be explained below. This approach can be considered as a simplified version of the XENON1T analysis where some cuts in cS1, cS2b space are applied to increase the signal/background ratio. Our simplified approach relies on the observation that in some of the subspaces where XENON1T reported signal and background best-fit values, XENON1T can be considered as a low background experiment. Specifically we will use the reference detector mass of 0.9 t. In this subspace illustrated in Figs. 2 and 3 of Ref. [1] both the electromagnetic and neutron background are suppressed. In this region, XENON1T reports two detected events and an estimated background, $n_B=1.62$ events.

We first have to check that relying only on partial data is a reasonable assumption to approximately reproduce the upper limit obtained in the full analysis. For this purpose we use the data for the best-fit point presented in Table 1 of [1]. For this reference point the DM mass, $M_R$, the cross section, $\sigma _R$, and the expected number of signal events after applying cuts, $n_R$, are given by

$$\begin{aligned} M_R=200\,\mathrm{GeV},\quad \sigma _R=4.7\times 10^{-47} \,\mathrm{cm}^2,\quad n_R=1.16. \end{aligned}$$

(10)

Using the Feldman–Cousins formula we can easily estimate the cross section required for a 90% exclusion, we find $\sigma = 1.65\times 10^{-46}\,\mathrm{cm}^2$, a value close to the one obtained by XENON1T, $\sigma =1.73\times 10^{-46}\,\mathrm{cm}^2$. Thus we conclude that the XENON1T data obtained after imposing cuts is suitable for obtaining upper limits on the DM-nucleon cross section.

In general to recast the result of a DD experiment while lacking the full information on signal events, cuts, backgrounds and the associated uncertainties, one needs at least to know the detection efficiency p(E) and the background distribution as function of the nucleus recoil energy after cuts. The efficiency of XENON1T, which we denote $p_{Xe}$, is shown in Fig. 1 of Ref. [1]. We use the efficiency of the second science Run, SR1. However, this efficiency does not include the effect of cS1, cS2b cuts. Indeed the number of signal events for the best-fit point obtained with this efficiency and for the full detector mass 1.3 t, is $n=3.56$ which corresponds to the number of DM signals before cuts cited in Table I [1]. The same table shows that this number is reduced by a factor 1.7/3.56 after cuts.¹ Moreover we note that using the efficiency $p_{Xe}$ for the excluded signal for a DM of 6 GeV ($\sigma =2.8 \times 10^{-8}\mathrm{pb}$) we obtain only 1.3 events, a number insufficient for a 90% exclusion. Thus we choose not to use directly $p_{Xe}$ and instead propose to reconstruct an effective efficiency by using the 90% exclusion cross section obtained by XENON1T from their complete analysis, this will then be validated by comparing with the XENON1T upper limit.

For a low background experiment it is reasonable to use non-binned likelihood

$$\begin{aligned}&L(p,\frac{dN_{\chi }(M_{\chi },\sigma )}{dE}) \nonumber \\&\quad = e^{- {\mathcal {L}} \int \limits _0^{E_{max}} \left( p(E)\left( \frac{dN_{\chi }(M_{\chi },\sigma )}{dE} + b(E)\right) +b_\gamma (E) \right) dE} \prod \limits _{k\in events}\left( b_k +s_k \right) \nonumber \\&\quad =\left[ e^{- {\mathcal {L}} \int \limits _0^{E_{max}} (p(E) b(E) +b_\gamma (E))dE} \prod \limits _{k\in events} {b_k}\right] \nonumber \\&\qquad \times e^{- {\mathcal {L}} \int \limits _0^{E_{max}} p(E)\frac{dN_{\chi }(M_{\chi },\sigma )}{dE} dE} \prod \limits _{k\in events}\left( 1+ \frac{s_k}{b_k} \right) \end{aligned}$$

(11)

where $dN_{\chi }/dE$ is the nuclei recoil energy distribution corresponding to the scattering of a DM of mass $M_\chi$ with a cross section $\sigma$, p(E) is the efficiency for the detection of signal events, ${\mathcal {L}}$ is the exposure, b(E) and $b_\gamma (E)$ are the neutron and electromagnetic background distribution, and $s_k= \frac{dN_{\chi }(M_{\chi },\sigma )}{dE_k}$ and $b_k$ are the signal and background probability distribution function (p.d.f.) for each detected event. Using a Bayesian approach with flat priors, we determine the credible interval for the cross section [$0,\sigma _{ex}$] corresponding to a fraction $1- \alpha$ of the posterior probability where,

$$\begin{aligned} \alpha (\sigma _{ex})=\frac{ \int _{\sigma _{ex}}^\infty L(p, {dN(M_{\chi }, \sigma )}/{dE})d\sigma }{ \int _0^{\infty } L(p, { dN(M_{\chi }, \sigma )}/{dE} ) d\sigma } \end{aligned}$$

(12)

Note that when nuisance parameters are taken into account when estimating the background, one has to integrate both the numerator and denominator in Eq. 12 with some prior. In the following we ignore such nuisance parameters, thus the term enclosed in squared brackets in Eq. 11 cancels out. In this approximation, the background contributes only via the ratio ${s_k}/{b_k}$ in Eq. 11. From Fig. 3 in Ref. [1], we deduce that the two events detected by XENON1T correspond to {cS1,cS2b} coordinates {17,400} and {50,1300} from which we estimate $E_r = 12 , 33~\mathrm{keV}$ and $s_k/b_k=$ 0.7, 0.2 respectively.

Fig. 2

Left: the reconstructed acceptances $p_{eff}^0(E)$ , $p_{eff}^1(E)$ and $p_{eff}^2(E)$ compared to the XENON1T total ($p_{Xe}$) efficiency corresponding to SR1 from Fig. 1 in Ref. [1] (grey band). Right: the 90% excluded cross section for SI interactions obtained with $p_{eff}^0(E)$ (red-dash), $p_{eff}^1(E)$ (green-dash) and $p_{eff}^2(E)$ (blue-dash) as compared with XENON1T (black)

For an approximate recast of XENON1T for all masses, in particular for low DM masses and small recoil energies, we consider Eq. 12 as an equation for the efficiency p(E) that has to be satisfied for all masses in the interval [6–1000] GeV, for $\alpha =0.1$ and the cross sections $\sigma =\sigma ^{90}(M_{\chi })$ obtained by XENON1T for SI interactions.

At first approximation we neglect the last term in Eq. 11, thus assuming that there is some effective subspace in the S1/S2 parameter space where no events were detected and which can be used to reproduce the exclusion cross section.² This approximation is motivated by a comparison of the recoil energy distributions corresponding to DM masses of 35 GeV and 200 GeV shown in Fig. 1-right for $\sigma ^{90}(M_{\chi })$ of XENON1T and using $p_{Xe}$. Clearly these two signals practically coincide for low energies while the signal for 200 GeV becomes much larger for $E\gtrsim 8~\mathrm{keV}$. These two signals leading to the same level of exclusion might indicate that the events with large recoil energies do not contribute significantly to the 90% exclusion. It is indeed expected that for a low background experiment the region where events are found (here at energies above 12 keV) does not contribute significantly to the exclusion. Our approach can be considered as a simplified version of the XENON1T analysis where some cuts in cS1/cS2b space are applied to increase the signal/background ratio.

We denote $p_{eff}^0(E)$ the effective detector efficiency after all cuts assuming no events were detected. Equation 12 leads to an integral equation for $p_{eff}^0(E)$ for all masses in the range $6~\mathrm{GeV}<M_{\chi }<1000~ \mathrm{GeV}$,

$$\begin{aligned} {\mathcal {L}} \int p_{eff}^0(E)\frac{dN(M_{\chi },\sigma ^{90}(M_{\chi })}{dE} dE = -\log (\alpha )=log(10). \end{aligned}$$

(13)

Equation 13 is a Fredholm equations of the first kind. The solution of such equations is not stable and leads to large oscillations in $p_{eff}(E)$. To smooth out these oscillations, rather than solving it directly we minimize the functional

$$\begin{aligned} J(p_{eff})= & {} \max _{M_{\chi }}\left( \left| {\mathcal {L}} \int p_{eff}(E)\frac{dN(M_{\chi },\sigma ^{90})}{ dE} dE +log(\alpha ) \right| \right. \nonumber \\&+ \left. \kappa \int \left| \frac{d^2}{dE^2}p_{eff}(E)\right| ^2 dE \right) \end{aligned}$$

(14)

with respect to the function $p_{eff}(E)$. The minimization covers all masses in the interval considered, thus allowing to obtain a good agreement for each mass. Note that the term with $\kappa$ damps oscillations only if $\kappa$ is large enough, while it spoils the solution to Eq. 13 when $\kappa$ becomes too large. The goal is therefore to find the minimal $\kappa$ which allows to obtain a solution without oscillations. To find the minimum of $J(p_{eff})$ we tabulate $p_{eff}(E)$ on a grid which extends from $E_0$ to some $E_{max}$ with a 1keV step size. The acceptance $p_{eff}(E)$ vanishes for $E\le E_0$, where $E_0$, the detection threshold, is taken as a free parameter as well. It is found to be $E_0=1$keV. The values of the function $p_{eff}(E)$ at each point on the grid except the first one are also free parameters and we impose the condition that $p_{eff}(E)\ge 0$ . Between grid points we use a cubic polynomial interpolation.

The solution $p_{eff}^0$ is shown in Fig. 2-left and is compared to $p_{{Xe}}$ the XENON1T total efficiency in Ref. [1]. The total efficiency corresponds to the second science run (SR1). We find that $p_{eff}^0$ nearly vanishes at the recoil energy of the detected event with the smallest recoil energy. This condition was not imposed in advance and testifies of the validity of our assumption. Indeed if we had found a non-negligible $p_{eff}^0$ in the region where events are detected, we would not be able to conclude that the results of XENON1T which were obtained using the full events space and including all 735 observed events in the likelihood, can be reproduced in the zero-event approximation.³ The fact that $p_{eff}^0$ and $p_{Xe}$ are comparable (allowing for the uncertainty in the XENON1T efficiency) at low energies is also consistent with the observation that for light DM (say 10GeV) for which the signal is concentrated at low recoil energies, the signal is located in a region without electronic recoil background events as can be seen by comparing Fig. 8 in [40] and Fig. 3 in [1]. We observe also that close to threshold, the efficiency $p_{eff}^{0}$ is slightly larger than $p_{Xe}(E)$, this is probably related to the fact that in this region the efficiencies rise sharply hence have larger uncertainties.

Some comments are in order. First, it is well known that the Bayesian credible intervals significantly depend on priors in the no-event case, for example, with the flat prior used the 90% exclusion corresponds to $\approx 2.3$ signal events while for Jeffreys prior it corresponds to $\approx 1.3$ events. Changing the prior in Eq. 12 would therefore rescale the efficiency obtained, however there would be no impact on the excluded cross section as the same prior is used to fit the efficiency and to calculate the exclusion. Second, we chose Bayesian statistics over the frequentist approach adopted by XENON1T because we have no information about the background distribution. In case of zero event and background, the upper limit for exclusion depends only on the number of predicted signal events, this in turn depends on the choice of the statistical method. Thus all approaches, whether Bayesian or frequentist, will lead to the same reconstructed efficiency $p_{eff}^{0}$ up to an overall scaling factor. This means that in the framework of inverse recasting, one will reproduce the same result for the 90% upper limit on the DM scattering cross-section. Moreover in the low recoil energy interval (1–14 keV) we estimate the background to be low, only 0.35 events, assuming a simple counting experiment with Feldman–Cousins statistics we can estimate that the number of events required for a 90% exclusion with such background changes from 2.44 to 2.08, hence a 15% correction. Thus including background would only slightly modify the efficiency $p_{eff}^{0}$, without affecting the exclusion limit.

Following the same procedure, we also derive the efficiency $p_{eff}^1$ taking into account one detected event in the subsample at $E_r \approx 12 \mathrm{keV}$. The corresponding solution for the efficiency is shown in Fig. 2. At last we find the efficiency $p_{eff}^2$ which allows to reproduce the XENON1T exclusion curve using the extended optimum interval method by Yellin [41]. In this case we assume that the background is uniformly distributed in the interval [1–50] keV. Note that for the best-fit signal of XENON1T, $p_{eff}^2$ leads to 1.04 events, only slightly smaller than $n_R$ in Eq. 10.

The recasted 90% excluded cross sections obtained with $p_{eff}^{0,1,2}$ are displayed in Fig. 2-right. All allow to reproduce the XENON1T exclusion within 10%. The largest difference is found near $M_{\chi } \approx 35GeV$, see Fig. 2. We have compared the results for the three different recasts for the applications in Sect. 5 and found no significant difference. Note that we did some approximations, for instance, we use only a sub-space of the full analysis where background is small, we ignored background uncertainty as well as the uncertainty in the energy of detected events. However the error introduced by these approximations are compensated by the fact that we fit the exclusion curve when solving for $p_{eff}$. In Sect. 4, we will see that despite these approximations, our recast works well for the slightly different recoil spectra that are expected for SD interactions.

To conclude this section we emphasize that our method of reconstruction of the XENON1T efficiency neglects the contribution from large recoil energies. Thus this recast cannot be used for models which produce a recoil energy spectrum dominant at large energies such as can be obtained with effective operators or with inelastic scattering [33,34,35], it would lead to too conservative limits.

Fig. 3

Comparison between the 90% excluded SI cross section on protons from DS-50 with a binomial distribution (full black) and from micrOMEGAs with the default option, a Poisson distribution using all bins with $n_{e-} \ge 4$ (full red). The left panel also shows the exclusion for light masses when using a binomial distribution in micrOMEGAs (red dash). The right panel shows the difference in the exclusion from DS-50 with (full) and without (dash) the binomial distribution as well as the impact of using only the bins with $n_{e-} \ge 7$ within micrOMEGAs (red dash)

3.2 DarkSide-50

The DarkSide-50 (DS-50) collaboration [3] provides the basic experimental data to allow to reproduce the experimental results using a standard procedure. In particular, the distribution for the number of ionizations $n_{e^-}$ in the Argon detector for an exposure ${\mathcal {L}}=6786~\mathrm{kg\cdot days}$ together with an estimation of the background and the ionization quenching are given. We use the numerical tables for the data and background provided by the DarkSide collaboration. We are thus able to construct a likelihood based on the Poisson formula

$$\begin{aligned} L=\prod _i \frac{(B_i+S_i)^{n_i}}{n_i!} \exp ^{-(B_i+S_i)} \end{aligned}$$

(15)

where $B_i$ and $S_i$ are the number of background and signal events in the ith bin where the bins are defined for the distribution of the number of ionizations. Following DS-50 analysis, we do not include the bins $n_{e^-}<4$ in the likelihood. For $4\le n_{e^-}<7$ there is a large difference between the data and the estimated background, hence, following the DS-50 analysis, we treat the additional background as a nuisance parameter when constructing the likelihood function [42]. It means that we include the contribution of the $4\le n_{e^-}<7$ only if the DM signal plus known background is larger than the experimental data.

The average number of ionizations is determined by quenching. The ionization quenching depends on the recoil energy and suffers from a large uncertainty [43]. We use the minimal quenching.⁴ We have checked that making a linear interpolation between the minimal and maximal values of the quenching for each energy and treating the parameter of interpolation as a nuisance parameter leads to very similar results. Moreover, the distribution of the number of ionizations around the average is not known. The assumption made for describing this distribution is essential for light DM, since one can find events with $n_{e^-}\ge 4$ that arise from the tail of the distribution.

DS-50 considers two cases, first a binomial distribution for the number of ionizations where the average number of ionizations is determined by quenching, while the maximal number of ionizations is determined by the minimal energy needed for one ionization. The minimal energy is set to $E_1=19.5$ eV. Second a $\delta$-like distribution is considered where the actual number of events equals the average one. The key feature of the binomial distribution is that it allows to naturally implement an energy threshold which cannot be done with the widely used Poisson distribution. However the Poisson distribution can be generalized in order to take into account an energy threshold, thus we use the following distribution for the number of ionized electron:

$$\begin{aligned} p(n) = \frac{C(\bar{n},E_1)}{n!} \left( \bar{n}(1-n\frac{E_1}{E})\right) ^n \end{aligned}$$

(16)

Here C is defined by the normalization condition, E is the energy of an atom after DM recoils, $E_1$ is the minimal ionization energy. The value of $\bar{n}$ is chosen in such a way to reproduce the quenching given by DS-50 [3]. This generalized Poisson distribution has a tail that decreases faster than the binomial distribution and thus provides a more conservative limit. To compute the signal in each bin, $S_i$, and construct the likelihood in Eq. 15, we convolute the predicted recoil energy distribution with the generalized Poisson distribution for the number of ionized electron, Eq. 16, as well as with a gaussian with 20% resolution. The latter is used to describe the uncertainty introduced when DS-50 reconstructs the number of ionizations from their counts of photons.

As the DS-50 collaboration, we use the profile likelihood ratio, $\lambda$, to calculate the confidence level for excluding a model [42]

$$\begin{aligned} C.L. = 1-\frac{1}{2} \int \limits _{-2\log (\lambda )}^{\infty } \frac{dx}{\sqrt{2\pi x}} e^{-\frac{x}{2}} \end{aligned}$$

(17)

Following the procedure of DS-50, we have implemented a nuisance parameter for the background which is represented by an overall factor of $\pm 15\%$.⁵

For $M_{\chi }>1.8GeV$ both the binomial and improved Poisson distributions lead to the same exclusion. In Fig. 3 (right) we compare our reconstructions of the 90% excluded cross section with the DS-50 exclusion. Note that DS-50 uses two different likelihoods, one using bins $n_{e^-}\ge 4$ for masses $M_{\chi }<2.9$ GeV and one using only the bins $n_{e^-}\ge 7$ for higher masses. Rather than splitting our analysis for different mass range and in order to have a smooth exclusion, we take into account all bins $n_{e^-}\ge 4$ for the whole DM mass range. We still reproduce well the DS-50 exclusion for masses $M_{\chi }>3.5$ GeV, since in this region the contribution from higher bins dominate. Around $M_{\chi }\approx 3$ GeV our exclusion is stronger since the bins $4\le n_{e^-}< 7$ give an important contribution to the likelihood. Finally our exclusion is more conservative at lower masses, by about 50% (200%) for $M_{\chi }\approx 1.8 (0.65)$ GeV.

In Fig. 3-left we compare the exclusion cross section obtained by DS-50 assuming a binomial distribution with the ones reconstructed by micrOMEGAs for both the binomial and generalized Poisson distributions for $M_\chi <2$ GeV, here we use all bins $n_{e^-}\ge 4$. Note however that for masses below 1.8 GeV there are large uncertainties in the DS-50 exclusion depending on the choice of quenching model [43].

3.3 PICO-60

PICO [4, 5] is a Bubble Chamber experiment which uses $C_3F_8$, with 1167 kg-day exposure at a thermodynamic threshold of 3.3 keV and 1404 kg$\cdot$ days at 2.45 keV. After the acoustic parameter cut, PICO reports 3 candidate events for the second run while no events were detected in the first run [4]. A combined analysis of both runs which includes a new efficiency for the first run was published in Ref.  [5].

To reconstruct the PICO-60 exclusion curve for SI interactions, we assume the central value of the acceptance shown in Fig. 3 of Ref. [5] for each run. In our statistical analysis we compare the total number of expected events for both runs combined with the total number of detected events. We estimate the total background for both runs as $B=1.47$ events assuming that the ratio of single to multiple bubble events caused by neutrons is 1/4. We use two statistical methods based on Feldman–Cousins  [45] and Neyman with one-side belt with the confidence level,

$$\begin{aligned} C.L.= \sum _{n=0}^{3} \frac{(S+B)^n}{n!} e^{-(S+B)} \end{aligned}$$

(18)

where S is the number of predicted events caused by DM, and B is the expected background. In both cases we reproduce the 90% exclusion for the SI cross section of PICO-60 [5] within 10%, see Fig. 4. For the applications in the following sections we will use the recasting based on Feldman-Cousins. The main result of PICO-60 however concerns limits on the DM-proton SD cross section, this result will be discussed in Sect. 4.

Fig. 4

Comparison of PICO SI 90% excluded cross section with our reconstructions based on Feldman–Cousins (FC) and Neyman methods (N)

3.4 CRESST

The CRESST-III detector uses CaWO$_4$ and the limits obtained correspond to data collected with a total exposure of 5.594 kg$\cdot$days or 3.64 kg$\cdot$days after cuts [10]. In this experiment, the background is not estimated and the Optimum Interval method of Yellin  [46] is used to set a limit on the DM cross section for unknown background. With its low nucleus recoil threshold of 30.1 eV, the CRESST-III detector is sensitive to DM masses larger than 188 MeV assuming the standard parameters for the DM velocity distribution, Eq. 5. Moreover, DM masses as low as 160 MeV can be probed when taking into account energy resolution.

Using the Optimum Interval method [41] and the data presented in [47], we have recasted the exclusion limit of CRESST-III. We use the total exposure of 5.594 kg$\cdot$days and take into account the cut-survival probability and the acceptance for each nucleus shown in Figs. 4 and 6 of [47]. The energy resolution was considered as a free parameter which was fitted to get the best agreement with CRESST-III low masses exclusion. Namely we use a Gaussian with $\sigma =5.5$ eV and with a cut at $2\sigma$. In this manner the 90% exclusion cross section of CRESST-III is reproduced with 10% precision, see Fig. 5.

Fig. 5

Comparison of the reconstructed exclusion for SI interactions from micOMEGAs with CRESST-III

Fig. 6

Comparison of the recasted 90% limit on $\sigma ^\mathrm{SD}_{\chi n}$ (left) and $\sigma ^{SD}_{\chi p}$ (right ) from micrOMEGAs with the XENON1T limits [2] (black) with different choices of form factors: SHELL (green/dot), SHELL-min (blue/dash-dot) [38] and EFT [39] (red/dash)

Fig. 7

Left: comparison of the recasted 90% limit on $\sigma ^\mathrm{SD}_{\chi p}$ from micrOMEGAs (red) with the PICO-60 limit [4, 5] (black) using the EFT form factors from Ref. [39]. The impact of the choice of form factor is illustrated for the SHELL (green-dot) and SHELL-min form factors (blue-dot-dash). Right: comparison of the recasted 90% limit on $\sigma ^\mathrm{SD}_{\chi n}$ from micrOMEGAs (red) with the CRESST-III limit [47] (black) with zero momentum form factors

4 Spin-dependent interactions: recasting experimental exclusions

In general, SI and SD interactions on a given atom lead to very similar recoil energy spectra. Their difference is typically around 5% and is due only to the small momentum dependences of the SI and SD nucleus form factors. Thus experimentalists use the same set of cuts and the same background estimation for both SI and SD interactions. It is therefore justified to use the recasting done for SI interactions and apply it directly to SD interactions. Because there is a strong dependence on the SD form factors, to perform the recasting we use the same set of form factors as each experiment. These were obtained in [38] and  [39] and we cite them here as SHELL and EFT respectively. Moreover, for the first the authors derive the theoretical uncertainty, we also compare our results with those obtained with the minimal form factors leading to the more robust exclusion, we cite this minimal set as SHELL-min, see Appendix A.3. When we derive the 90% limit on SD cross sections using the same SD form factors used in each experiment, we find that our limit agrees with the experimental result with the same level of accuracy found for SI interactions as will be described below.

First we derive the 90% limit on SD cross sections on neutrons and protons for XENON1T, for this we take the SHELL SD form factors [38] which are also used by XENON1T. We find that an agreement below the 15% level with the limits on both $\sigma ^{SD}_{\chi n}$ and $\sigma ^{SD}_{\chi p}$, see Fig. 6. Taking into account the uncertainty on these form factors has little impact on $\sigma ^{SD}_{\chi n}$, but weakens the limit on $\sigma ^{SD}_{\chi p}$ by roughly a factor 2. The form factors EFT lead to a more stringent limit on $\sigma ^{SD}_{\chi n}$ while the limit on $\sigma ^{SD}_{\chi p}$ weakens by more than one order of magnitude. Note however that XENON1T has a much lower sensitivity to $\sigma ^{SD}_{\chi p}$, indeed Xenon has an even number of protons and their spins nearly cancel each other leading to small SD proton form factors.

Fig. 8

Left: comparison of the nucleus recoil distributions after folding in the XENON1T acceptance, $p_{eff}$, for a heavy mediator $M_M=100$ GeV (full) and $M_{\chi }=10,15$ GeV with that of a light mediator, $M_M=10$ MeV (dash) and $M_{\chi }=15$ GeV. Right: dependence of the excluded cross section on the mediator mass for $M_\chi =10,30,90$ GeV

Using the PICO acceptance described in Sect. 3.3 we derive the 90% limit on SD cross section on protons and compare it with the limit presented by the PICO collaboration [4, 5], see Fig. 7-left. For this we choose the form factors EFT also used by the experiment. Our reconstruction reproduces the PICO-60 exclusion within 10%, which is roughly the same precision that was obtained for SI interactions. To check the impact of the choice of form factors, we have also derived the exclusion using the SHELL and SHELL-min form factors. This weakens significantly the limit at low DM masses, up to a factor 2 at 4 GeV, while the effect is much more moderate for DM masses above 100 GeV. The difference with the EFT set remains below 10% (35%) for the SHELL (minimal) form factors.

CRESST-III is sensitive to spin-dependent DM-neutron interactions through the $^{17}O$ isotope despite its small abundance of 0.0367%. For this isotope, the SD form factor is only known in the zero momentum limit [37], we take the spin expectation $\langle S_n\rangle =0.5$. Following the same procedure as for SI interactions, we derive the recasted 90% limit on $\sigma ^\mathrm{SD}_{\chi n}$ and in Fig. 7-right, we make a comparison with the results of CRESST-III [47], the discrepancy is below $10\%$.⁶ The agreement with CRESST-III for the exclusion is at the same level as for the SI case.

5 Applications

In this section we show how to exploit our reconstruction of DD experimental limits to obtain limits on specific DM models while taking into account uncertainties from astrophysical and nuclear physics parameters. All numerical results presented below can be easily reproduced with the micrOMEGAs code. The corresponding code is stored in mdlIndep/dd_exp.c of micrOMEGAs.

5.1 The case of a light mediator

When DM-nucleus interactions are due to the exchange of a light mediator in t-channel, the standard formula that relates the DM-nucleon cross section at zero momentum with the recoil energy distribution cannot be applied. Indeed it rests on the assumption that the mass of the mediator is much larger than the Mandelstam variable $t=- 2 M_A E_R$ where $E_R$ is the nucleus recoil energy and $M_A$ the mass of the recoiling nucleus. For the typical minimal recoil energy $E_R\approx 2 \mathrm{keV}$ and $M_{Xe}$=130GeV this corresponds to $t=-(22\mathrm{MeV})^2$. Thus for mediator masses significantly below 1 GeV, an additional factor describing the t-dependence should be included. The recoil energy distribution from DM-nucleus elastic scattering is then replaced with

$$\begin{aligned} \frac{dN_{A}^{light}}{dE} = \frac{M_M^4}{(M_M^2 +2 M_{A} E)^2} \frac{dN_{A}^{std}(\sigma _0)}{dE} \end{aligned}$$

(19)

where $N_{A}^{std}$ is the standard expression for the number of recoil events for a point-like interaction, Eq. 1, [36] with elastic scattering cross section $\sigma _0$, $M_M$ is the mass of the t-channel mediator. Taking into account the contribution of the transfer momentum in the propagator of the light mediator leads to an overall decrease of the recoil signal and to a shift towards lower energies. This can be seen in Fig. 8 (left) where the signals for a DM with mass of 15 GeV are compared in the case of a light mediator $M_M=10~\mathrm{MeV}$ and a heavy mediator, $M_M=100~\mathrm{GeV}$. Moreover the recoil spectrum with the light mediator is shown to be very similar to the one for $M_\chi =10$ GeV and $M_M=100~\mathrm{GeV}$. These signals include the reconstructed acceptance of XENON1T, $p_{eff}$, and are obtained for $\sigma ^{90}$, to ease the comparison the distribution for the light mediator includes a normalisation factor.

In any model with a light mediator, we can use Eq. 19 to calculate the recoil energy signal and extract the dependence of the 90% excluded cross section on the mediator mass. The zero velocity excluded cross section ($\sigma _0$) is displayed in Fig. 8 (right) for XENON1T and for different DM masses. As expected, the mediator mass dependence comes into play at $M_M=100$MeV and the effect is significant at 50MeV. For very small mediator masses, all DM masses have a similar dependence on $M_M$, the reason is that the key ingredient in setting the limit is the detector threshold. The model independent limits on SI interactions in the case of a light mediator obtained from the micrOMEGAs recasting are compared in Fig. 9 for different experiments. Moreover the limit derived by the XENON1T collaboration using a S2 only analysis that allows to extend the sensitivity to lower masses is also shown for comparison [9].

Fig. 9

Limits on the spin-independent DM nucleon point-like cross section for a light mediator, $M_M=10$ MeV, using the micrOMEGAs recast of XENON1T, DarkSide-50, PICO-60 and CRESST-III. The limit derived by XENON1T using a ionization-only analysis is also displayed, XENON1T-S2 [9]

Fig. 10

Limits on the Z’ coupling, $g=g_f g_\chi$, from DarkSide-50, PICO-60, and XENON1T, for the Z’ model with pure vector couplings (left) and from PICO-60 and XENON1T for pure axial-vector couplings (right) for a light mediator, $M_Z'= 1$ MeV

To illustrate the effect of the light mediator on the direct detection exclusion in a specific model we consider the case of a Z’ mediator with a universal coupling to SM fermions,

$$\begin{aligned} {\mathcal {L}}= & {} - Z'_\mu \left( g_\chi \bar{\chi }\gamma ^\mu \chi + g'_\chi \bar{\chi }\gamma ^\mu \gamma ^5 \chi \right) \nonumber \\&- \sum _f Z'_\mu \left( g_f \bar{f} \gamma ^\mu f +g'_f \bar{f} \gamma ^\mu \gamma _5 f \right) \end{aligned}$$

(20)

We assume either pure vector couplings ($g'_\chi =g'_f=0$) or axial-vector ($g_\chi =g_f=0$) couplings which give rise respectively to SI and SD interactions. We further assume identical couplings to all fermions f. The results are displayed in Fig. 10 for both SI interactions and SD interactions. For SI interactions, the region that is compatible with the measured value of the relic density is excluded by XENON1T for $M_\chi > 8$ GeV, this region corresponds to $g=g_\chi g_f \approx 1.4 \times 10^{-12}$ with less than 10% variation over the mass range considered. For SD interactions, the current experiments cannot yet probe the preferred value for the relic density, as the couplings probed are roughly three orders of magnitude larger than the ones required by the relic density ($g \approx 1.4 \times 10^{-12}$). Moreover CRESST-III only probes values of couplings $g > 10^{-7}$, thus the corresponding limits are not displayed. Note that Fig. 10 shows the best limit whether it comes from SD interactions with protons (PICO) or neutrons (XENON1T).

5.2 Millicharged Dark matter

Millicharged DM which interacts with the SM through photons provides another example of a light mediator, the massless photon in this case. Typically a kinetic mixing between a new gauge boson and the hypercharge leads to DM interacting with the photon with a millicharge, $q_{\chi }$, [48]

$$\begin{aligned} {\mathcal {L}}= q_{\chi } e \bar{\chi }\gamma ^\mu \chi A_\mu \end{aligned}$$

(21)

where we have omitted the terms that describe interactions with the new gauge boson. The recoil energy distribution for DM nucleus elastic scattering is similar to the one for the light mediator, Eq. 19,

$$\begin{aligned} \frac{dN_{A}^{m}}{dE} = \frac{M_{ph}^4}{(2 M_{A} E)^2} \frac{dN_{A}^{std}(\sigma _0)}{dE} \end{aligned}$$

(22)

where

$$\begin{aligned} \sigma _0= 16 \pi \alpha ^2_{EM} q^2_\chi Z_A^2 \frac{\mu ^2_{\chi p}}{M^4_{ph}} \end{aligned}$$

(23)

and $M_{ph}$ is a parameter with mass dimension which does not enter the final result.

The 90% lower limits on $q_{\chi }$ obtained after imposing the DarkSide-50 and XENON1T limits are presented in Fig. 11.

Fig. 11

The 90% exclusion on the DM millicharge $q_{\chi }$ as a function of the DM mass using recasted results of DarkSide-50 (left) and of XENON1T (right). The region above the top curve cannot be probed by underground DD experiments

Direct detection experiments cannot probe large values of the charge, $q_\chi$, since a millicharged DM will loose energy through its interaction with rocks before it reaches the detector. Elastic scattering of DM particles with atomic nuclei is the main process responsible for energy loss. The cross section for elastic scattering reads [49]

$$\begin{aligned} \frac{d\sigma }{d\cos {\theta _{cm}}}=2\pi \left| \frac{2\mu }{q} \int \limits _0^\infty V(r)r \sin (rq)dr\right| ^2 \end{aligned}$$

(24)

where $\mu$ is the reduced mass of colliding particles and q is the transfer momentum. For a nucleus charge screened by electrons, the potential is given by

$$\begin{aligned} V(r)=\frac{q_\chi Z_A e^2}{r}e^{-r/R_A} \end{aligned}$$

(25)

where the atomic radius is approximated by $R_A \approx 0.8853 Z^{\frac{1}{3}}_A/m_e \alpha _{em}$. Note that this rough approximation is sufficient since the energy loss depends only logarithmically on $R_A$. The energy loss, $E_{loss}$ of one millicharged particle in an elastic collision with a nucleus is obtained after integrating Eq. 24,

$$\begin{aligned} \langle E_{lost}\sigma _A \rangle= & {} 2\pi \frac{(q_\chi Z_A e^2)^2}{v^2 M_A} \left( \log (1+(2v\mu _{\chi A} R_A)^2) \right. \nonumber \\&-\left. \frac{(2v\mu _{\chi A} R_A)^2}{1+(2v\mu _{\chi A} R_A)^2}\right) \end{aligned}$$

(26)

The energy, $E_\chi$ of a DM particle passing through the Earth is then given by

$$\begin{aligned} \frac{dE_\chi }{dx} = -\sum _{A \in Earth} \langle E_{lost}\sigma _A \rangle n_A \end{aligned}$$

(27)

where $n_A$ is the number density of the element A in the Earth and x is the distance from the surface. If the DM mass is above the GeV scale then $2v\mu _{\chi A} R_A \gg 1$ and $\frac{dE_\chi }{dx}\approx \frac{C}{E_\chi }$. Thus, at some finite distance from the Earth surface, the DM will stop and drift towards the Center of the Earth driven by gravitational interactions. For DM to be detectable its energy must be above the detector threshold $E_{tr}$, thus the condition that a DM with maximum velocity $v=v_{esc}+v_{Earth}$ will reach the detector located at a distance H below the surface of the Earth with $E_\chi >E_{tr}$ leads to a linear equation in $q_\chi ^2$

$$\begin{aligned}&\int \limits ^{E_{min}}_{E_{max}} \frac{dE_\chi }{{dE_\chi }/dx} = H \quad \mathrm{where}\quad E_{min} = \frac{E_{tr}M_A M_{\chi }}{4 \mu _{\chi A}^2} \nonumber \\&E_{max} = \frac{M_{\chi }}{2}(v_{Esc}+v_{Earth})^2 \end{aligned}$$

(28)

The corresponding upper limit on the millicharge excluded by either DarkSide-50 or XENON1T is at least three orders of magnitude above the respective lower limits, see Fig. 11. Here we used H = 1400 m and $E_{tr}= 0.1(1.6)$ keV for DarkSide-50 (XENON1T).

Fig. 12

Influence of the uncertainty on the DM velocity distribution on XENON1T (left) and DarkSide-50 (right) 90% excluded SI cross sections. The allowed region for the Maxwell distribution (shaded black) and the SHM++ distribution (shaded red) together with the full curves showing the standard central values for both distributions

Fig. 13

Impact of the velocity distribution on exclusion for the Z’ model with vector couplings from DarkSide-50 (left) and XENON1T (right) in the $M_{Z'}-M_\chi$ plane. Here $g_{Z'}=1\times 10^{-7}$, $g_\chi =5.5\times 10^{-5}$. The exclusion obtained with the Maxwell distribution and standard parameters is compared with the ones obtained with SHM++ with standard values for the parameters and with the parameters leading to the most severe constraint (SHM++(max)). The regions below the curves are excluded

5.3 Dependence on DM distributions

As mentioned previously, most experiments publish their results assuming that the DM velocity distribution in the neighbourhood of the Sun is a Maxwell distribution with parameters given in Eq. 5. However, the recent estimates for $\rho _{\chi }$ point to a slightly larger value [50, 51]

$$\begin{aligned} \rho _{\chi }= (0.39 \pm 0.03)(1.2 \pm 0.2)(1\pm \delta _{triax})\mathrm{GeV/cm^3} \end{aligned}$$

(29)

where $\delta _{triax}< 0.2$. Clearly, since $\rho _{\chi }$ is just an overall factor, changing its value will amount to simply rescaling the 90% excluded cross section by a factor of $\rho _{\chi }/0.3$.

To estimate the impact of the DD limits on the parameters of the velocity distribution we have varied the parameters of the Maxwell distribution within the range [52,53,54]

$$\begin{aligned} v_{Rot}= & {} 220\pm 18 \frac{km}{s},\quad v_{Earth}=232-252 \frac{km}{s} \nonumber \\ v_{Esc}= & {} 580\pm 63 \frac{km}{s},\quad \rho _\chi =0.468\pm 0.202 \end{aligned}$$

(30)

The strongest and weakest 90% excluded cross sections for the XENON1T experiment for these intervals are shown in Fig. 12 together with the exclusion corresponding to the standard parameters in Eq. 5. For DM masses above roughly 10 GeV, most of the variations in the exclusion limit is due to $\rho _{\chi }$, in particular the upper $1\sigma$ range leads to a more aggressive limit by about a factor 2 while the limit is weakened by around 10% when using the lowest value for $\rho _{\chi }$. Roughly another 10% shift in the limit is due to the variation of other parameters. For low DM masses, corrections can be much larger. For example for $M_{\chi }\approx 6~\mathrm{GeV}$ the excluded cross section increases by more than a factor 2. This is mainly due to a decrease in $v_{Esc}$ which requires a heavier DM to pass the threshold for nuclear recoils. For the same reason an increase in $v_{Esc}$ leads to a more aggressive limit. An alternative DM distribution which is compatible with Gaia data was suggested in Ref. [55], it leads to more stringent limits at all masses since the main difference with the Maxwell distribution is in the much larger central value for $\rho _{\chi }$. We have also varied the parameters of the SHM++ distribution within their $1\sigma$ range defined in Eq. 36 and found a near overlap of the most stringent exclusion with that of the Maxwell distribution. Again, $\rho _\chi$ and $v_{esc}$ are the parameters that have the largest impact on the exclusion limit. After factoring out the linear dependence on $\rho _\chi$, we still find that the limit shifts by more than a factor 2 for $M_{\chi }\approx 7~\mathrm{GeV}$ and by about 20% for $M_{\chi }> 200~\mathrm{GeV}$. Similar conclusions are obtained for the DarkSide-50 exclusions in the low mass region, see Fig. 13-right.

We also examine the impact of the velocity distribution on the exclusion limit for the simplified Z’ model with vector couplings introduced in the previous section. Fixing the values of the couplings to $g_{Z'}=1\times 10^{-7}$ and $g_\chi =5.5\times 10^{-5}$ we show how much the exclusion on the Z’ mass from DarkSide-50 and XENON1T can be reinforced assuming an aggressive exclusion with the SHM++ distribution. The latter, labelled SHM++(max) corresponds to the upper value of the $1\sigma$ range for the parameters $\rho _\chi ,{\text{ vrot }}, {\text{ vesc }}$ in Eq. 36. With this choice the lowest limit on $M_{Z'}$ increases by more than a factor 2 for $M_\chi =1.8$ GeV to about 40% for $M_{\chi }> 100$ GeV as compared to the Maxwell distribution with standard parameters, Eq. 5. This confirms our expectations that the impact of the velocity distribution is more important for spectra peaked at low energies.

6 Conclusion

In this paper we demonstrate how the results from recent DM direct detection experiments can be applied to DM models with features that can somewhat differ from the ones assumed when deriving the experimental limits. After validating the recast of experimental exclusions, we illustrated how these can be applied to specific DM models, in particular models with a light mediator or a millicharged DM for which the spectrum of nuclear energy recoil is shifted towards low energies from the one of a heavy mediator. We also illustrated the impact of the choice of nuclear form factor for spin dependent interactions and of the choice of velocity distributions. These recasts can also be used to derive direct detection limits on multicomponent DM. These recasts are available in micrOMEGAs which contains new routines that provide the exclusion cross section for the direct detection experiments that provide the best exclusion for spin independent and spin dependent interactions for DM masses from 160 MeV upto the TeV range. Note that these recasts can be used with any of the generic models implemented in micrOMEGAs and that as for all direct detection routines apply to models where the direct detection cross section can be described by the low-energy Lagrangians for fermion, scalar or vector DM listed in Ref. [27] extended to the case of light mediators. These routines will be extended to include future experimental limits as they become available.

Appendix A: micrOMEGAs routines

We describe the micrOMEGAs routines that can be used to extract constraints on DM models based on the results of the direct detection experiments. Examples on how to use these routines can be found in mdlIndep/dd_exp.c of micrOMEGAs. All results presented here can also be reproduced with this code.

1.1 A.1: Experimental data

The SI 90% DD limits tabulated from the results presented by XENON1T [1], DarkSide-50 [3], PICO-60 [5] and CRESST-III [10] are accessible through the following functions

• XENON1T_90(Mdm) for $6< M_{DM} < 1000$ GeV [1],

• DS50_90(Mdm) for $0.7< M_{DM} < 15$ GeV [3],

• PICO60_90(Mdm) for $3< M_{DM} < 10000$ GeV [5],

• CRESST_III_90(Mdm) for $0.35< M_{DM} < 12$ GeV [10].

The corresponding SD 90% exclusion limits are contained in the functions

• PICO60_SDp_90(Mdm) for $3< M_{DM} < 10000$ GeV [5],

• XENON1T_SDp_90(Mdm) for $6< M_{DM} < 1000$ GeV [2],

• XENON1T_SDn_90(Mdm) for $6< M_{DM} < 1000$ GeV [2],

• CRESST_III_SDn_90(Mdm) for $0.35< M_{DM} < 12$ GeV [10].

These functions give the excluded cross sections in $\mathrm{cm}^2$. For a DM mass outside the range specified the function returns NaN.

1.2 A.2: Recasting the experimental limits with micrOMEGAs

• DD_pvalCS(expCode, $f_v$, $\sigma _{SI_P},\sigma _{SI_N},\sigma _{SD_P}, \sigma _{SD_N}$,&expName)

calculates the value $\alpha = 1-C.L.$ for a model with DM-nucleon cross sections $\sigma _{SI_P},\sigma _{SI_N},\sigma _{SD_P}, \sigma _{SD_N}$. Cross sections are specified in [pb] units. The return value 0.1 corresponds to a 90% exclusion. The expCode parameter can be any of the codes XENON1T_2018,DarkSide_2018, CRESST_2019,PICO_2019 or their combination concatenated with the symbol $\mid$. There is also a predefined parameter that currently combines these experiments

The parameter char* expName is used to indicate the experiment that provides the best exclusion among those specified in expCode. The function DD_pvalCS calculates the exclusion for each experiment independently, returns the smallest $\alpha$, and assigns the name of the corresponding experiment to expName if it is not NULL.

The $f_v$ parameter specifies the DM velocity distribution in the detector frame. For example, one can use Maxwell or SHMpp which are included in micrOMEGAs, otherwise the user can define another distribution. The DM velocity distribution has to be normalized as in Eq. 4. The units are km/s for v and s/km for $f_v(v)$. DD_pvalCS implicitly depends on the global parameters Mcdm and rhoDM which specify the DM mass and DM local density respectively.

For XENON1T one can chose between $p_{eff}^{q}$ with $q=0,1,2$, see Sect. 3.1. The flag Xe1TnEvents=q allows to choose the corresponding recasting, otherwise and by default the code uses $p_{eff}^1$. For PICO-60, the user can choose between the recasting based on Feldman-Cousins statistics, PICO60Flag=0 which is the default value, or the one based on Neyman one side belt exclusion, PICO60Flag=1.

• DD_factorCS(expCode, $\alpha$, $f_v$, $\sigma _{SI_P},\sigma _{SI_N},\sigma _{SD_P},\sigma _{SD_N}$,&expName)

returns the overall factor which should be applied to the cross sections, $\sigma _{SI_P},\sigma _{SI_N},\sigma _{SD_P}, \sigma _{SD_N}$ to reach the exclusion level $\alpha$. All parameters are the same as in DD_pvalCS above.

• *dNdEFact(Enr_kev, A)

is the address of the function which modifies the nucleus recoil distribution for DD_pvalCS and DD_factorCS to take into account a t-channel propagator with small or zero mass. By default dNdEfact=NULL and this function does not contribute to the calculation of the direct detection cross sections. Otherwise it is taken as an additional factor in the nucleus recoil distribution, see Eq. 19. The parameter Enr_kev is the recoil energy in [keV] units, A is the atomic mass number of the nucleus. This function should be defined by the user, an example is given in mdlIndep/dd_exp.c.

• DD_pval(expCode, $f_v$,&expName)

• DD_factor(expCode, $\alpha$, $f_v$,&expName)

These functions are similar to DD_pvalCS and DD_factorCS described above but use the cross section calculated from the DM model under consideration in micrOMEGAs. The necessary corrections for a light mediator are implemented automatically, these functions do not use dNdEFact.

The routines described above require SD form factors when considering SD limits, by default they use the same form factors as each experiment. The SD form factors can be replaced using the command

• setSpinDepFF(ExperimentID, setID)

where the choice for ExperimentID is given at the beginning of this section, and setID can be

• EFT: corresponding to the form factors in [39], Eq. 8.

• SHELL: corresponding to the average form factors in [38], Eq. 31.

• SHELLm: corresponding to the minimal form factor of [38], Eq. 32.

See below.

1.3 A.3: Spin-dependent form factors

The minimal and maximal values for the SD form factors, $S_{00}(q),S_{01}(q),S_{11}(q)$ , are computed in Ref. [38] within the shell model. The flag spinDepFlag=SHELL corresponds to the average

$$\begin{aligned} S_{ab}=(S_{ab}^{min} +S_{ab}^{max})/2 \end{aligned}$$

(31)

which are obtained from the minimum and maximum fitted values in Table VI in [38].

Note that $S_{ab}^{min}$ often lead to a negative value for the subdominant component to the form factor. Since this has no physical meaning, to define the form factors that lead to the most robust exclusion we rather use the minimum value of the proton-only, $S_p^{min}$, and neutron-only, $S_n^{min}$, form factors also given in [38]. These correspond to the minimal form factor for the case when only one type of interaction (with proton or neutron) is included. With this we construct the nucleus form factors

$$\begin{aligned} S_{00}= & {} \frac{1}{4} \left( S_p^{min} +S_n^{min} \pm 2 \sqrt{S_p^{min} S_n^{min}}\right) \nonumber \\ S_{11}= & {} \frac{1}{4} \left( S_p^{min} +S_n^{min} {\mp } 2 \sqrt{S_p^{min} S_n^{min}}\right) \nonumber \\ S_{01}= & {} \frac{1}{2} (S_p^{min}-S_n^{min}) \end{aligned}$$

(32)

The sign in Eq. 32 is chosen to reproduce the ratio $S_{00}(0)/{S}_{11}(0)$ for the central value of the form factors in Ref. [38].

1.4 A.4: Velocity distribution

Ignoring the direction of motion of DM particles and the small effect of DM acceleration by the gravitational field of the Sun, the DM velocity distribution in the vicinity of the direct detection experiment is given by

$$\begin{aligned} f(\text{ v})= \int \limits _{|\vec {v}|<\mathrm{vEsc}} d^3\vec {v} F_G(\vec {v} -\vec {v}_{Earth}) \delta (\text{ v } -|\vec {v}|) \end{aligned}$$

where $F_G$ is the DM velocity distribution in the frame, of the galaxy, $\vec {v}_{Earth}$ is the velocity of the Earth in the Galaxy and vEsc is the maximal velocity in our Galaxy due to its finite gravitationnal potential. vEsc and vEarth=$|\vec {v}_{Earth}|$ are global parameters of micrOMEGAs.

The velocity distributions that are available in micrOMEGAs  are the following

• Maxwell(v)

returns

$$\begin{aligned} F_G^M(\text{ v})= c_{\mathrm{norm}} \frac{1}{(2\pi \mathrm{vRot}^2)^{3/2}} \exp \left( -\frac{(\vec {v})^2}{\mathrm{vRot}^2}\right) \theta (\text{ vEsc } -|\vec {v}|) \end{aligned}$$

which corresponds to the isothermal model. Here vRot is the orbital velocity of stars in the Milky Way, it is also a global parameter of micrOMEGAs. $c_\mathtt{norm}$ is the normalization factor,

$$\begin{aligned} c_{\mathtt{norm}}^{-1}= {\text{ erf }}\left( \frac{\text{ vEsc }}{\mathrm{vRot}} \right) -\frac{2}{\sqrt{\pi }} \frac{\text{ vEsc }}{\mathrm{vRot}} \exp \left( -\frac{\text{ vEsc}^2}{\mathrm{vRot}^2}\right) \end{aligned}$$

• SHMpp(v)

returns the velocity distribution SHM++ proposed in [55].

$$\begin{aligned} F_G(\vec {v}) = (1-\eta ) F_G^M(\text{ v}) +\eta F_G^S(\text{ v}) \end{aligned}$$

(33)

This distribution consists of two components. The first, $F_G^M(\vec {v})$, is the standard Maxwell velocity distribution described above. The second component is the velocity distribution from the Gaia sausage [56, 57], it is not spherically symmetric and is defined by the anisotropy parameter $\beta$ with

$$\begin{aligned} F_G^S(\vec {v})= & {} \frac{c_{\mathrm{norm}} }{(2\pi )^{3/2} \Delta v_r \Delta v_\theta \Delta v_\phi } \exp \left( -\left( \frac{v_r}{\Delta v_r}\right) ^2 -\left( \frac{v_\theta }{\Delta v_\theta }\right) ^2\right. \nonumber \\&-\left. \left( \frac{v_\phi }{\Delta v_\phi }\right) ^2 \right) \theta (\text{ vEsc } -|\vec {v}|) \end{aligned}$$

(34)

where

$$\begin{aligned} \Delta v_r= \frac{\mathrm{vRot}}{\sqrt{1-\frac{2}{3}\beta }}\;,\;\; \Delta v_\phi =\Delta v_\theta =\frac{\mathrm{vRot} \sqrt{1-\beta }}{\sqrt{1-\frac{2}{3}\beta }} \end{aligned}$$

(35)

and

$$\begin{aligned} c_{\mathtt{norm}}^{-1}= & {} {\text{ erf }}\left( \frac{\text{ vEsc }}{\mathrm{vRot}} \right) - \left( \frac{1-\beta }{\beta } \right) ^{1/2}\\&\times \exp \left( -\frac{\text{ vEsc}^2}{\mathrm{vRot}^2}\right) {\text{ erfi }}\left( \frac{\text{ vEsc }}{\mathrm{vRot}}\frac{\beta ^{1/2}}{(1-\beta )^{1/2}} \right) \end{aligned}$$

where ${\text{ erfi }}$ is the imaginary error function.

The central values and uncertainties of the SHM++ parameters are

$$\begin{aligned} \text {rhoDM}= & {} 0.55\pm 0.17~\mathrm{GeV}/\mathrm{cm}^3\nonumber \\ \mathrm{vRot}= & {} 233\pm 3~\mathrm{km/s}\nonumber \\ \text{ vEsc }= & {} 580\pm 63~ \mathrm{km/s}\nonumber \\ \beta =\text {betaSHMpp}= & {} 0.9\pm 0.05\nonumber \\ \eta = \text {etaSHMpp}= & {} 0.2\pm 0.1 \end{aligned}$$

(36)

Note that these central values for the global parameters, vRot, vEsc and rhoDM are different from the ones in Eq. 5.