Inelastic dark matter nucleus scattering

TL;DR

This paper explores how inelastic dark matter nucleus scattering can improve detection prospects and particle physics insights in xenon-based detectors, beyond traditional elastic scattering methods.

Contribution

It introduces a low-energy effective field theory framework to analyze inelastic scattering signals in dark matter detection experiments.

Findings

Inelastic transitions increase the discovery potential once a dark matter signal is detected.

Inelastic scattering helps better determine the underlying particle physics.

Inelastic signals provide additional detectable signatures in xenon detectors.

Abstract

Direct detection experiments aim at the detection of dark matter in the form of weakly interacting massive particles (WIMPs) by searching for signals from elastic dark matter nucleus scattering. Additionally, inelastic scattering in which the nucleus is excited is expected from nuclear physics and provides an additional detectable signal. In the context of a low-energy effective field theory we investigate the experimental reach to these inelastic transitions for xenon-based detectors employing a dual-phase time projection chamber. We find that once a dark matter signal is established, inelastic transitions enhance the discovery reach and we show that they allow a better determination of the underlying particle physics.

Tables

Table 1. Table 1 : Non-relativistic effective operators for DM-nucleon interactions. We follow the numbering scheme of  [ 31 ] .

${𝒪}_1=1_{\chi }{1}_N$	${𝒪}_3=i{\overrightarrow{S}}_N\cdot \left(\frac{\overrightarrow{q}}{m_N}\times {\overrightarrow{v}}^{⟂}\right)1_{\chi }$
${𝒪}_4={\overrightarrow{S}}_{\chi }\cdot {\overrightarrow{S}}_N$	${𝒪}_5=i{\overrightarrow{S}}_{\chi }\cdot \left(\frac{\overrightarrow{q}}{m_N}\times {\overrightarrow{v}}^{⟂}\right)1_N$
${𝒪}_6=\left({\overrightarrow{S}}_{\chi }\cdot \frac{\overrightarrow{q}}{m_N}\right)\left({\overrightarrow{S}}_N\cdot \frac{q}{m_N}\right)$	${𝒪}_7={\overrightarrow{S}}_N\cdot {\overrightarrow{v}}^{⟂}{1}_{\chi }$
${𝒪}_8={\overrightarrow{S}}_{\chi }\cdot {\overrightarrow{v}}^{⟂}{1}_N$	${𝒪}_9=i{\overrightarrow{S}}_{\chi }\cdot \left({\overrightarrow{S}}_N\times \frac{\overrightarrow{q}}{m_N}\right)$
${𝒪}_{10}=i{\overrightarrow{S}}_N\cdot \frac{\overrightarrow{q}}{m_N}{1}_{\chi }$	${𝒪}_{11}=i{\overrightarrow{S}}_{\chi }\cdot \frac{\overrightarrow{q}}{m_N}{1}_N$
${𝒪}_{12}={\overrightarrow{S}}_{\chi }\cdot \left({\overrightarrow{S}}_N\times {\overrightarrow{v}}^{⟂}\right)$	${𝒪}_{13}=i\left({\overrightarrow{S}}_{\chi }\cdot {\overrightarrow{v}}^{⟂}\right)\left({\overrightarrow{S}}_N\cdot \frac{\overrightarrow{q}}{m_N}\right)$
${𝒪}_{14}=i\left({\overrightarrow{S}}_{\chi }\cdot \frac{\overrightarrow{q}}{m_N}\right)\left({\overrightarrow{S}}_N\cdot {\overrightarrow{v}}^{⟂}\right)$	${𝒪}_{15}=-\left({\overrightarrow{S}}_{\chi }\cdot \frac{\overrightarrow{q}}{m_N}\right)\left(\left({\overrightarrow{S}}_N\times {\overrightarrow{v}}^{⟂}\right)\cdot \frac{\overrightarrow{q}}{m_N}\right)$

Table 2. Table 2 : Ratio of the inelastic scattering rate on natural xenon R i ​ n subscript 𝑅 𝑖 𝑛 R_{in} relative the elastic rate R e ​ l subscript 𝑅 𝑒 𝑙 R_{el} for the 14 Operators of the non-relativistic effective theory for m χ = 300 subscript 𝑚 𝜒 300 m_{\chi}=300\, GeV and c k p = c k n subscript superscript 𝑐 𝑝 𝑘 subscript superscript 𝑐 𝑛 𝑘 c^{p}_{k}=c^{n}_{k} . For operators marked by " − - " the ratio is < 10 − 3 absent superscript 10 3 <10^{-3} .

Operator	1	3	4	5	6	7	8	9	10	11	12	13	14	15
$R_{in}/R_{el}$	$-$	$-$	0.13	0.04	0.06	62.	0.009	0.07	0.11	$-$	$-$	27.	$-$	$-$

Table 3. Table 3 : Overview on detector specific parameters in the most sensitive dual-phase LXe detectors to date.

Parameter	XENON1T	LUX	PandaX-II
Target mass [kg]	2000 [2]	250 [3]	580 [4]
Exposure [t$\cdot$yr]	1 [2]	0.09 [3]	0.15 [4]
$E$ [V/cm]	81 [2]	181 [3]	400[4]
$g_1$ [PE/$\gamma$]	$\sim 0.14$ [53]	$\sim 0.12$ [54]	$\sim 0.11$ [4]
$\omega$ [PE/e-]	$\sim 28$ [53]	$\sim 25$ [54]	$\sim 24$ [4]
${\sigma }_{\omega }$ [PE/e-]	$\sim 7$ [53]	$\sim 6$ [54]	$\sim 7$ [55]
$p_{\mathrm{extract}}$ [%]	$\sim 93$ [53]	$\sim 50$ [54]	$\sim 55$ [4]
${\sigma }_{\mathrm{PMT}}$ [PE/$\gamma$]	$\sim 0.4$ [56]	$\sim 0.3$ [57]	$\sim 0.4$ [56]

Table 4. Table 4 : List of parameters used in the simulation of liquid xenon microphysics. See text for further explanations.

Description	Parameter	ER	NR	Reference
Lindhard quenching factor	$L$	1	$\frac{k\cdot g\left(\kappa \right)}{1+k\cdot g\left(\kappa \right)}$	[44]
Fano factor	$F$	0.03	0.03	[45]
Exciton-ion-ratio	$n_{\mathrm{ex}}/n_i$	0.15	$1.240\cdot E^{-0.0472}\cdot \left(1-e^{-239\kappa }\right)$	[44]
Recombination probability	$r$	$\frac{\left(n_{\gamma }/n_e\right)-n_{\mathrm{ex}}/n_e}{\left(n_{\gamma }/n_e\right)+1}$	$1-\frac{\ln \left(1+n_i\xi \right)}{n_i\xi }$	[46, 43]
Recombination fluctuation	${\sigma }_r$	0.06	$n_i\sqrt{0.0056\left(1-r\right)}$	[44]
Penning quenching	$f_{\text{l}}$	1	$\frac{1}{1+\eta {\kappa }^{\lambda }}$	[44]

Equations

$$\mathcal{H}=\sum _{\tau =0,1}\sum _{k=1}^{15}{c}_k^{\tau }{\mathcal{O}}_k{t}^{\tau },$$

$$\frac{\text{mbox}d{\sigma }_T(v^2,E_R)}{\text{mbox}d{E}_R}=\frac{m_T}{2\pi v^2}⟨|{\mathcal{M}}_{\mathrm{eff}}{|}^2⟩$$

$$⟨|{\mathcal{M}}_{\mathrm{eff}}{|}^2⟩$$

$$+R_{{\Sigma }^{{}^{\prime \prime }}}^{\tau {\tau }^{{}^{\prime }}}⟨j_i||{\Sigma }_{L;\tau }^{{}^{\prime \prime }}||j_f⟩⟨j_f||{\Sigma }_{L;{\tau }^{{}^{\prime }}}^{{}^{\prime \prime }}||j_i⟩+R_{{\Sigma }^{{}^{\prime }}}^{\tau {\tau }^{{}^{\prime }}}⟨j_i||{\Sigma }_{L;\tau }^{{}^{\prime }}|j_f⟩⟨j_f||{\Sigma }_{L;{\tau }^{{}^{\prime }}}^{{}^{\prime }}||j_i⟩$$

$$+\frac{q^2}{m_N^2}{R}_{{\Phi }^{{}^{\prime \prime }}}^{\tau {\tau }^{{}^{\prime }}}⟨j_i||{\Phi }_{L;\tau }^{{}^{\prime \prime }}||j_f⟩⟨j_f||{\Phi }_{L;{\tau }^{{}^{\prime }}}^{{}^{\prime \prime }}||j_i⟩+2\frac{q^2}{m_N^2}{R}_{M{\Phi }^{{}^{\prime \prime }}}^{\tau {\tau }^{{}^{\prime }}}⟨j_i||{\Phi }_{L;\tau }^{{}^{\prime \prime }}||j_f⟩⟨j_f||M_{L;{\tau }^{{}^{\prime }}}||j_i⟩$$

$$+\frac{q^2}{m_N^2}{R}_{{\tilde{\Phi }}^{{}^{\prime }}}^{\tau {\tau }^{{}^{\prime }}}⟨j_i||{\tilde{\Phi }}_{L;\tau }^{{}^{\prime }}||j_f⟩⟨j_f||{\tilde{\Phi }}_{L;{\tau }^{{}^{\prime }}}^{{}^{\prime }}||j_i⟩+\frac{q^2}{m_N^2}{R}_{\Delta }^{\tau {\tau }^{{}^{\prime }}}⟨j_i||{\Delta }_{L;\tau }||j_f⟩⟨j_f||{\Delta }_{L;{\tau }^{{}^{\prime }}}||j_i⟩$$

$$+\frac{q^2}{m_N^2}{R}_{\Delta {\Sigma }^{{}^{\prime }}}^{\tau {\tau }^{{}^{\prime }}}⟨j_i||{\Sigma }_{L;\tau }^{{}^{\prime }}||j_f⟩⟨j_i||{\Delta }_{L;{\tau }^{{}^{\prime }}}||j_i⟩$$

$$+\frac{q^2}{m_N^2}{R}_{\Sigma }^{\tau {\tau }^{{}^{\prime }}}⟨j_i||{\Sigma }_{L;\tau }||j_f⟩⟨j_f||{\Sigma }_{L;{\tau }^{{}^{\prime }}}||j_i⟩+\frac{q^2}{m_N^2}{R}_{{\Delta }^{{}^{\prime }}}^{\tau {\tau }^{{}^{\prime }}}⟨j_i||{\Delta }_{L;\tau }^{{}^{\prime }}||j_f⟩⟨j_f||{\Delta }_{L;{\tau }^{{}^{\prime }}}^{{}^{\prime }}||j_i⟩$$

$$+\frac{q^2}{m_N^2}{R}_{\tilde{\Omega }}^{\tau {\tau }^{{}^{\prime }}}⟨j_i||{\tilde{\Omega }}_{L;\tau }|j_f⟩⟨j_f||{\tilde{\Omega }}_{L;{\tau }^{\prime }}|j_i⟩+\frac{q^2}{m_N^2}{R}_{{\tilde{\Phi }}^{{}^{\prime }}}^{\tau {\tau }^{{}^{\prime }}}⟨j_i||{\tilde{\Phi }}_{L;\tau }||j_f⟩⟨j_f||{\tilde{\Phi }}_{L;{\tau }^{{}^{\prime }}}||j_i⟩$$

$$+\frac{q^2}{m_N^2}{R}_{{\tilde{\Delta }}^{{}^{\prime \prime }}}^{\tau {\tau }^{{}^{\prime }}}⟨j_i||{\tilde{\Delta }}_{L;\tau }^{{}^{\prime \prime }}||j_f⟩⟨j_f||{\tilde{\Delta }}_{L;{\tau }^{{}^{\prime }}}^{{}^{\prime \prime }}||j_i⟩+\frac{q^2}{m_N^2}{R}_{{\Delta }^{{}^{\prime }}\Sigma }^{\tau {\tau }^{{}^{\prime }}}⟨j_i||{\Delta }_{L;\tau }^{{}^{\prime }}||j_f⟩⟨j_f||{\Sigma }_{L;{\tau }^{{}^{\prime }}}||j_i⟩],$$

$$⟨J_f;T_f{M}_{T_f}||A_{L,\tau }||J_i;T_i{M}_{T_i}⟩=(-1{)}^{T_f-M_{T_f}}\left(\begin{array}{ccc}{T}_f& \tau & T_i\\ M_{T_f}& M_{\tau }& M_{T_i}\end{array}\right)⟨J_f;T_f⋮⋮A_{L,\tau }⋮⋮J_i;T_i⟩.$$

$$⟨J_f;T_f⋮⋮A_{L,\tau }⋮⋮J_i;T_i⟩=\sum _{|\alpha ||\beta |}{\psi }_{L\tau }^{f,i}(|\alpha |,|\beta |)⟨|\alpha |⋮⋮A_{L,\tau }⋮⋮|\beta |⟩$$

$${\psi }_{L;\tau }(|\alpha |,|\beta |)=\sqrt{2(2J+1)(2T+1)(2j_{\alpha }+1)}{\delta }_{|\alpha ||\beta |}{\delta }_{L0}.$$

$$\frac{dR}{d{E}_R}=\frac{1}{m_T}\frac{{\rho }_{\chi }}{m_{\chi }}{\int }_{v\ge v_{min}}{d}^{3}vvf(\vec{v}+{\vec{v}}_E)\frac{d\sigma }{d{E}_R}$$

$$v_{min}=\sqrt{\frac{m_T{E}_R}{2{\mu }_T^2}}+\frac{E^{\ast }}{\sqrt{2m_T{E}_R}},$$

$$n_{\mathrm{q}}$$

$$N_{\mathrm{i}}$$

$$N_{\mathrm{ex}}$$

$$N_{\mathrm{reco}}=\mathrm{Gaus}(r{N}_{\mathrm{i}},{\sigma }_{\mathrm{r}}{N}_{\mathrm{i}}).$$

$$N_{\gamma }$$

$$N_{\mathrm{e}}$$

$$N_{\gamma }^{\det }=\mathrm{Binom}(N_{\gamma },g_1),$$

$$N_{\mathrm{e}}^{\det }=\mathrm{Binom}(N_{\mathrm{e}},p_{\mathrm{ex}}).$$

$$S1=\mathrm{Gaus}(N_{\gamma }^{\det },{\sigma }_{\mathrm{PMT}}\sqrt{N_{\gamma }^{det}}).$$

$$S2=\mathrm{Gaus}(\omega N_{\mathrm{e}},{\sigma }_{\omega }\sqrt{N_{\mathrm{e}}}).$$

$$f_{{\mathcal{O}}_{\alpha },m_{\chi }}$$

$$f_{{\text{R}}_{\beta }}$$

$${\mu }_i(R_{\chi },{\text{R}}_{\text{BG}})=\text{Exposure}\times \left(R_{\chi }\cdot f_{{\mathcal{O}}_{\alpha },m_{\chi }}(i)+\sum _{\beta }{R}_{\beta }\cdot f_{{\text{R}}_{\beta }}(i)\right).$$

$$\mathcal{L}(\mathcal{D}|R_{\chi },{\mathbf{R}}_{\text{BG}}):=\prod _{i=1}^{N_{\text{bins}}}\left[\frac{({\mu }_i(R_{\chi },{\text{R}}_{\text{BG}}){)}^{n_i(\mathcal{D})}}{n_i(\mathcal{D})!}\right]e^{-{\mu }_i(R_{\chi },{\text{R}}_{\text{BG}})}\cdot \prod _{\beta }{\mathcal{L}}_{\beta }(R_{\beta }),$$

$${\mathcal{L}}_{\beta }(R_{\beta }):=\frac{1}{\sqrt{2\pi {\sigma }_{\beta }^2}}\cdot \exp \left\{\frac{{\left(R_{\beta }-{\overline{R}}_{\beta }\right)}^2}{2{\sigma }_{\beta }^2}\right\}$$

$$\lambda (0):=\frac{\underset{{\text{R}}_{\text{BG}}}{\max }\mathcal{L}(\mathcal{D}|R_{\chi }=0,{\text{R}}_{\text{BG}})}{\underset{R_{\chi },{\text{R}}_{\text{BG}}}{\max }\mathcal{L}(\mathcal{D}|R_{\chi },{\text{R}}_{\text{BG}})}$$

$$q_0:=\{\begin{array}{ll}-2\ln \lambda (0)& R_{\chi }\ge 0\\ 0& R_{\chi }<0\end{array}.$$

$$q_A:=-2\ln \left[\frac{\underset{{\Theta }_A,{\text{R}}_{\text{BG}}}{\max }\mathcal{L}({\mathcal{D}}_A|{\Theta }_A,{\text{R}}_{\text{BG}})}{\underset{{\Theta }_B,{\text{R}}_{\text{BG}}}{\max }\mathcal{L}({\mathcal{D}}_A|{\Theta }_B,{\text{R}}_{\text{BG}})}\right]\text{and}{q}_B:=-2\ln \left[\frac{\underset{{\Theta }_A,{\text{R}}_{\text{BG}}}{\max }\mathcal{L}({\mathcal{D}}_B|{\Theta }_A,{\text{R}}_{\text{BG}})}{\underset{{\Theta }_B,{\text{R}}_{\text{BG}}}{\max }\mathcal{L}({\mathcal{D}}_B|{\Theta }_B,{\text{R}}_{\text{BG}})}\right],$$

$$p:={\int }_{q_{\text{z}}}^{\infty }f(q|{\mathcal{D}}_A)\text{d}q,$$

$$g=3{\kappa }^{0.15}+0.7{\kappa }^{0.6}+\kappa ,$$

Full text

1 Introduction

The search for dark matter (DM) in the form of weakly interacting massive particles (WIMPs) with direct detection experiments has made great progress in recent years [1, 2, 3, 4]. With the next generation of experiments on the way it is a timely question to ask what we can expect to learn about the properties of WIMPs once a signal has been established. We do know neither the mass of the DM particle nor other properties such as the spin. In addition, the interactions between WIMPs and Standard Model (SM) particles are not understood either. We hope to be able to shed more light on these parameters once a detection is established. In the following we will use the terms DM and WIMP interchangeably.

Direct detection aims to observe the recoil induced by DM scattering off nuclei. The typical velocity of DM bound to the Milky Way is limited by the galactic escape velocity $v_{esc}\approx 550\,\mbox{km}/\mbox{s}$ and, consequently, the highest nuclear recoil energy for an earth-based target assuming WIMPs with a mass of $100\,\mbox{GeV}$ is $\approx 300\;\mbox{keV}$. Therefore, direct detection experiments are testing the non-relativistic limit of the DM nucleus interaction and it is possible to describe the scattering process in terms of a non-relativistic effective theory of DM and nucleons. The general operator basis for non-relativistic DM-nucleon interactions and the formalism for the description of DM scattering off nuclei has been developed by [5]. The non-relativistic effective theory can be used for various applications [6, 7, 8, 9, 10, 11, 12, 13, 14]. Experimental constraints on the effective operators from dedicated analyses are reported in [15, 16].

Different combinations of WIMP mass and WIMP-nucleus interactions can lead to very similar signals in direct detection experiments. Therefore, a large exposure or additional signals are needed to gain information on the properties of DM. Typically, a large exposure is achieved with bigger and more sensitive detectors than available to date and is therefore very costly. However, direct detection experiments are sensitive to inelastic WIMP-nucleus scattering in addition to the standard signature of elastic WIMP nucleus scattering [17, 18, 19]. In this process, the nucleus is excited and due to the different nuclear structure of the ground and the excited state the response depends crucially on the DM-nucleon interaction. Exploiting the inelastic signal enhances the experimental sensitivity to DM properties without necessitating any modifications on the detector setup. This possibility is particularly attractive in the context of direct detection experiments employing dual-phase time projection chambers (TPC) with liquid xenon (LXe) targets [20, 21, 22]. In addition to being the leading technology for the standard channel of direct WIMP searches [2, 3, 4], these detectors are well suited for inelastic searches [23, 24, 25] since two common isotopes, 129Xe and 131Xe, have unusually low-energetic exited states with an excitation energy of $39.6\,\mbox{keV}$ and $80.2\,\mbox{keV}$, respectively. It is well known that inelastic scattering can help to distinguish the spin-dependent (SD) and spin-independent (SI) interactions traditionally considered in direct detection [18]. A detailed nuclear physics calculation of the inelastic SD DM nucleon scattering in xenon has been presented in [26] and the sensitivity of a XENONnT-like experiment with a $15\,\mbox{tonne}\times\mbox{year}$ exposure has been studied in [27]. The absence of signals of this type can be re-interpreted as an upper bound on conventional SD interactions with the strongest being of $\sim 10^{-38}\,{\mbox{GeV}}^{-2}$ for a WIMP mass of 200 GeV [25].

The aim of this work is three fold: First, we include in our study the full set of non-relativistic effective operators introduced in [5] rather then just the conventional spin-independent and spin-dependent interactions. Second, we analyze the discovery reach to these operators for future experiments in light of recent null results and derive the prospects for an ultimate future direct detection experiment such as DARWIN [28]. Third, we investigate the impact of inelastic scattering on our ability to pin down the DM-nucleon interaction once a signal has been detected.

The structure of this paper is as follows: We briefly introduce the non-relativistic effective DM-nucleon theory and present the modifications required to describe inelastic scattering in Sec. 2. Next, we discuss the general properties of the signal and outline the process of signal creation and detection in LXe TPCs in Sec. 3. The detector model employed for the following studies is described in detail. In Sec. 4 the prospects for discovering an inelastic signal with future experiments are discussed while Sec. 5 focuses on the use of the additional information gained from an analysis of inelastic scattering for the determination of the DM-nucleon interaction. Finally, we present our conclusions in Sec. 6. Additional information regarding the description of inelastic scattering and the detector model are presented in the Appendix.

2 Inelastic scattering in non-relativistic effective theory

In our analysis we use the parametrization of DM-nucleon interactions in terms of non-relativistic effective operators [5]111The operators entering at the nucleon-DM level can be mapped to the operators of a theory of relativistic DM [6, 29, 30]. Since we are interested in low-energy observables we prefer to remain agnostic regarding the UV-completion.. In this section we briefly review the formalism for calculating the DM nucleus scattering rates and comment on the extensions that are required to incorporate inelastic scattering. For a detailed introduction to the non-relativistic effective field theory of DM and nucleons we refer to [5, 31]. A more concise summary can be found for example in [9]. A Mathematica package which facilitates the computation of the different cross sections for elastic scattering has been published [31, 32]. We base our implementation of inelastic scattering on this package and extend it where necessary.

In the non-relativistic limit relevant for direct detection experiments, the form of the interaction between DM and nucleons is constrained by energy and momentum conservation and by Galilean invariance, i.e. invariance with respect to shifts in the three-dimensional particle velocity. Interactions respecting these constraints can be described by combinations of five operators, $1,$ $\vec{S}_{\chi}$, $\vec{S}_{N}$, $\vec{v}^{\perp}=\vec{v}+\frac{\vec{q}}{2\mu_{N}}$, and $i\vec{q}$, which act on the two-particle Hilbert space spanned by the nucleon and DM states. In the following, we limit ourselves to fermionic DM and to operators that are at most linear in $\vec{S}_{\chi}$, $\vec{S}_{N}$ and $\vec{v}^{\perp}$. With these constraints 14 linearly independent operators can be constructed222 Note that a further operator, $\mathcal{O}_{2}=\vec{v}^{\perp}\cdot\vec{v}^{\perp}$, cannot arise as the leading operator from a complete UV-theory and has therefore been omitted in [31]. For spin-1 DM four additional operators can be generated [11, 33]. We will not consider them here but they can be straightforwardly included in an analysis of inelastic scattering., see Tab. 1.

Linear combinations of these operators form the Hamiltonian density

[TABLE]

where the operators $t^{\tau}$ act on isospin states with $t^{0}=1$ and $t^{1}=\tau_{3}$ (third Pauli matrix). The coefficients $c_{k}^{t}$ parametrize the strength of the isoscalar and isovector interactions and can be recast in terms of the interaction strength with protons and neutrons via $c_{k}^{p}=(c^{0}_{k}+c^{1}_{k})/2$ and $c_{k}^{n}=(c^{0}_{k}-c^{1}_{k})/2$. Note that $c_{k}^{i}$ have dimension of mass*-2* due to the normalization convention chosen in [31]. Consequently, we are implicitly encoding in the $c_{k}^{p,n}$ coefficients a New Physics scale $\Lambda$ such that $c_{k}^{i}=\frac{1}{\Lambda^{2}}$.

The differential recoil cross section for DM scattering off a target nucleus of mass $m_{T}$ and spin $J$ is given by

[TABLE]

where $\langle|\mathcal{M}_{\rm eff}|^{2}\rangle$ denotes the spin-summed and spin-averaged square matrix element of DM-nucleus scattering and $v$ is the speed of the incoming DM particle in the rest frame of the target nucleus. Using standard techniques from nuclear physics the averaged matrix element can be expressed as

[TABLE]

where $R_{A}^{\tau\tau^{\prime}}$ are DM response functions which are known analytically and are given in Appx. B. The nuclear response functions $W^{\tau\tau^{\prime}}_{AB}=\sum_{L}\langle j_{i}||A_{L,\tau}||j_{f}\rangle\langle j_{f}||B_{L,\tau^{\prime}}||j_{j}\rangle$ are given by angular momentum reduced nuclear matrix elements of definite angular momentum $\langle j_{f}||A||j_{i}\rangle$. They depend on nuclear matrix elements and have to be evaluated numerically. The operators $A\in\{M,\Sigma^{\prime},\Sigma^{\prime\prime},\Phi^{\prime\prime},\tilde{\Phi}^{\prime},\Delta,\Delta^{\prime},\Sigma,\tilde{\Omega},\tilde{\Phi}$ ,$\tilde{\Delta}^{\prime\prime}\}$ arise from the multipole expansion of the nucleon operators in Tab. 1 and are well established from the treatment of semileptonic weak interactions with nuclei [34, 35, 36]333Following the notation of [5] the operator $\tilde{A}$ denotes the symmetrized version of operator $A$.. The first five lines of eq. 2 collect terms which are present in elastic and inelastic DM nucleus scattering and match the squared matrix element given in [31]. The terms in the last three lines are not considered in the standard non-relativistic theory of DM-nucleus interactions since they vanish for ground state to ground state transitions and are thus irrelevant for elastic scattering. However, they need to be included in an analysis of inelastic scattering.

Expanding the suppressed indices, the angular momentum reduced matrix elements $\langle j_{f}||A||j_{i}\rangle$ can be expressed in terms of the Wigner 3-j symbol and the isospin reduced matrix elements as

[TABLE]

Considering only one-body currents, the nuclear matrix elements can be further simplified. Expressed as the sum of products of single particle matrix elements with the elements of the one-body density matrix (OBDM) they read

[TABLE]

where $|\alpha|$ and $|\beta|$ represent the set of nonmagnetic single-nucleon spatial and spin quantum numbers and the sum extends over the complete set of $|\alpha|$ and $|\beta|$. In the harmonic oscillator basis the one-body nuclear matrix elements $\langle|\alpha|\vdots\vdots A_{J,T}\vdots\vdots|\beta|\rangle$ of the relevant operators are known analytically, see for example [9]. The one-body density matrix elements $\psi^{f,i}_{L;\tau}(|\alpha|,|\beta|)$ parametrize the overlap of the nucleon wave functions in the nucleus. They need to be determined in a nuclear physics computation. Following the literature we perform a nuclear shell model calculation. The contribution from closed shells in the core is given by  [37]

[TABLE]

Treating the outer shells is more difficult and their contribution can only be determined numerically. We use the public shell model code Nushell@MSU [38]. We perform calculations for 129Xe and 131Xe using an $3s_{1/2}2d_{3/2}2d_{5/2}1g_{7/2}1h_{11/2}$ model space above a 100Sn core, the “jj55pn” model space of Nushell@MSU, and the “sn100pn” interaction Hamiltonian [39]. For 131Xe we perform an unrestricted diagonalization of the valence nucleon system. For 129Xe the basis is too large and we were not able to perform an unrestricted diagonalization with our computational resources. Therefore, we need to truncate the model space. Following [26] we limit the number of neutron excitations into the energetically disfavored $1h_{11/2}$ orbital to three. A table with the values of OBDMEs in the format suitable for use with the DMFormFactor package is provided as supplementary material. In order to assess the validity of our shell-model computation we perform a comparison with the literature. We compare our results for the (in-)elastic nuclear structure functions for interactions with neutrons $S_{n}$ with [26]. As can be seen in Fig. 1 the overall agreement is good. Nevertheless, we would like to caution that our nuclear structure calculation is exploratory and would like to encourage the community to improve it. A more detailed study would be highly beneficial for a future search for inelastic scattering in direct detection experiments.

3 Inelastic scattering in LXe experiments

With the formalism for the description of nuclear excitation induced by DM scattering at hand we move forward by considering the experimental aspects in more detail. We will first analyze the signal rates associated with the different operators. After a brief introduction into the working principle of LXe TPCs we present the Monte Carlo (MC) simulation that converts the deposited energy of incident particles into measurable detector signals. Finally, we will introduce the relevant background sources and their rates and conclude by presenting the signal-background discrimination in LXe TPCs.

3.1 Signal rates

The scattering rates in a detector at earth are given by

[TABLE]

where $m_{T}$ and $m_{\chi}$ are the target and DM mass, respectively. The local DM density is denoted $\rho_{\chi}$ while $\vec{v}$ is the DM velocity. The velocity distribution of DM in the galactic rest frame is given by $f(\vec{v})$ and $\vec{v}_{E}$ is the velocity of Earth. We use the Standard Halo model which assumes a truncated Maxwell-Boltzmann distribution for $f(\vec{v})$. The input parameters are the solar circular speed, which we take to be $v_{0}=220\,\mbox{km}/\mbox{s}$, the speed of Earth $v_{E}=232\,\mbox{km}/\mbox{s}$ and the galactic escape speed $v_{esc}=550\,\mbox{km}/\mbox{s}$.

The lower boundary of the integration, $v_{min}$, is the minimal speed required to induce a nuclear recoil with energy $E_{R}$. It is given by

[TABLE]

where $\mu_{T}$ is the reduced mass of the DM-nucleus system and $E^{\ast}$ is the nuclear excitation energy.

As mentioned in the introduction, two xenon isotopes, 129Xe and 131Xe, have a low lying first excited state. Both isotopes are very common and have a natural abundance of $26.4\%$ and $21.2\%$ in xenon, respectively. The excitation energy of 129Xe is $39.6\,\mbox{keV}$ and the first excited state of 131Xe lies at $80.2\,\mbox{keV}$. While this $\mathcal{O}(10\,\mbox{keV})$ energy difference is small on nuclear scales it is sizable compared to the typical amount of energy available in DM-nucleus scattering and leads to a significant increase of $v_{min}$ compared to elastic scattering. Therefore, we should expect at least a moderate suppression of the inelastic rate with respect to the elastic rate.

It is instructive to take a look at the total rates of inelastic and elastic scattering. We report the ratio of the inelastic rate $R_{in}$ to the elastic rate $R_{el}$ in Tab. 2 for the representative choice of $m_{\chi}=300$ GeV and the same coupling to neutrons and protons (isoscalar interaction).

Six of the 14 basic operators predict an inelastic scattering rate which is less than $10^{-3}$ of the elastic scattering rate. Given that the elastic signal has not been detected yet, the observation of these inelastic signals is challenging and is unlikely to take place in near future. In contrast, the operators, $\mathcal{O}_{4}$, $\mathcal{O}_{5}$, $\mathcal{O}_{6}$, $\mathcal{O}_{9}$ and $\mathcal{O}_{10}$ feature an inelastic scattering rate that is only mildly suppressed by $\approx 5-15\%$. For those interactions, a future detection of inelastic transitions is conceivable. These very different outcomes for the various operators can be understood from the basic properties of elastic and inelastic scattering. As it is well known, the elastic scattering amplitude of the standard SI interaction, which corresponds to $\mathcal{O}_{1}$, is coherently enhanced and scales as $A^{2}$ where $A$ is the number of nucleons in the nucleus. Consequently, inelastic scattering is strongly suppressed for $\mathcal{O}_{1}$. In contrast, the standard SD interaction, which is described by $\mathcal{O}_{4}$, does not profit from the $A^{2}$ enhancement since the interaction leads to a spin-flip. By definition inelastic transitions change the structure of the nucleus and are therefore incoherent processes. Hence, the inelastic scattering cross section of $\mathcal{O}_{4}$ is not suppressed compared to the elastic one. This pattern can be generalized to the other operators and we find that basically every DM-nucleus interaction that is sensitive to the nuclear spin predicts interesting inelastic transition rates while the other operators do not. For most of these operators $R_{in}/R_{el}=\mathcal{O}(0.1)$ which agrees with the naive reasoning that the higher momentum transfer required by the inelastic process leads to a mild suppression. However, two operators, $\mathcal{O}_{7}$ and $\mathcal{O}_{13}$, predict that the inelastic rate exceeds the elastic one by one order of magnitude. This effect is more subtle and can ultimately be traced back to the nuclear structure; the elastic scattering rate of $\mathcal{O}_{7}$ and $\mathcal{O}_{13}$ is suppressed by $v_{T}^{\perp}$ while the inelastic rate receives additional contributions from the purely inelastic nuclear response functions $\tilde{\Omega}$ and $\tilde{\Phi}$ that do not suffer from this suppression. In the following, we will not consider $\mathcal{O}_{8}$, which is a borderline case with $R_{in}/R_{el}\approx 10^{-2}$, and focus on the operators with a clear detection potential.

Due to the even proton number in xenon, the dominating contribution to the interaction comes from unpaired neutrons in the two considered isotopes 129Xe and 131Xe. In addition, nuclear structure factors for protons suffer from substantial nuclear uncertainties which are known to impact the interpretation of experimental data [40, 1]. In the following, we assume a similar interaction strength of DM with neutrons and protons and will therefore only include the neutron coupling in our analysis.

3.2 Liquid xenon dual-phase detectors

In this section we assess the experimental observability of the introduced DM-nucleus interaction types.

The most sensitive direct detection DM detectors use the dual-phase time projection chamber (TPC) technology and feature a scalable LXe target that is contained in a cylindrical vessel. Radiation penetrating the target can deposit energy in the form of nuclear recoils (NRs) by scattering off the xenon nucleus (in case of WIMPs or neutrons) or in the form of electronic recoils (ERs) by interacting with atomic electrons (in case of $\gamma$-rays and $\beta$-electrons). The recoils excite and ionize neighbouring xenon atoms resulting in two distinct measurable quantities proportional to the recoil energy. For NRs some energy is lost to heat. The partition into measurable channels depends on the recoil energy and the recoil type. Consequently, the independent measurement of the excitation and ionization signals provides a discrimination between WIMP signals and ER backgrounds. Additionally, the two signals allow for a 3-dimensional reconstruction of the interaction’s position which is used to select only events from the inner radio-pure volume of the detector and therefore exploit the excellent self-shielding properties of LXe.

Xenon atoms excited by the recoil form dimers which de-excite on time scales $\mathcal{O}((1-10)\,\mathrm{ns})$ by emitting scintillation light with an average wavelength of $178\,$nm. The light signal (S1 signal) is typically observed by photomultiplier tubes (PMTs) that cover the top and bottom surfaces of the cylindrical volume. The ionization charges are extracted from the interaction site by an electric field ($\mathcal{O}(0.1\,\mathrm{kV/cm})$) applied across the target. The electrons are drifted upwards to the liquid-gas interface where they are extracted into the gaseous xenon phase by a second field ($\mathcal{O}(1\,\mathrm{kV/cm})$). The extraction field accelerates the electrons such that they excite and ionize the gas atoms inducing a second scintillation signal (S2 signal) that is proportional to the number of electrons [41].

The two lowest energy states of 129Xe and 131Xe feature lifetimes $\lesssim 1\,$ns and de-excite by emitting photons with energies of $39.6\,\mbox{keV}$ and $80.2\,\mbox{keV}$, respectively. In an inelastic scattering process of WIMPs, the NR and the ER induced by the de-excitation photon cannot be resolved by current detectors due to the short lifetime. Hence, both signals are measured simultaneously resulting in an event characteristic different from the standard elastic NR signal.

3.3 Detector Simulation

We use a MC simulation in order to convert an energy deposition in LXe into the detectable signals S1 and S2 and follow the same approach as the NEST model [42, 43, 44]. In the following, we use the symbol $n$ as an average number of quanta and $N$ as the number of quanta sampled from a probability density function. $N$ is therefore impacted by statistical fluctuations. Details on the parameter values used in the simulation can be found in Appx. A.

The average number of quanta $n_{\mathrm{q}}$ freed by a recoil of the energy $\epsilon$ divides into the number of excitons $n_{\mathrm{ex}}$ and the number of ions $n_{\mathrm{i}}$ and is given by

[TABLE]

with the average energy required per quantum $W=(13.7\pm 0.2)\,$eV and the Lindhard quenching factor $L$ that accounts for the energy lost to atomic motion. In order to consider fluctuations in the particle tracks, the number of quanta $N_{\mathrm{q}}$ is sampled from a Gaussian with the mean $n_{\mathrm{q}}$ and a standard deviation of $\sqrt{F\cdot n_{\mathrm{q}}}$, where $F$ is the Fano factor [45]. The partition into excitons and ions is determined by

[TABLE]

with the exciton-ion-ratio $n_{\mathrm{ex}}/n_{\mathrm{i}}$. All created excitons lead to the emission of a scintillation photon. However, ion-electron pairs can recombine, creating additional excitons and at the same time reducing the number of free electrons. The number of recombinations $N_{\mathrm{reco}}$ depends on the recombination probability $r$ and respective statistical fluctuations:

[TABLE]

The parameter $r$ is described by the Thomas-Imel box model [46] and depends on the applied electric field. In order to calculate the number of photons $N_{\gamma}$ emitted from the interaction site, another quenching factor $f_{\text{l}}$ has to be applied. This factor takes into account that two excitons can interact and produce only a single photon (Penning effects). Consequently, the final measurable quantities $N_{\gamma}$ and $N_{\mathrm{e}}$ are given by:

[TABLE]

After modelling the LXe microphysics, the quanta produced at the interaction site are propagated through the detector. The conversion of a scintillation photon into a measurable signal $N^{\mathrm{det}}_{\gamma}$ in a PMT is a binomial process

[TABLE]

with the detection probability $g_{1}$ that depends on detector specific parameters such as the geometry, reflectivity of materials, quantum efficiency of the PMTs, etc. On the side of the electron signal, the number of detected electrons $N_{\mathrm{e}}^{\mathrm{det}}$ depends on the probability $p_{\mathrm{ex}}$ of extracting the electrons from the liquid into the gaseous xenon phase:

[TABLE]

The measured S1 signal in number of photo-electrons (PE) emitted from the PMT photo-cathode is given after taking into account the resolution $\sigma_{\mathrm{PMT}}$ of the PMT response to photons:

[TABLE]

For S2 signals, the PMT response is not limiting the resolution due to the larger number of photons ($\mathcal{O}(10)$) generated during the proportional scintillation process. However, this process is also a statistical one and the electron amplification gain $\omega$ with its standard deviation $\sigma_{\omega}$ is determining the measured S2 signal:

[TABLE]

For NRs the procedure is followed exactly in the described way and the parameters are set to the values given in Tab. 4. For ERs the expression for the recombination probability is somewhat more complicated and depends on both, the electric field and the deposited energy $\epsilon$. Furthermore, incident gamma and beta particles show small differences in the partition of the deposited energy into the mean number of photons and electrons emitted and extracted from the interaction site, $n_{\gamma}$ and $n_{\text{e}}$, respectively. These differences are caused by the fact that beta particles only induce recoils on the outer shell electrons while gamma particles can also interact with inner shell electrons. The latter process results in additional Auger electrons. NEST implements the complete process and predicts $n_{\gamma}$ and $n_{\mathrm{e}}$ in dependence of $\epsilon$ and for various electric field strengths in [43] (Figure 1). For simplicity the provided data is utilized in this work and is smeared by the statistical fluctuations of the LXe microphysics processes. The total number of quanta produced at the interactions site is given by $n_{\mathrm{q}}=n_{\mathrm{e}}+n_{\gamma}$ which is valid for ERs due to the non-existing quenching ($L=1$). Furthermore, $n_{\gamma}$ and $n_{\mathrm{e}}$ enter into the calculation of $r$ (see Tab. 4). Given $n_{\mathrm{q}}$ and $r$ the full simulation including all statistical processes can be run for ERs as described above.

During inelastic scattering of WIMPs, a NR is induced together with an ER from the de-excitation photon. The signals from both recoils cannot be resolved in time. Hence, the S1 and S2 signals from the two processes are added at the end of this simulation resulting in S2/S1-ratios characteristic for inelastic scattering (see Fig. 3).

During the final phase of this work, the new version NESTv2 became available with updated charge and light yield curves for ERs [47]. The predicted new curves agree with the previous version within 20%. Due to time constraints we have not been able to implement the NEST update but we do not expect a large impact on our conclusions.

3.4 Detector configuration

Beside the reliance on the knowledge on LXe microphysics, the above described detector simulation depends on the detector configuration given by the electric field $E$, photon detection probability $g_{1}$, electron extraction probability $p_{\mathrm{extract}}$, electron amplification gain $\omega$ with its width $\sigma_{\omega}$, and the PMT single photon resolution $\sigma_{\mathrm{PMT}}$.

The most sensitive xenon detectors constructed to date are the XENON1T [2], LUX [3] and PandaX-II [4] experiments. An overview on the respective values for the mentioned parameters is given in Tab. 3. While the anticipated electric drift field across the TPC for XENON1T and LUX was $1$ and $2\,$kV/cm, respectively [48, 21], the achieved values were lower. However, higher drift fields have been used by the PandaX-II ($400\,$V/cm) and previous experiments such as XENON100 ($530\,$V/cm) [49] and ZEPLIN-III ($3400\,$V/cm) [50]. The NEST model suggests that the discrimination between NRs and ERs is improved for larger electric fields, especially at high energies [47]. Hence, future detectors are still aiming at high electric drift fields of several hundred V/cm [51, 48]. We choose for our benchmark detector model a drift field of $500\,$V/cm. All other parameters are fixed to the values achieved in XENON1T since this experiment has shown the highest sensitivity to WIMP-nucleon interactions among LXe direct detection experiments to date [2, 52].

It is common within the community to use only the S2 signal S2b measured in the bottom PMT array for the WIMP analysis, instead of the full S2 signal. This is due to the production of the S2 signal in the xenon gas phase located a few centimeters below the top PMT array. While the signal is only detected by a few photo sensors in the top array, it distributes more uniformly on the bottom array. Hence, S2b is less affected by localized effects from non-functional PMTs. We follow this approach in our simulation and multiply $\omega$ and its width $\sigma_{\omega}$ by a factor of 0.4 which is the fraction of light observed in the bottom PMT array in XENON1T [53].

We simulated the detector threshold by removing all simulated events with $\mathrm{S1}<3\,$PE and $\mathrm{S2}<200\,$PE which is compatible with the region of interest for DM searches with recent LXe experiments [2, 3, 4]. Even though the signal detection efficiency also highly depends on the performance of the data selection, we do not consider any cut acceptances in the simulation. The combined cut acceptance usually has no strong energy dependence and takes values of around $90\,$%. Therefore, it could be also treated as a universal scale factor on the signal rate which has small impact on the findings in this work.

3.5 Backgrounds

The background in large scale LXe direct detection experiments is composed of ERs and NRs induced by ambient radioactivity and solar neutrinos. The following background components are the most relevant and challenging ones for present and upcoming DM searches:

• •

Materials: $\gamma$ and $\beta$ radiation from the detector materials induce ERs in the target. Tonne scale direct detection experiments are able to exploit the good self-shielding properties of the LXe by selecting an inner radio-pure fiducial volume, reducing this background to a negligible level.

• •

85Kr: A small remnant of natural Krypton with the $\beta$-emitting isotope 85Kr remains in xenon that has been commercially produced by air liquefaction. However, its concentration can be further reduced to a sub-dominant level by cryogenic distillation [58].

• •

222Rn: The radio-active noble gas 222Rn is produced within the uranium decay chain, emanates from the detector materials and mixes within the LXe target due to its long half-life of 3.8 days. The daughter isotope 214Pb decays via the emission of a $\beta$-electron contributing to the ER background in the low-energy region. During detector construction, materials are screened and selected for their low intrinsic radioactivity and Radon emanation [59, 60], and 222Rn can be removed from LXe by means of cryogenic distillation [61, 62] or adsorption on charcoal [63].

• •

136Xe: The isotope 136Xe is a double $\beta$ emitter with a half-life of $2.17\cdot 10^{21}\,$y [64] and a $Q$-value of $2458\,$keV. Its abundance in natural xenon yields 8.9% and its depletion offers one possibility to reduce this background.

• •

Solar $\nu$s: Solar neutrinos can scatter elastically off electrons, inducing ERs and presenting an irreducible background which is homogeneously distributed over the target.

• •

Neutrons: Neutrons are produced by spontaneous fission and ($\alpha$,n)-reactions within the Uranium and Thoron decay chains. The particles are emitted from the detector materials and propagate through the target. In large scale detectors the vast majority of neutrons ($\sim 80$%) scatter twice within the volume [53]. Hence, this background can be reduced by single scatter event selection. Furthermore, future detectors will employ neutron vetos in order to further reduce this background component [65]. Muon-induced neutrons are reduced to a negligible level in large scale detectors by active Cherenkov muon vetos.

• •

CE$\nu$NS: Coherent scattering of $pp$ and 7Be solar neutrinos off nuclei induces NRs in the region of interest and represents an irreducible background which ultimately limits the sensitivity of future LXe DM detectors [66].

The energy spectra of the discussed background components and respective uncertainties are extracted from the XENON1T prediction [48]. An exception is the uncertainty of the 136Xe background rate where we neglect shape uncertainties of the spectrum potentially induced by the phase space factor that were considered in the reference. Instead, we assume that the dominating uncertainty is coming from the half-life measurement which yields 3%. The relative rate uncertainties for materials, 85Kr, 222Rn, solar neutrinos, neutrons and CE$\nu$NS yield 10%, 20%, 10%, 2%, 17% and 17%, respectively.

In this article we derive prospects for future LXe DM detectors like XENONnT, LZ [65] and DARWIN [28] which will have reduced background levels achieved by exploiting the reduction methods explained above. In particular, we consider for a XENONnT-like experiment a reduction of the 222Rn induced background by a factor of 100 and for the 85Kr component a factor of 10 w.r.t. the XENON1T prediction, as done in [48]. Additionally, the background from neutrons is scaled down by a factor of 20 assuming a good performance of the neutron veto and the double scatter rejection. Fig. 2 shows the simulated S1 spectrum for all background components. In the standard WIMP search region (up to  100 PE) solar neutrinos will be the dominating background while the leading contribution is coming from 136Xe at higher energies. For a DARWIN-like detector, radiogenic neutrons are assumed to be negligible while the 222Rn and 85Kr contributions are conservatively set to the XENONnT level.

Backgrounds from detector artifacts such as events from accidentally paired lone S1 and lone S2 signals from light and charge insensitive regions or radio-activity at the TPC surfaces that suffer from charge loss and are mis-reconstructed inside the FV are neglected in this study.

3.6 Signal and background discrimination

The advantage of LXe TPCs is their ability to measure both the light and charge signal of the interaction. This permits a three dimensional position reconstruction which is a powerful tool for background reduction by fiducialization as explained above. Furthermore, the ratio between S2 and S1 signals of NRs is smaller than for ERs which allows to differentiate between the two processes and therefore provides a further background discrimination.

Fig. 3 shows the background and signal distribution for the effective operators $\mathcal{O}_{4}$ (left) and $\mathcal{O}_{6}$ (right) in log10(S2b/S1) versus S1, denoted as analysis space. We assume a 100 tonne$\times$year exposure and an effective scale of $\Lambda=1900\,\mbox{GeV}$ (cross section of $\sigma_{SD}=1.5\times 10^{-41}\,\mbox{cm}^{2}$) and $\Lambda=230\,\mbox{GeV}$. The WIMP mass is set to $m_{\mathrm{\chi}}=300\,$GeV. This corresponds to a signal of $\simeq 100$ events in the detector. The density of the background population is indicated by the color scale. Regions where NR and ER events are expected to be located are indicated by 10%-50%-90% (dashed-solid-dashed) contours lines and form band-shaped structures. Elastic WIMP signal events (olive markers) populate the NR band that is separated from the ER background. The mono-energetic de-excitation signals of 129Xe and 131Xe (dashed blue and cyan ellipses) are shifted towards smaller S2b values compared to the ER band that is populated by $\beta$-radiation induced events. This is due to different photon and electron yields for $\beta$ and $\gamma$ interactions in LXe [44]. For WIMPs scattering inelastically off the two xenon isotopes, the de-excitation signal is added to the NR, resulting in signal regions that are stretched towards higher S1 signals (solid blue and cyan lines). The expected number of inelastic signal events is indicated by blue and cyan markers. The displacement of the signal regions from the ER band allows to suppress the bulk of the background population and therefore enhances the sensitivity. It is interesting to note that the irreducible background for inelastic DM-nucleus scattering, that is inelastic neutrino-nucleus scattering, is highly suppressed [67]. Therefore, a detector that allows for a clear separation between electronic backgrounds and the inelastic DM signal, which consists of a coincident photon and nuclear recoil, could in principle operate in background free mode. Hence, sensitivities beyond the conventional neutrino floor [66] could be achieved.

The recoil spectrum for the effective operator $\mathcal{O}_{4}$ falls off exponentially towards higher energies while the spectrum for $\mathcal{O}_{6}$ features a maximum at $\sim$130 keV. Hence, elastic events for $\mathcal{O}_{4}$ populate only the low-S1 region while the signal population extends more towards higher S1 values for $\mathcal{O}_{6}$ making this signal well distinguishable from $\mathcal{O}_{4}$ and background also in the purely elastic scattering channel.

4 Discovery reach for inelastic scattering

In this section we study the observability of the inelastic signal contributions for the effective operators $\mathcal{O}_{\alpha}$, $\alpha\in\{4,5,6,7,9,10,13\}$, and a given detector configuration. We use frequentist statistics to derive the DM scattering rate that allows to discriminate the signal from backgrounds with a significance of at least $3\sigma$ in $90\%$ of experiments.

4.1 Statistical treatment

We employ a binned Log-Likelihood procedure for the detection of a positive signal as described in [68] and [69]. Using the results of Sec. 3 we derived the binned probability distributions $f$ for the DM signal, which depend on $m_{\chi}$ and the operator $\mathcal{O}_{\alpha}$, and the backgrounds $\beta$

[TABLE]

with $\beta\in\{\textrm{Solar }\nu\textrm{s},\,^{85}\textrm{Kr},\,^{136}\textrm{Xe},\,^{222}\textrm{Rn},\,\textrm{Neutrons},\,\textrm{CE}\nu\textrm{NS}\}$. We simulate $10^{6}$ events to determine each probability distribution and employ $100\times 150$ bins in $S1-\log_{10}(S2/S1)$ space. The theoretically expected number of events $\mu_{i}$ in the $i$th bin, $i=1,...,N_{\textrm{bins}}$, is given by

[TABLE]

Here $f(i)$ is the probability of finding an event in the $i$th bin and $R_{\beta}$ and $R_{\chi}$ are the background and the DM rates, respectively. We denote the different background rates collectively as $\textbf{R}_{\textrm{BG}}:=(R_{\textrm{Solar }\nu\textrm{s}},R_{{}^{85}\textrm{Kr}},R_{{}^{136}\textrm{Xe}},R_{{}^{222}\textrm{Rn}},R_{\textrm{Neutrons}},R_{\textrm{CE}\nu\textrm{NS}})$ in the following.

For a given set of data $\mathcal{D}$ the binned likelihood function is given by

[TABLE]

where $n_{i}(\mathcal{D})$ is the observed number of events in bin $i$. The background likelihood function

[TABLE]

parametrizes our knowledge about the expected background rates $\overline{R}_{\beta}$ and allows for deviations $\sigma_{\beta}$ according to the background uncertainties (see Sec. 3.5). The likelihood function is maximized both under the null hypothesis of no DM signal and by leaving all rates unconstrained. The ratio of the maxima

[TABLE]

defines the frequentist test statistic for discovery of a positive signal

[TABLE]

According to Wilks’ theorem [70] the test statistic approaches a chi-square distribution for an infinite amount of experimental trials if the null hypothesis is true. Combined with the results found by Wald [71] this allows to reject the background only hypothesis with a significance of $Z=\sqrt{q_{0}}$ for data $\mathcal{D}$ including a DM signal.

We define the discovery reach for each DM mass as the value of the DM rate $\overline{R}_{\chi}$ for which $90\%$ of experiments find a $q_{0}$-value with a statistical significance of $Z\geq 3$. For fixed exposure and DM mass we conduct $2500$ pseudo-experiments and generate mock data which we use to derive the $q_{0}$ distribution. We vary the rate to determine the rate $\overline{R}_{\chi}$ for which $Z\geq 3$ for $90\%$ of the $2500$ $q_{0}$-values in the set $\{q_{0}\}$. By repeating this procedure for different DM masses we are able to map out the discovery reach. In order to ensure the comparability of the elastic and inelastic case we derive both discovery reaches ourselves and do not rely on results from literature.

4.2 Results

We summarize the main results of this section in Fig. 4, Fig. 5 and Fig. 6. We start by reviewing the discovery potential of inelastic transitions induced by operator $\mathcal{O}_{4}$, which corresponds to DM with the usual spin-dependent interaction, before turning to $\mathcal{O}_{5}$, $\mathcal{O}_{6}$, $\mathcal{O}_{9}$ and $\mathcal{O}_{10}$, all of which have a substantial but subdominant inelastic scattering rate. Finally, we discuss $\mathcal{O}_{7}$ and $\mathcal{O}_{13}$ which lead to a signal that is dominated by inelastic scattering.

The left panel of Fig. 4 displays the discovery reach for elastic (inelastic) scattering of DM on natural xenon assuming the standard SD operator $\mathcal{O}_{4}$ for future direct detection experiments, i.e. a XENONnT-like (DARWIN-like) detector with 4 tonne$\times$years ($200$ tonne$\times$years) exposure. The current limit on this interaction type from XENON1T  [72] is depicted in purple. Alternatively, we present the discovery reach in terms of the suppression scale of the effective theory $\Lambda$ in the right panel. As can be seen, the current bounds allow for the detection of a sizable inelastic signal with a DARWIN-like detector. However, if elastic DM-nucleus scattering is not observed by XENONnT the room for a detectable inelastic signal shrinks substantially.

For the same experimental benchmark scenarios we show in Fig. 5 the discovery reach for the operators with a dominant elastic signal, i.e. $\mathcal{O}_{5}$, $\mathcal{O}_{6}$, $\mathcal{O}_{9}$ and $\mathcal{O}_{10}$, as a limit on the scale 444The differential scattering cross section for these operators does not have the same properties as the one for SD-interactions and a recast of the bound as a limit on the zero momentum transfer cross section is highly misleading.. The curve for the elastic reach plateaus at $m_{\chi}\approx 200\,\mbox{GeV}$ and remains fairly stable up to $1$ TeV. The inelastic reach is stronger than the elastic one at high mass for all operators and, therefore, a clear detection of the inelastic signal in DARWIN is possible. For lower masses the scale of new physics that can be probed in inelastic scattering decreases rapidly and the competition with the elastic search channel becomes very challenging.

$\mathcal{O}_{7}$ and $\mathcal{O}_{13}$ lead to an inelastic scattering rate that exceeds the elastic one substantially. Depending on the DM mass, inelastic scattering could actually constitute the discovery channel for DM. Hence, we compare the discovery reach of both the elastic and the inelastic signal for a XENONnT-like detector and an exposure of $4$ tonne$\times$years in Fig. 6. A search for an inelastic signal can easily outstrip the elastic search for $m_{\chi}\geq 60\,(50)$ GeV for $\mathcal{O}_{7}$ ($\mathcal{O}_{13}$). An analysis searching for inelastic scattering in already collected data could therefore produce the world’s best limits on $\mathcal{O}_{7}$ and $\mathcal{O}_{13}$ or, more optimistically, lead to the first detection of DM in direct detection experiments.

5 Discriminating DM nucleus interactions

A DM signal in a direct detection experiment will be the starting point for numerous investigations into the properties of DM and its interactions with SM particles. Initially, the statistics of only a few events will not allow for a detailed analysis of the structure of DM interactions. However, more precise statements will become possible once more data are collected. Searching for inelastic DM-nucleus scattering is a potentially powerful way to exploit the data from ongoing experiments and learn more about DM. In this section we will investigate whether the additional information encoded in inelastic scattering events can make a difference in our quest for the nature of DM and facilitate the identification of the true DM interaction. For related analyses without inelastic scattering see for example [73, 74, 75, 76].

We aim to assess whether it is possible to reject a default hypothesis in favor of the effective operator responsible for the interaction. To be more specific, we focus on interactions with inelastic scattering, take $\mathcal{O}_{4}$, i.e. the standard SD interaction, as our default hypothesis and determine the expected rejection power of this assumption in light of a hypothetical detection.

5.1 Statistical treatment

Given a signal generated by an operator $\mathcal{O}_{B}$ and an experimental background we aim to determine the statistical level with which we can reject the hypothesis that the signal is due to a different operator $\mathcal{O}_{A}$ as a function of experimental exposure. The operators will either describe only elastic scattering processes or additionally include inelastic interactions. By comparing the rejection power of these two hypotheses we can determine whether the inelastic operator contributions enhance the discrimination prospects.

We define two frequentist test statistics for hypotheses testing:

[TABLE]

where $\Theta_{i}:=(R_{i},m_{i})$ with $i\in\{A,B\}$. The likelihood function $\mathscr{L}$ remains as in Eq. (4.4) with the exception that the DM mass is treated as a free parameter. Both test statistics $q_{A}$ and $q_{B}$ test the hypothesis $\mathcal{H}_{A}$, which we take as null or default hypothesis ($A=\mathcal{O}_{4}$ in this work), against the alternative hypothesis $\mathcal{H}_{B}$ for different mock data. Each hypothesis is characterized by the binned probability distributions $f_{\mathcal{O}_{A/B},m_{A/B}}$ which are associated with the respective operators $\mathcal{O}_{A/B}$ and include either elastic or elastic plus inelastic effects. The two mock data sets $\mathcal{D}_{A/B}$ are generated from $f_{\mathcal{O}_{A/B},m_{A/B}}$ with fixed parameter values. A repetition of the calculation provides us with sets $\{q_{A}\}$ and $\{q_{B}\}$ from which we derive the distributions $f(q|\mathcal{D}_{A})$ and $f(q|\mathcal{D}_{B})$. The discrimination between different operators $\mathcal{O}_{A}$ and $\mathcal{O}_{B}$ requires the test statistics to be derived using mock data from both hypotheses since Wilk’s and Wald’s theorems apply in neither of them555In contrast to the rejection of the background-only hypothesis, the parameter spaces for the hypotheses testing of two distinct operators are not nested which violates a requirement of Wilk’s theorem..

In detail our procedure is as follows: We start by simulating the data $\mathcal{D}_{B}(\overline{\mathcal{R}}_{B},\overline{m}_{B})$ from the background distributions using Poisson random sampling and $f_{\mathcal{O}_{B},m_{B}}$ with a fixed DM mass $\overline{m}_{B}=300\,\textrm{GeV}$ and a fixed DM rate $\overline{\mathcal{R}}_{B}=7\cdot 10^{-6}\,\mbox{kg}^{-1}\,\mbox{days}^{-1}$. Then we derive $q_{B}$ using the DM rate $\mathcal{R}_{B}$, the background rates $\textbf{R}_{\textrm{BG}}$ and the DM mass $m_{B}$ as fit parameters. Repeating this calculation $2500$ times enables us to derive the distribution $f(q|\mathcal{D}_{B})$ and the median best fit mass $\hat{m}_{B}$ for the alternative hypothesis. The best fit mass $\hat{m}_{B}$ in turn is used as mass $\overline{m}_{A}$ for the production of $10000$ mock data $\mathcal{D}_{A}(\overline{\mathcal{R}}_{A},\overline{m}_{A})$ with $\overline{\mathcal{R}}_{A}=\overline{\mathcal{R}}_{B}$ which again are produced by Poisson random sampling the distribution over the bins. From that we find $\{q_{A}\}$ and in consequence $f(q|\mathcal{D}_{A})$. This strategy is applied to our two cases in which we either consider only elastic scattering or include additionally the inelastic channels to both operators.

Our next goal is to derive the $p$-value as a measure of rejection. The hypothesis $\mathcal{H}_{B}$ can be rejected with probability $p$ against hypothesis $\mathcal{H}_{A}$ by

[TABLE]

where $q_{\textrm{z}}\in\{q_{B}\}$, such that $\int_{q_{\textrm{z}}}^{\infty}f({q|\mathcal{D}_{B}})\,\textrm{d}q=z$. We define the confidence level as $z=0.95$. In the left panel of Fig. 7 we illustrate the definition of the $p$-value for two example hypotheses. The right panel of the figure shows how the overlapping area underneath the two distributions $f(q|\mathcal{D}_{A})$ and $f(q|\mathcal{D}_{B})$ shrinks in consequence of an increasing exposure and thus leads to a decreasing $p$-value.

As a result we find for each pair of DM operators the discrimination $p$-value as a function of the detector exposure for elastic only and elastic plus inelastic interactions. Finally, we compare the exposure values that allow to reject the hypothesis $\mathcal{H}_{A}$ with a statistical level of $2\,\sigma$ for the elastic only and elastic plus inelastic cases and deduce whether the inelastic signal contribution improves the discrimination between the two operators.

5.2 Results

For the calculation of the discrimination power between two operators under the assumption of purely elastic or additional inelastic DM scattering we assume an elastic scattering rate of $\mathcal{R}_{\chi}=7.0\times 10^{-6}\,\mbox{kg}^{-1}\,\mbox{days}^{-1}$ leading to $\approx 10$ events in the region of interest. This would allow for a $\approx 3\sigma$ detection in an experiment with $4$ tonne$\times$years exposure assuming a DM mass of $m_{\chi}=300$ GeV as a representative example of heavy WIMP. We consider the six possibilities $\mathcal{H}_{B}\in\{\mathcal{H}_{\mathcal{O}_{5}},\mathcal{H}_{\mathcal{O}_{6}},\mathcal{H}_{\mathcal{O}_{7}},\mathcal{H}_{\mathcal{O}_{9}},\mathcal{H}_{\mathcal{O}_{10}},\mathcal{H}_{\mathcal{O}_{13}}\}$ while $\mathcal{H}_{A}=\mathcal{H}_{\mathcal{O}_{4}}$. As observed in the previous section for the discovery reach, a clear separation between $\mathcal{O}_{7,13}$ and the other operators emerges. With only $10$ events, elastic signals from different operators are basically indistinguishable and the discrimination $p$-values yield $\approx 0.95$ and $\approx 0.66$ for $\mathcal{O}_{7}$ and $\mathcal{O}_{13}$, respectively. However, once the inelastic signal is considered the $p$-value drops well below $10^{-5}$. This behavior is expected for $\mathcal{O}_{7}$ and $\mathcal{O}_{13}$ from the enhanced scattering rate in the inelastic channel w.r.t. the elastic channel.

For the other four operators the inelastic signal rate is suppressed compared to the elastic one. $10$ elastic events correspond to $\approx 1$ inelastic scattering events. Therefore, the inelastic and the elastic signal lead to the same discrimination $p$-values of $\mathcal{O}(0.1)$ for this chosen exposure. In Fig. 8 we show the evolution of the $p$-values calculated for purely elastic DM scattering or an additional inelastic interaction as a function of exposure. The inelastic signal substantially improves the $p$-value for $\mathcal{O}_{5}$ (top left) and the exposure required to exclude the default hypothesis at the $p=0.05$ level, $E_{p=0.05}$, decreases by more than a factor of 2. $E_{p=0.05}$ is in reach of an extended run of the discovery device while reaching the same level of confidence only based on the elastic signal would probably require to scale up the detector’s target mass.

In the case of $\mathcal{O}_{9}$ the consideration of inelastic signals decreases $E_{p=0.05}$ from $\approx 11.5\;\mbox{tonne}\times\mbox{years}$ to $9\;\mbox{tonne}\times\mbox{years}$. For $\mathcal{O}_{6}$ and $\mathcal{O}_{10}$, the discrimination $p$-values in the vicinity of $E_{p=0.05}$ for the elastic-only case are $\approx 50\%$ weaker than those from the elastic plus inelastic channel.

Summarizing our findings, inelastic scattering events improve the discrimination power for four of the six considered operators and add a moderate improvement for the other two. Therefore, a search for inelastic scattering is going to provide valuable additional information. Given that the required data are acquired simultaneously to the elastic scattering data, the inelastic channel offers a potential to improve our understanding of DM without constructing additional experiments.

6 Conclusion

Direct detection experiments offer one of the most promising ways to search for DM. In light of the substantial gains in sensitivity of future experiments, a detection of the DM particle is an exciting possibility. Once the detection has been established, the next question will be what about the DM particle’s nature that can be extracted from the observations. Subleading signatures, which will become accessible once a detection has been established, might play a major role in answering this question. Direct detection experiments are not only sensitive to the standard signature, i.e. the energy deposit for a nucleus recoiling against the DM in an elastic scattering process, but can also probe inelastic scattering in which the nucleus is excited. Xenon-based direct detection experiments are particularly interesting in this context since they have shown leading sensitivities to DM and the two common isotopes 129Xe and 131Xe have low-lying first excited states. In this paper we quantified the physics reach of future xenon-detectors w.r.t. the full set of non-relativistic operators. We extend the non-relativistic effective theory to account for inelastic processes and perform a nuclear shell model calculation for the ground and the first excited state. Combined with a full Monte-Carlo simulation of a realistic two-phase xenon-detector this allows us to derive the discovery reach for inelastic scattering.

We find encouraging prospects for the detection of inelastic scattering with future xenon-based detectors for eight of the $14$ non-relativistic effective operators. In the case of two operators, $\mathcal{O}_{7,13}$, the prospects for discovery of DM induced nuclear excitations are even better than for the discovery of the conventional elastic signal. Consequently, an analysis of current data with inelastic scattering in mind would lead to the world’s best limits on these operators or might even reveal the first direct detection of DM. The additional information encoded in the inelastic scattering events can be used to learn more about the theory of DM once a signal is established. We have analyzed the power with which the standard hypothesis of DM interactions with inelastic interactions, i.e. the standard spin-dependent interaction $\mathcal{O}_{4}$, can be rejected if the data follows one of the other operators that lead to inelastic scattering. We found that including the inelastic signal improves the $p$-value significantly in four of the six cases under consideration. In the other cases the elastic signal itself is already distinct enough to shift the $p$-value to $0.05$ rather fast with increasing exposure and the benefit of using inelastic scattering events is less prominent.

In light of the excellent prospects for $\mathcal{O}_{7,13}$ we would like to encourage the experimental community to revisit the collected data with inelastic scattering in mind and perform a dedicated analysis. We also want to stress the importance of using all available data once a DM signal is discovered.

Acknowledgements

We thank M. Lindner, T. Marrodan Undagoitia and C. McCabe for helpful discussions about xenon-based dark matter detectors in general and inelastic scattering in particular. We also thank B. A. Brown for providing the NushellX@MSU code for shell model calculations.

Appendix A Detector simulation

In this section we give further details on the parameters employed in the detector simulation that converts energy deposition in LXe into the measured signals S1 and S2b. A list of the parameters and how they are treated for ERs and NRs is given in Tab. 4:

• •

The Lindhard quenching factor $L$ quantifies the fraction of the energy that is converted into visible signals, i.e. light or ionization charges. For ERs this fraction is assumed to be 100% while part of the energy is lost to atomic motion for NRs. The parameter $k$ in the formula give in Tab. 4 relates the electronic stopping power and the velocity of NRs and yields 0.1394. The function $g(\kappa)$ is parametrized by

[TABLE]

with $\kappa=11.5\,E\cdot Z^{-7/3}$ that depends on the electric field $E$ and the proton number $Z$, i.e. 54 for xenon.

• •

The Fano factor $F$ is a measure for the fluctuation of the number of ions and scintillation photons produced by an energy deposition and has been measured to yield a value of 0.03 for ionizing radiation. Since there are no dedicated measurements for $F$ in LXe for NRs, we assume the same value for NRs as for ERs.

• •

The ratio $n_{\mathrm{ex}}/n_{\mathrm{i}}$ of the number of excitons to ions created in the interaction is fixed for ERs while it follows a field dependent function for NRs (see Tab. 4).

• •

The probability $r$ for an electron-ion pair to recombine after its creating is determined by the Thomas-Imel box model [46]. For NRs it takes the parameters $\xi=\gamma\cdot E^{-\delta}$, $\gamma=0.014$ and $\delta=0.062$ which have been estimated from experimental data. The expression for $\xi$ in the case of ERs is more complicated and depends on both, the deposited energy $\epsilon$ and $E$. For simplicity, $n_{\gamma}$ and $n_{\mathrm{e}}$, i.e. the mean number of photons and electrons emitted and extracted from the interaction site, are obtained from the average light and charge yield curves for ERs given in [43]. The numbers are used to re-calculated the recombination fraction. During this work the new version NESTv2 became available [47]. The predicted light and charge yield values agree with the previous version within 20%.

• •

The standard deviation of recombination fluctuations $\sigma_{\mathrm{r}}$ yields a constant value of 0.06 for ERs while for NRs it has been shown to be dependent on $r$.

• •

In the case of NRs two excitons can interact to produce a single photon leading to an additional quenching factor $f_{\text{l}}$ (Penning quenching) that is parametrized by $\kappa$ (see Lindhard quenching), $\eta=3.32$ and $\lambda=1.14$.

Appendix B DM response function

[TABLE]

[TABLE]