capbtabboxtable[][]

# Prospects for Dark Matter signal discovery and model selection via timing
information in a low-threshold experiment

Riccardo Catena  and Vanessa Zema

###### Abstract

In the recent years, many low-threshold dark matter (DM) direct detection
experiments have reported the observation of unexplained excesses of events at
low energies. Exemplary for these, the experiment CRESST has detected
unidentified events below an energy of about 200 eV – a result hampering the
detector performance in the search for GeV-scale DM. In this work, we test the
impact of nuclear recoil timing information on the potential for DM signal
discovery and model selection on a low-threshold experiment limited by the
presence of an unidentified background resembling this population of low-
energy events. Among the different targets explored by the CRESST
collaboration, here we focus on Al2O3, as a sapphire detector was shown to
reach an energy threshold as low as 19.7 eV [1]. We test the ability of a low-
threshold experiment to discover a signal above a given background, or to
reject the spin-independent interaction in favour of a magnetic dipole
coupling in terms of $p$-values. We perform our $p$-value calculations: 1)
taking timing information into account; and 2) assuming that the latter is not
available. By comparing the two approaches, we find that under our assumptions
timing information has a marginal impact on the potential for DM signal
discovery, while provides more significant results for the selection between
the two models considered. For the model parameters explored here, we find
that the $p$-value for rejecting spin-independent interactions in favour of a
magnetic dipole coupling is about 0.11 when the experimental exposure is 460
g$\times$year and smaller (about 0.06) if timing information is available. The
conclusion on the role of timing information remains qualitatively unchanged
for exposures as large as 1 kg$\times$5 year. At the same time, our results
show that a 90% C.L. rejection of spin-independent interactions in favour of a
magnetic dipole coupling is within reach of an upgrade of the CRESST
experiment [2].

##  1 Introduction

While the identity of our Universe’s invisible mass component, Dark Matter
(DM), remains unknown, the experimental search for its microscopic
constituents progressed rapidly in recent years [3]. In particular, the search
for nuclear recoils induced by the non-relativistic scattering of Milky Way DM
particles in low-background deep underground detectors – the so-called DM
direct detection technique [4, 5] – played, and will continue playing, a
central role in this context [6]. For example, the null result of the
operating experiments led to stringent constraints on the DM-nucleon and
-electron coupling, e.g. [7], and the effort placed on assessing the expected
performance of next-generation direct detection experiments has been
remarkable in recent years [8, 9, 10, 11, 12, 13, 14, 15]. In this assessment,
the performance of an experiment is measured by the exposure it must reach to
reject the background-only hypothesis for a specific background model, or to
reject a particle physics model in favour of an alternative one at a given
significance. This assessment is crucial, as it contributes determining which
experiments to prioritise in the coming years. One important component of this
assessment is to clarify quantitatively whether the simultaneous measurement
of nuclear recoil energy spectrum and associated time distribution of the
observed recoils is advantageous for a given experimental setup, or it is
preferable to focus on the measurement of the energy spectrum alone. This
question must be addressed proactively, as building ultra-low background
experiments that are stable over the time scale of one year (the period of the
Earth revolution around the Sun) is a difficult experimental approach posing
highly non trivial challenges in particular for cryogenic detectors.

The impact of timing information on the ability of direct detection
experiments to distinguish DM models with nearly degenerate recoil spectra has
been investigated for single-target experiments employing fluorine, germanium
or xenon targets and with a relatively high energy threshold in [16]. Using a
Bayesian approach to model selection, Ref. [16] finds that including timing
information may enhance this ability, but only with exposures beyond the
expectation of the experiments currently running (indicated as Generation 2
experiment in [16]). The results of Ref. [16] apply to weakly interacting
massive particles of mass around 50 GeV and cannot a priori be extended to DM
at the GeV scale. This extension would require lowering the assumed energy
threshold below about 0.3 keV.

Similarly, the impact of timing information on the ability of direct detection
experiments to identify a DM signal over an experimental background has not
been investigated for low-threshold experiments. Consequently, it remains
unclear whether it is worth taking the challenge of operating cryogenic
experiments under stable conditions over a period of one year or longer.
Importantly, this is a timely question, as the low-threshold experiment CRESST
has recently reported unexplained low-energy events in one of the detector
modules [17]. These start at the detector energy threshold, 30.1 eV, and
extend up to 200 eV, inside the acceptance region in the light yield - recoil
energy plane. The CRESST collaboration reports that the excess is present in
all the operated modules and that its energy distribution varies from module
to module, observation which disfavours a single common origin of this effect.
However, this excess limits CRESST sensitivity. Consequently, it is reasonable
to ask whether a time-dependent analysis improves the significance for signal
identification and model selection in the presence of such background.

In this work, we investigate the impact of timing information on the potential
for DM signal discovery and model selection of a low-threshold experiment. The
analysis is divided into two parts. In the first one, we test the impact of
timing information on the ability of a low-threshold experiment to discover a
DM signal over an experimental background that resembles the population of
low-energy events found in [17]. The second part of the analysis focuses on DM
model selection. Here, we test the ability of a low-threshold experiment to
reject a model where DM couples to nucleons via spin-independent interactions
in favour of an alternative hypothesis where DM couples to nuclei via magnetic
dipole interactions. Our choice of alternative hypothesis is motivated by the
fact that the time of maximum DM-induced nuclear recoil rate as a function of
the minimum velocity required to transfer a given momentum in the scattering
can be target-dependent in the case of magnetic dipole interactions, whereas
it does not depend on the target for spin-independent interactions [18].
Consequently, a multi-target experiment like CRESST should be able to exploit
this feature when comparing different DM models. Remarkably, we find that
timing information has not a significant impact on DM signal discovery, but on
model selection for experimental exposures that are within reach of next-
generation, low-threshold direct detection experiments. For example, we find
that with an exposure of 460 g$\times$year the $p$-value for rejecting spin-
independent interactions in favour of a magnetic dipole coupling is around
0.11, while it is about 0.06 when the time distribution of the event rate is
available. At the same time, our results show that a 90% C.L. rejection of
spin-independent interactions is within reach of a future upgrade of the
CRESST experiment [2].

This work is organised as follows. In Sec. 2 we outline the theoretical
framework used in our analysis. In Sec. 3, we critically review the timing
information that is available to a low-threshold, multi-target experiment like
CRESST. In Sec. 4, we introduce the statistical framework that we apply to
assess the prospects for DM signal discover and model selection with a low-
threshold experiment. We present our results on DM signal discovery in Sec.
5.1 and on DM model selection in Sec. 5.2. We finally summarise and conclude
in Sec. 6. In an Appendix, we list useful expressions to compute scattering
rates and cross sections in the case of spin 1/2 DM.

##  2 Theoretical framework

The differential rate of DM-nucleus scattering events per unit detector mass
in a direct detection experiment can be expressed as an integral over DM
particle velocities in the detector rest frame,

$\displaystyle\frac{{\rm d}R}{{\rm
d}E_{R}}=\sum_{T}\xi_{T}\frac{\rho_{\chi}}{m_{\chi}m_{T}}\int_{|\boldsymbol{v}|\geq
v_{\rm min}}{\rm d}^{3}v\,|\boldsymbol{v}|f(\boldsymbol{v},t)\,\frac{{\rm
d}\sigma_{T}}{{\rm d}E_{R}}\,,$ (2.1)

where $v_{\rm min}=\sqrt{2m_{T}E_{R}}/(2\mu_{T})$ is the minimum velocity a DM
particle must have to deposit an energy $E_{R}$ in the detector, $m_{T}$ and
$\mu_{T}$ are the target nucleus and DM-nucleus reduced mass, respectively,
$m_{\chi}$ is the DM particle mass and $\rho_{\chi}$ is the local DM density.
In multi-target detectors, like CRESST, the contribution to Eq. (2.1) of each
element in the detector material is weighted by the corresponding mass
fraction, $\xi_{T}$. The differential rate of DM-nucleus scattering events
also depends on the differential cross section for DM-target nucleus
scattering, ${\rm d}\sigma_{T}/{\rm d}E_{R}$, and on the local DM velocity
distribution expressed as a function of the DM velocity in the detector
reference frame, $f(\boldsymbol{v},t)$. The velocity distribution
$f(\boldsymbol{v},t)$ is a periodic function of time, $t$, with period of one
year. As a result, the expected rate of nuclear recoils is a function of time
with the same periodicity. We model $f(\boldsymbol{v},t)$ by taking the
gravitational focusing of the Sun into account, which implies [19],

$f(\boldsymbol{v},t)=\tilde{f}(\boldsymbol{v_{\odot}}+\boldsymbol{v_{\infty}}[\boldsymbol{v}+\boldsymbol{V}_{\oplus}(t)])\,,$
(2.2)

where $\tilde{f}$ is the DM velocity distribution in the Galactic reference
frame (where, by construction, the mean velocity of DM particles is zero),
$\boldsymbol{v_{\odot}}\simeq(11,232,7)$ km s-1 is the relative velocity of
the Sun with respect to the Galactic centre in Galactic coordinates,
$\boldsymbol{v}$ is the DM particle velocity in the detector rest frame and
$\boldsymbol{V}_{\oplus}(t)$ is the relative velocity of the Earth with
respect to the Sun. Consistently, we denote by
$\boldsymbol{v_{s}}=\boldsymbol{v}+\boldsymbol{V}_{\oplus}$ the DM particle
velocity in the Sun rest frame. With this notation, the solar frame velocity
$\boldsymbol{v_{\infty}}$ a DM particle must have at infinity to move with
velocity $\boldsymbol{v_{s}}$ with respect to the Sun at the detector position
is given by [19],

$\boldsymbol{v_{\infty}}[\boldsymbol{v_{s}}]=\frac{v_{\infty}^{2}\boldsymbol{v_{s}}+v_{\infty}(GM_{\odot}/r_{s})\boldsymbol{\hat{r}_{s}}-v_{\infty}\boldsymbol{v_{s}}(\boldsymbol{v_{s}}\cdot\boldsymbol{\hat{r}_{s}})}{v_{\infty}^{2}+(GM_{\odot}/r_{s})-v_{\infty}(\boldsymbol{v_{s}}\cdot\boldsymbol{\hat{r}_{s}})}\,,$
(2.3)

where $r_{s}$ is Earth-Sun distance, $\boldsymbol{\hat{r}_{s}}$ is the time
dependent unit vector that points from the Sun to the Earth,
$v_{\infty}^{2}=v^{2}-2GM_{\odot}/r_{s}$, $M_{\odot}$ is the mass of the Sun,
and $G$ is Newton’s constant. Explicit expressions for
$\boldsymbol{V_{\oplus}}$ and $\boldsymbol{\hat{r}_{s}}$ in Galactic
coordinates can be found in [20] and [21, 22], respectively. For $\tilde{f}$,
we assume a Gaussian distribution truncated at the escape velocity $v_{\rm
esc},$

$\displaystyle\tilde{f}(\boldsymbol{\tilde{v}})=\frac{1}{N_{\rm
esc}(2\pi\sigma_{v}^{2})^{3/2}}\,\exp\left(-\frac{|\boldsymbol{\tilde{v}}|^{2}}{2\sigma_{v}^{2}}\right)\Theta\left(v_{\rm
esc}-|\boldsymbol{\tilde{v}}|\right)\,,$ (2.4)

where $\boldsymbol{\tilde{v}}=\boldsymbol{v_{\odot}}+\boldsymbol{v_{\infty}}$.
Here, we assume the one-dimensional velocity dispersion $\sigma_{v}=164$ km
s-1 and a local escape velocity of $v_{\rm esc},=544$ km s-1. The
normalisation factor, $N_{\rm esc}$, in Eq. (2.4) is given by,

$N_{\rm esc}=\mathop{\mathrm{erf}}\left(\frac{v_{\rm
esc}}{\sqrt{2}\sigma_{v}}\right)-\sqrt{\frac{2}{\pi}}\frac{v_{\rm
esc}}{\sigma_{v}}\exp\left(-\frac{v_{\rm esc}^{2}}{2\sigma_{v}^{2}}\right)\,.$
(2.5)

$\mathcal{O}_{1}=\mathds{1}_{\chi N}$ | $\mathcal{O}_{9}=i\mathbf{S}_{\chi}\cdot\left(\mathbf{S}_{N}\times\frac{\mathbf{q}}{m_{N}}\right)$
---|---
$\mathcal{O}_{3}=i\mathbf{S}_{N}\cdot\left(\frac{\mathbf{q}}{m_{N}}\times\mathbf{v}^{\perp}\right)$ | $\mathcal{O}_{10}=i\mathbf{S}_{N}\cdot\frac{\mathbf{q}}{m_{N}}$
$\mathcal{O}_{4}=\mathbf{S}_{\chi}\cdot\mathbf{S}_{N}$ | $\mathcal{O}_{11}=i\mathbf{S}_{\chi}\cdot\frac{\mathbf{q}}{m_{N}}$
$\mathcal{O}_{5}=i\mathbf{S}_{\chi}\cdot\left(\frac{\mathbf{q}}{m_{N}}\times\mathbf{v}^{\perp}\right)$ | $\mathcal{O}_{12}=\mathbf{S}_{\chi}\cdot\left(\mathbf{S}_{N}\times\mathbf{v}^{\perp}\right)$
$\mathcal{O}_{6}=\left(\mathbf{S}_{\chi}\cdot\frac{\mathbf{q}}{m_{N}}\right)\left(\mathbf{S}_{N}\cdot\frac{\hat{{\bf{q}}}}{m_{N}}\right)$ | $\mathcal{O}_{13}=i\left(\mathbf{S}_{\chi}\cdot\mathbf{v}^{\perp}\right)\left(\mathbf{S}_{N}\cdot\frac{\mathbf{q}}{m_{N}}\right)$
$\mathcal{O}_{7}=\mathbf{S}_{N}\cdot\mathbf{v}^{\perp}$ | $\mathcal{O}_{14}=i\left(\mathbf{S}_{\chi}\cdot\frac{\mathbf{q}}{m_{N}}\right)\left(\mathbf{S}_{N}\cdot\mathbf{v}^{\perp}\right)$
$\mathcal{O}_{8}=\mathbf{S}_{\chi}\cdot\mathbf{v}^{\perp}$ | $\mathcal{O}_{15}=-\left(\mathbf{S}_{\chi}\cdot\frac{\mathbf{q}}{m_{N}}\right)\left[\left(\mathbf{S}_{N}\times\mathbf{v}^{\perp}\right)\cdot\frac{\mathbf{q}}{m_{N}}\right]$
Table 1: Interaction operators defining the non-relativistic effective theory
of spin 1/2 DM-nucleon interactions [23, 24]. $\mathbf{S}_{N}$
($\mathbf{S}_{\chi}$) is the nucleon (DM) spin,
$\mathbf{v}^{\perp}=\mathbf{v}-\mathbf{q}/(2\mu_{N})$, where $\mu_{N}$ is the
DM-nucleon reduced mass, is the transverse relative velocity and
$\mathds{1}_{\chi N}$ is the identity in the DM-nucleon spin space.

The differential cross section for DM-nucleus scattering, ${\rm
d}\sigma_{T}/{\rm d}E_{R}$ encodes all particle and nuclear physics inputs
required to evaluate Eq. (2.1). Within the non relativistic effective theory
of DM-nucleon interactions, ${\rm d}\sigma_{T}/{\rm d}E_{R}$ can be expressed
as follows,

$\displaystyle\frac{{\rm d}\sigma_{T}}{{\rm
d}E_{R}}=\frac{2m_{T}}{(2J_{T}+1)v^{2}}\sum_{k}\sum_{\tau,\tau^{\prime}}\left(\frac{q^{2}}{m_{N}^{2}}\right)^{\ell_{k}}R^{\tau\tau^{\prime}}_{k}\left(v_{T}^{\perp
2},{q^{2}\over m_{N}^{2}}\right)W_{k}^{\tau\tau^{\prime}}(q^{2})\,,$ (2.6)

where $J_{T}$ is the target nucleus spin, $v=|\mathbf{v}|$, $v_{T}^{\perp
2}=v^{2}-q^{2}/(4\mu_{T}^{2})$, $q=\sqrt{2m_{T}E_{R}}$ is the momentum
transfer and $m_{N}$ is the nucleon mass. At most eight DM and nuclear
response functions, $R^{\tau\tau^{\prime}}_{k}$ and
$W^{\tau\tau^{\prime}}_{k}$,
$k=M,\Sigma^{\prime},\Sigma^{\prime\prime},\Phi^{\prime\prime},\Phi^{\prime\prime}M,\tilde{\Phi}^{\prime},\Delta,\Delta\Sigma^{\prime}$,
respectively, can appear in Eq. (2.6) [24]. The $R^{\tau\tau^{\prime}}_{k}$
functions depend on $q^{2}/m_{N}^{2}$, $v_{T}^{\perp 2}$ and on the coupling
constants for DM-nucleon interactions, $c_{j}^{\tau}$. They are known
analytically, and we list them for spin 1/2 DM in the Appendix. Here, the
index $j$ labels the type of DM-nucleon interaction [24]. Assuming one-body
DM-nucleon interactions, there are $4+20J_{\chi}$ interaction types that are
invariant under Galilean transformations and spatial rotations for a DM
particle of spin $J_{\chi}$ [25]. We list the 14 independent interaction
operators arising for spin 1/2 DM in Tab. 1. The eight nuclear response
functions $W_{k}^{\tau\tau^{\prime}}$ in Eq. (2.6) are quadratic in reduced
matrix elements of nuclear charges and currents, and must be computed
numerically. This calculation is performed within the nuclear shell model for
heavy targets, such as xenon. However, for light nuclei ab initio methods have
recently been applied [26]. The indexes $\tau$ and $\tau^{\prime}$ run from 0
to 1: 0 corresponds to “isoscalar” interactions and 1 to “isovector”
interactions [24]. In terms of, e.g. $c_{1}^{\tau}$, the neutron to proton
coupling ratio can be expressed as
$r=c_{p}/c_{n}=(c_{1}^{0}-c_{1}^{1})/(c_{1}^{0}+c_{1}^{1})$. Finally,
$\ell_{k}=0$ for $k=M,\Sigma^{\prime},\Sigma^{\prime\prime}$, and $\ell_{k}=1$
otherwise.

The way the cross section ${\rm d}\sigma_{T}/{\rm d}E_{R}$ depends on $v$
determines how the nuclear recoil rate depends on $t$. Inspection of Eq.
(LABEL:eq:R) shows that there are four ways the velocity $v$ can enter the
cross section ${\rm d}\sigma_{T}/{\rm d}E_{R}$. The interactions labelled by
$\mathcal{O}_{1}$, $\mathcal{O}_{7}$, $\mathcal{O}_{8}$ and $\mathcal{O}_{11}$
in the literature [24] and here associated with the coupling constants
$c_{1}^{\tau}$, $c_{7}^{\tau}$, $c_{8}^{\tau}$ and $c_{11}^{\tau}$ are
representative of these four classes. The corresponding DM-nucleus scattering
cross section reads as,

$\displaystyle\frac{{\rm d}\sigma_{T}}{{\rm
d}E_{R}}\Bigg{|}_{\mathcal{O}_{1}}$
$\displaystyle=\frac{2m_{T}}{(2J_{T}+1)v^{2}}\sum_{\tau\tau^{\prime}}c_{1}^{\tau}c_{1}^{\tau^{\prime}}W_{M}^{\tau\tau^{\prime}}(q^{2})\equiv\frac{a(q^{2})}{v^{2}}\,,$
$\displaystyle\frac{{\rm d}\sigma_{T}}{{\rm
d}E_{R}}\Bigg{|}_{\mathcal{O}_{7}}$
$\displaystyle=\frac{2m_{T}}{(2J_{T}+1)v^{2}}\sum_{\tau\tau^{\prime}}c_{7}^{\tau}c_{7}^{\tau^{\prime}}\frac{(v^{2}-v_{\rm
min}^{2})}{8}W_{M}^{\tau\tau^{\prime}}(q^{2})\equiv
b(q^{2})\left(1-\frac{v_{\rm min}^{2}}{v^{2}}\right)\,,$
$\displaystyle\frac{{\rm d}\sigma_{T}}{{\rm
d}E_{R}}\Bigg{|}_{\mathcal{O}_{8}}$
$\displaystyle=\frac{2m_{T}}{(2J_{T}+1)v^{2}}\sum_{\tau\tau^{\prime}}\frac{c_{8}^{\tau}c_{8}^{\tau^{\prime}}}{4}\left[v^{2}W_{M}^{\tau\tau^{\prime}}(q^{2})-v_{\rm
min}^{2}\left(W_{M}^{\tau\tau^{\prime}}(q^{2})-4\frac{\mu^{2}_{T}}{m^{2}_{N}}W_{\Delta}^{\tau\tau^{\prime}}(q^{2})\right)\right]\,,$
$\displaystyle\equiv c(q^{2})\left(1-\frac{v_{\rm
min}^{2}}{v^{2}}\right)+d(q^{2})\frac{v_{\rm min}^{2}}{v^{2}}$
$\displaystyle\frac{{\rm d}\sigma_{T}}{{\rm
d}E_{R}}\Bigg{|}_{\mathcal{O}_{11}}$
$\displaystyle=\frac{2m_{T}}{(2J_{T}+1)v^{2}}\sum_{\tau\tau^{\prime}}\frac{c_{11}^{\tau}c_{11}^{\tau^{\prime}}}{4}\frac{q^{2}}{m^{2}_{N}}W_{M}^{\tau\tau^{\prime}}(q^{2})\equiv
e(q^{2})\frac{v_{\rm min}^{2}}{v^{2}}\,,$ (2.7)

where we assumed that the DM particle has spin $J_{\chi}=1/2$. In Eq. (2.7),
we also emphasised the dependence on $v$ by introducing the functions $a$,
$b$, $c$, $d$, and $e$ which are implicitly defined via the above equation. We
refer to [23, 24] for further details on the effective theory of DM-nucleon
interactions and to [27, 28] for a discussion on its limitations and
extensions. Here, we just mention that the $\mathcal{O}_{1}$ interaction is
the familiar spin-independent interaction [24] and that the interactions
$\mathcal{O}_{7}$, $\mathcal{O}_{8}$ and $\mathcal{O}_{11}$ arise from the
non-relativistic reduction of simplified models. They can be the leading DM-
nucleon interactions in specific ranges for the DM mass [29, 30, 31, 32].

While in the non-relativistic effective theory of DM-nucleon interactions the
coupling constants $c_{j}^{\tau}$ are considered as independent, in concrete
models one generically expects to generate specific combinations of operators
in Tab. 1. In particular, this applies to the DM magnetic dipole model. The
model assumes that DM is made of spin 1/2 particles and is characterised by
the interaction Lagrangian
$\mathscr{L}=\frac{1}{2}\lambda_{\chi}\overline{\chi}\sigma^{\mu\nu}\chi
F_{\mu\nu}$, where $\chi$ is the DM spinor and $F_{\mu\nu}$ the
electromagnetic field strength tensor. Here, the coupling constant
$\lambda_{\chi}$ has dimension GeV-1. From this Lagrangian one finds the non-
relativistic amplitude for DM-nucleon scattering [33],

$\displaystyle\mathcal{M}$
$\displaystyle=e\lambda_{\chi}\left[2m_{N}Q_{N}\langle\mathcal{O}_{1}\rangle+4m_{\chi}g_{N}\langle\mathcal{O}_{4}\rangle+\frac{8m^{2}_{N}m_{\chi}Q_{N}}{|\mathbf{q}|^{2}}\langle\mathcal{O}_{5}\rangle-\frac{4m^{2}_{N}m_{\chi}g_{N}}{|\mathbf{q}|^{2}}\langle\mathcal{O}_{6}\rangle\right]\,,$
(2.8)

where $Q_{N}$, $N=p,n$, is the nucleon electric charge (0 for neutrons and 1
for protons), $g_{N}$, $N=p,n$, is the nucleon $g$-factor, and angle brackets
denote matrix elements between two-component nucleon and DM spinors. Eq. (2.8)
shows that the amplitude for DM-nucleon scattering in the DM magnetic dipole
model is a linear combination of matrix elements of the operators in Tab. 1
where some of the coefficients are constant while others scale like
$1/|\mathbf{q}|^{2}$ with the momentum transfer. Consequently, if one promotes
the coupling constants $c_{j}^{\tau}$ to functions of the momentum transfer,
the DM magnetic dipole model arises from a specific combination of operators
in Tab. 1. The associated cross section for DM-nucleus scattering reads as
follows,

$\displaystyle\frac{{\rm d}\sigma_{T}}{{\rm d}E_{R}}$
$\displaystyle=\frac{8m_{T}}{(2J_{T}+1)v^{2}}\alpha\pi\lambda_{\chi}^{2}\Bigg{[}\Bigg{(}\frac{1}{m_{\chi}^{2}}-\frac{1}{\mu_{T}^{2}}+\frac{1}{\mu_{T}^{2}}\frac{v^{2}}{v^{2}_{\rm
min}}\Bigg{)}W_{M}^{pp}+\frac{1}{m_{N}^{2}}\Bigg{(}\tilde{\mu}_{p}^{2}W_{\Sigma^{\prime}}^{pp}+2\tilde{\mu}_{p}\tilde{\mu}_{n}W_{\Sigma^{\prime}}^{np}$
$\displaystyle+\tilde{\mu}_{n}^{2}W_{\Sigma^{\prime}}^{nn}\Bigg{)}+4W_{\Delta}^{pp}-4\tilde{\mu}_{p}W_{\Delta\Sigma^{\prime}}^{pp}-4\tilde{\mu}_{n}W_{\Delta\Sigma^{\prime}}^{pn}\Bigg{]}\,,$
(2.9)

where $\alpha$ is the fine structure constant, while
$\tilde{\mu}_{N}=g_{N}/2$, with $\tilde{\mu}_{p}=2.8$ and
$\tilde{\mu}_{n}=-1.9$, is the dimensionless magnetic moment of the nucleon
[34]. The nuclear response functions in Eq. (2.9) are related to the ones in
Eq. (2.7) by $W_{M}^{pp}=W_{M}^{00}+2W_{M}^{01}+W_{M}^{11}$ and analogous
expressions [24].

We conclude this section by specifying our choices of nuclear targets and
response functions. So far, the CRESST experiment has exploited CaWO4, Li2MoO4
and Al2O3 targets. Here, we primarily focus on detectors employing Al2O3
crystals. This choice is motivated by the fact that Al2O3 is made of two light
elements, aluminium and oxygen, and Al2O3 detectors can operate with energy
thresholds as low as 19 eV [1]. Both properties make Al2O3 a promising target
to search for light DM particles of mass of a few GeV. Furthermore, aluminium
has spin 5/2 and the corresponding $W_{\Sigma^{\prime}}^{\tau\tau^{\prime}}$,
$W_{\Delta}^{\tau\tau^{\prime}}$ and
$W_{\Delta\Sigma^{\prime}}^{\tau\tau^{\prime}}$ functions are different from
zero. Consequently, Al2O3 detectors can be used to probe a wide range of DM-
nucleon interactions. In the calculations reported in Sec. 3.1, Sec. 3.2 and
Sec. 3.3, as well as in Secs. 5.1 and 5.2, we use the nuclear form factors of
[35] for aluminium and oxygen (as well as calcium in some investigations). For
the same elements, we use Helm form factors when treating the background
amplitude as a nuisance parameter (see Sec. 4), as they are faster to
evaluate. Results obtained by using Helm form factors lead to overestimate the
expected number of DM-induced nuclear recoil events by about 10%. Finally,
when comparing the timing information available to a Al2O3 multi-target
experiment with the one accessible with a fluorine, germanium, iodine, or
xenon detector, for the latter we use the nuclear form factors provided by the
DMFormFactor code [36].

##  3 Timing information in multi-target detectors

In this section, we critically review the timing information that is available
to low-threshold multi-target DM direct detection experiments when DM couples
to nuclei via one of the interactions in Tab. 1. This timing information is
contained in the coefficients of the Fourier series expansion of the
differential rate of nuclear recoil events in Eq. (2.1). The time-averaged
rate, $A_{0}$, and the annual modulation amplitude, $A_{1}$, are the first two
coefficients in this expansion. Starting from this observation, in Sec. 3.1 we
characterise the interactions in Tab. 1 in terms of the predicted
$A_{1}(v_{\rm min})$ curve. In Sec. 3.2, we show that most of the timing
information is contained in $A_{0}$ and $A_{1}$ for all interactions in Tab.
1, and hence higher-order harmonics can be neglected. Finally, in Sec. 3.3, we
characterise the interactions in Tab. 1 in terms of the predicted $t_{\rm
max}(v_{\rm min})$ curve introduced in [18], where $t_{\rm max}$ is the time
of maximum nuclear recoil rate. At the end of Sec. 3.3, we also comment on the
expected amplitude vs $v_{\rm min}$ and $t_{\rm max}$ vs $v_{\rm min}$ curves
in the case of DM-nucleus magnetic dipole interactions. The timing information
that we review in this section will then be used in the following two sections
to assess the prospects for DM signal discovery and model selection with next-
generation low-threshold detectors.

###  3.1 Annual modulation amplitude

Since the vectors $\boldsymbol{V_{\oplus}}$ and $\boldsymbol{\hat{r}_{s}}$ are
periodic functions of $t$ with period of one year, the rate of DM-nucleus
scattering events, Eq. (2.1), can be expanded in Fourier series around a
reference time $t_{0}>0$,

$\frac{{\rm d}R}{{\rm d}E_{R}}=A_{0}+\sum_{n=1}^{\infty}A_{n}\cos
n\omega(t-t_{0})+\sum_{n=1}^{\infty}B_{n}\sin n\omega(t-t_{0})\,,$ (3.1)

where $\omega=2\pi/\textrm{year}$ and for $t_{0}$ we choose the value of $t$
that maximises the vector $\boldsymbol{v_{\odot}}+\boldsymbol{V_{\oplus}}$.
The expansion coefficients $A_{0}$, $A_{n}$ and $B_{n}$ depend on the nuclear
recoil energy or, equivalently, on $v_{\rm min}$. They also depend on the
differential cross section for DM-nucleus scattering, which in turn depends on
how DM couples to nuclei. If the DM velocity distribution in the Galactic
reference frame is isotropic, then $B_{n}=0$ [37]. Furthermore, when
$f(\boldsymbol{v},t)$ changes slowly over velocity variations of the order of
$|\boldsymbol{V}_{\oplus}|$, then $A_{0}\gg A_{1}\gg A_{n\geq 2}$ [22]. The
existence of this hierarchy is well-known in the case of the familiar spin-
independent and spin-dependent interactions. In Sec. 3.2, we show that it also
applies to all DM-nucleon interactions considered in Tab. 1. In the “single
cosine approximation”, all expansion coefficients in Eq. (3.1) but $A_{0}$ and
$A_{1}$ are set to zero. Within this approximation, the modulation amplitude,
$A_{1}$, can be written as

$\displaystyle A_{1}(v_{\rm min})=\frac{1}{2}\left[\frac{{\rm d}R}{{\rm
d}E_{R}}\left(v_{\rm min},t_{0}\right)-\frac{{\rm d}R}{{\rm
d}E_{R}}\left(v_{\rm min},t_{0}+\mbox{year}/2\right)\right]\,,$ (3.2)

where we emphasised the dependence of the differential rate ${\rm d}R/{\rm
d}E_{R}$ on $v_{\rm min}$ and the time of the year.

Figure 1: Modulation amplitude $A_{1}$ as a function of the minimum velocity
$v_{\rm min}$ for $m_{\chi}$=10 GeV and the four interactions
$\mathcal{O}_{1}$, $\mathcal{O}_{7}$, $\mathcal{O}_{8}$ and
$\mathcal{O}_{11}$. Different colours refers to distinct targets. For the
interactions $\mathcal{O}_{1}$ and $\mathcal{O}_{11}$ (top panels), $A_{1}$
changes sign for $v_{\rm min}$ around 200 km s-1. For the interactions
$\mathcal{O}_{7}$ and $\mathcal{O}_{8}$ (bottom panels) $A_{1}>0$ for every
$v_{\rm min}$. Here, we assume a positive reference time $t_{0}$ (see Eq.
(3.2)). Figure 2: Same as Fig. 1, but now for $m_{\chi}=100$ GeV.

We now characterise the interactions in Tab. 1 in terms of the predicted
$A_{1}(v_{\rm min})$ curve. Since the dependence on $v$ of the DM-nucleus
scattering cross section determines the time-dependence of the nuclear recoil
rate in Eq. (2.1), we can restrict our analysis to the four representative
interactions in Eq. (2.7). For the latter, the four panels in Fig. 1 show
$A_{1}$ as a function of $v_{\rm min}$. We set the DM mass to $m_{\chi}=10$
GeV and one coupling constant at the time to $1/m_{V}^{2}$, where $m_{V}=246$
GeV. Each panel contains seven lines of different colours corresponding to
distinct target materials: oxygen, aluminium, sodium, calcium, germanium,
iodine and xenon. In our calculations, we take the isotopic abundance of the
different elements into account. For the interactions $\mathcal{O}_{1}$ and
$\mathcal{O}_{11}$ (top panels), the modulation amplitude changes sign for
$v_{\rm min}$ around 200 km s-1. This is consistent with having assumed
$t_{0}>0$ and $A_{1}$ unconstrained. On the other hand, we could have assumed
$A_{1}>0$ and $t_{0}$ unconstrained. This second convention would have implied
a change of sign for $t_{0}$, i.e. an “inversion of phase”, for $v_{\rm
min}\sim$ 200 km s-1. For the interactions $\mathcal{O}_{7}$ and
$\mathcal{O}_{8}$ (bottom panels) there is no inversion of phase, i.e.
$A_{1}>0$ for every $v_{\rm min}$. Fig. 1 also shows that for the
$\mathcal{O}_{1}$, $\mathcal{O}_{8}$ and $\mathcal{O}_{11}$ interactions, the
heavier the target the larger $A_{1}$. This is expected, as the corresponding
scattering cross section depends on $W_{M}^{\tau\tau^{\prime}}$, and,
therefore, scales with the number of nucleons squared in the small momentum
transfer limit. In contrast, for the spin-dependent $\mathcal{O}_{7}$
interaction, aluminium exhibits the largest modulation amplitude. Fig. 2 shows
$A_{1}$ as a function of $v_{\rm min}$ for the same interactions and targets
as Fig. 1, but now for $m_{\chi}=100$ GeV. Changing the DM particle mass, our
results remain qualitatively the same. Both Fig. 1 and Fig. 2 assume an
exposure of 1 kg$\times$day.

Figure 3: Ratios $|A_{1}|/A_{0}$ (left panel) and $|A_{2}|/|A_{1}|$ (right
panel) as a function of $v_{\rm min}$. Red dot-dashed lines correspond to a DM
particle of mass $50$ GeV scattering on a xenon target via $\mathcal{O}_{1}$.
Black solid lines refers to a DM particle of mass 10 GeV scattering on an
aluminium target via $\mathcal{O}_{7}$. The hierarchy $A_{0}\gg A_{1}\gg
A_{2}$ is valid at all $v_{\rm min}$ for interactions where $A_{1}$ does not
change sign, like $\mathcal{O}_{7}$. It is valid at all $v_{\rm min}$ but
around 200 km s-1 for interactions where $A_{1}$ changes sign around 200 km
s-1, like $\mathcal{O}_{1}$.

###  3.2 Higher-order harmonics

As we have seen in the previous subsection, the effective theory of DM-nucleon
interactions predicts two families of modulation amplitude. In the first one,
$A_{1}$ changes sign as a function of $v_{\rm min}$. In the second one,
$A_{1}>0$ for every $v_{\rm min}$. In this subsection we investigate the
validity of the single cosine approximation, i.e. the hierarchy $A_{0}\gg
A_{1}\gg A_{n\geq 2}$, for the two families of interactions. We focus on
$\mathcal{O}_{1}$ and $\mathcal{O}_{7}$ as representatives of the first and
second family, respectively.

Fig. 3, left panel (right panel), shows the ratio $|A_{1}|/A_{0}$
($|A_{2}|/|A_{1}|$) as a function of $v_{\rm min}$. Here, we consider two
benchmark scenarios: 1) A DM particle with interactions of $\mathcal{O}_{1}$
type and mass of $50$ GeV scattering on a xenon target (red, dot-dashed line).
2) A DM particle of mass 10 GeV interacting via the $\mathcal{O}_{7}$
interaction with an aluminium target (black, solid line). The red lines in
Fig. 3 reproduce the results of [22], which we use to validate our
calculations. The black lines correspond to results obtained here and of
interest for a low-threshold detector operating Al2O3 crystals. The spike
around 200 km s-1 in the $|A_{1}|/A_{0}$ curve for $\mathcal{O}_{1}$ reflects
the change of sign of $A_{1}$ for this interaction. Similarly, the smooth
$|A_{1}|/A_{0}$ curve associated with $\mathcal{O}_{7}$ follows from the
constant sign of $A_{1}$ in this case. As far as the $|A_{2}|/|A_{1}|$ curve
is concerned, we find two zeros at, respectively, 160-180 km s-1 and 500-540
km s-1 for the $\mathcal{O}_{1}$ interaction, and at 415-420 km s-1 and
690-700 km s-1 for the $\mathcal{O}_{7}$ interaction. Consequently, we
conclude that the single cosine approximation is valid at all $v_{\rm min}$
for interactions where $A_{1}$ does not change sign, like $\mathcal{O}_{7}$.
We also find that for interactions like $\mathcal{O}_{1}$ the single cosine
approximation is valid at all $v_{\rm min}$ values but around 200 km s-1, as
expected. While we present our results focusing on two specific combinations
of model parameters, it is important to stress that at a given $v_{\rm min}$,
and for single interactions, the ratios $|A_{1}|/A_{0}$ and $|A_{2}|/|A_{1}|$
do not depend on DM particle and target mass, as these are either reabsorbed
in the definition of $v_{\rm min}$ or factored out in time-independent pre-
factors; see Eq. (3.5) below111This is not exactly true for the category
$\mathcal{O}_{8}$ ($\mathcal{O}_{5}$ ), since $\eta$ and $\tilde{\eta}$
present different target dependent coefficients. However, Fig. 4 shows that
the term proportional to $\tilde{\eta}$ dominates, therefore treating the
ratios between higher order harmonics as independent from model parameters is
a good approximation..

###  3.3 Target dependence of the time of maximum rate

We now characterise the interactions in Tab. 1 in terms of the predicted
$t_{\rm max}(v_{\rm min})$ curve, extending previous results to target nuclei
of interest for the CRESST experiments. We start by noticing that to evaluate
the rate in Eq. (2.1) with differential cross sections for DM-nucleus
scattering given in Eqs. (2.7) and (2.9), we have to compute one or two of the
integrals,

$\displaystyle\eta(v_{\rm min},t)$ $\displaystyle=\int_{v\geq v_{\rm min}}{\rm
d}^{3}v\,\frac{f(\boldsymbol{v},t)}{v}\,,$ (3.3)
$\displaystyle\tilde{\eta}(v_{\rm min},t)$ $\displaystyle=\int_{v\geq v_{\rm
min}}{\rm d}^{3}v\,vf(\boldsymbol{v},t)\,.$ (3.4)

In terms of the $\eta(v_{\rm min},t)$ and $\tilde{\eta}(v_{\rm min},t)$
functions, the expected rate of nuclear recoils can be written as

$\displaystyle\frac{{\rm d}R}{{\rm d}E_{R}}=\mathscr{A}(v_{\rm
min},m_{T})\,\eta(v_{\rm min},t)+\mathscr{B}(v_{\rm
min},m_{T})\,\tilde{\eta}(v_{\rm min},t)\,.$ (3.5)

When $\mathscr{A}(v_{\rm min},m_{T})\neq 0$ and $\mathscr{B}(v_{\rm
min},m_{T})=0$, as in the case of the familiar spin-independent interaction,
the time of maximum rate, $t_{\rm max}$, is entirely determined by
$\eta(v_{\rm min},t)$ and, as such, it is a universal, i.e. target-independent
function of $v_{\rm min}$. The same applies to models where
$\mathscr{B}(v_{\rm min},m_{T})\neq 0$ and $\mathscr{A}(v_{\rm min},m_{T})=0$.
However, if both $\mathscr{A}(v_{\rm min},m_{T})$ and $\mathscr{B}(v_{\rm
min},m_{T})$ are different from zero and depend on $m_{T}$, and if the two
terms in Eq. (3.5) have comparable size, then the function $t_{\rm max}(v_{\rm
min})$ can depend on the target the DM particle scatters on. In this latter
case, there might be a range of $v_{\rm min}$ where the rate of DM-induced
nuclear recoils observed by a multi-target detector depends on time via the
superposition of two annual modulations of different phase. Therefore, it is
important to understand whether the target dependence of $t_{\rm max}(v_{\rm
min})$ is a generic feature of DM-nucleus interactions, or it can only occur
in specific models. Here, we address this question within the effective theory
of DM-nucleon interactions. First, we compute $t_{\rm max}(v_{\rm min})$ for
single interactions in Tab. 1. Then we focus on the specific linear
combination of interaction operators corresponding to the DM magnetic dipole
interaction model. This latter calculation extends the results of Ref. [18] to
aluminium, calcium and oxygen.

Figure 4: Time of maximum, $t_{\rm max}$, for the differential rate in Eq.
(2.1) as a function of $v_{\rm min}$ for the interactions $\mathcal{O}_{1}$
and $\mathcal{O}_{5}$. In the former case Eq. (2.1) is proportional to the
$\eta$ function, in the latter to the $\tilde{\eta}$ function. In both cases,
we assume aluminium as a target material and a DM particle mass of 10 GeV.

Inspection of Eq. (2.7) and Eq. (LABEL:eq:R) shows that $\mathscr{A}(v_{\rm
min},m_{T})$ and $\mathscr{B}(v_{\rm min},m_{T})$ are simultaneously different
from zero for the $\mathcal{O}_{5}$ and $\mathcal{O}_{8}$ interactions only.
At the same time, by explicitly evaluating the $\mathscr{A}(v_{\rm
min},m_{T})$ and $\mathscr{B}(v_{\rm min},m_{T})$ functions we find that in
both cases the first term in Eq. (3.5) always dominates over the second one.
As a result, in models where DM couples to nuclei via one of the operators in
Tab. 1 at the time, $t_{\rm max}(v_{\rm min})$ is a universal, target-
independent curve. We also find that all operators with differential cross
section scaling like $a(q^{2})/v^{2}$ and $e(q^{2})v_{\rm min}^{2}/v^{2}$
predict a $t_{\rm max}(v_{\rm min})$ curve modulating from about January 1st
at small $v_{\rm min}$ (when gravitational focusing is included) to June 1st
at large $v_{\rm min}$. For these interactions $A_{1}$ changes sign at small
$v_{\rm min}$. For the remaining interactions, the curve $t_{\rm max}(v_{\rm
min})$ is flat and $A_{1}$ has the same sign for all $v_{\rm min}$ values.
Fig. 4, left panel, illustrates our results on $t_{\rm max}(v_{\rm min})$
focusing on the $\mathcal{O}_{1}$ and $\mathcal{O}_{5}$ interactions.

Figure 5: Left panels. Time of maximum differential rate of nuclear recoil
events, $t_{\rm max}$, as a function of $v_{\rm min}$ for the DM magnetic
dipole interaction model. The top panel assumes calcium as a target material,
whereas the bottom panel refers to aluminium. Right panels. Differential rate
at different $v_{\rm min}$ values for calcium (top panel) and aluminium
(bottom panel). We set $m_{\chi}$ to the values in the legends.

Let us now focus on the DM magnetic dipole interaction model. Here, we extend
the results of [18] to aluminium, calcium and oxygen, as these targets are of
interest for CRESST detectors. Fig. 5 shows the $t_{\rm max}(v_{\rm min})$
curve for calcium (top panels) and aluminium (bottom panels). The
corresponding curve for oxygen is identical to the one associated with
calcium. The right panels in Fig. 5 report differential nuclear recoil rates
as a function of time for different values of $v_{\rm min}$ and, in the case
of aluminium, for different DM particle masses. In the case of calcium, the
curve $t_{\rm max}(v_{\rm min})$ is approximately flat, as the term
proportional to the $\tilde{\eta}$ function in Eq. (3.5) dominates over the
term proportional to $\eta$. This is due to the fact that 16O and 40Ca have
spin zero, which implies $W_{\Sigma^{\prime}}=0$. As a result, the surviving
term in Eq. (2.9) is proportional to $W^{pp}_{M}$ for both the $\eta$ and the
$\tilde{\eta}$ contributions and the coefficient of $\tilde{\eta}$ is
proportional to $1/v_{min}^{2}$, which explains why the $\tilde{\eta}$
contribution dominates over the $\eta$ one at small $v_{min}$.. In contrast,
the curve $t_{\rm max}(v_{\rm min})$ has a minimum above 100 km s-1 in the
case of aluminium. This arises from the interplay of the $\eta$ and
$\tilde{\eta}$ terms in Eq. (3.5). However, this minimum is pronounced only
for $m_{\chi}>10$ GeV. Therefore, for magnetic dipole DM-nucleus interactions,
aluminium, calcium and oxygen are characterised by the same, approximately
flat $t_{\rm max}(v_{\rm min})$ curve for small $m_{\chi}$ values. On the
other hand, for the familiar spin-independent interaction, the $t_{\rm
max}(v_{\rm min})$ curve span a period of about five months. Consequently, it
is reasonable to expect that timing information on the observed nuclear
recoils can be used to discriminate spin-independent from magnetic dipole DM-
nucleus interactions. Since the phenomenology of the two models significantly
differs in the small $v_{\rm min}$ limit, low-threshold experiments should
then be optimal to discriminate between them. While we here focused on the DM
magnetic dipole, we verified that the same conclusions apply to the anapole DM
model, which predicts a similar ${\rm d}\sigma_{T}/{\rm d}E_{R}$.

##  4 Statistical framework

In this section, we introduce the statistical framework used to assess the
prospects for DM signal discovery and model selection via timing (and energy)
information in a low-threshold experiment. As anticipated, our goal is to
determine whether a time dependent analysis can improve the sensitivity of
CRESST in the presence of an unidentified background resembling the low-energy
events recently reported by the CRESST collaboration. Here, we assume that $n$
nuclear recoil events have been observed in a low-threshold detector. The data
sample is distributed over a signal region, $\Delta S$, that we divide into
$N$ bins. $\Delta S$ and the $N$ bins are one-dimensional energy intervals,
when we only analyse the energy spectrum of the nuclear recoils (1D analysis).
They are two-dimensional, energy$\times$time intervals, when we assume that
timing information on the nuclear recoils is also accessible (2D analysis).

###  4.1 Signal discovery

Focusing on DM signal discovery, we compare the background-only hypothesis
with an alternative, signal-plus-background hypothesis. In the background-only
hypothesis the nuclear recoil events observed in the detector can be explained
in terms of experimental backgrounds alone. In the signal-plus-background
hypothesis, the observed events include a DM signal contribution. We model the
energy spectrum of the background events by assuming that they are distributed
in energy as the energy spectrum observed in CRESST-III, which presents an
excess below 200 eV, as reported in [17] and described in the introduction. We
also assume that the background events are homogeneously distributed in time.
In our calculations, the number of background events, $\lambda_{b}$, in a pre-
defined signal region $\Delta S$ is an input. We set $\lambda_{b}$ to the
desired value by varying the overall background normalisation, $\mathscr{B}$.
This model for the experimental background is only taken as a reference for
the order of magnitude of the expected background level, as the energy
spectrum reported in [17] was observed using CaWO4, whereas here we consider
Al2O3 as target material. Specifically, for the experimental background we
assume the energy spectrum [38],

$\displaystyle\frac{{\rm d}N_{b}}{{\rm
d}E_{R}}=\mathscr{B}\left(p_{0}+p_{1}E_{R}+p_{2}e^{-E_{R}/p_{3}}\right)\,,$
(4.1)

where we set $p_{0}$, $p_{1}$, $p_{2}$ and $p_{3}$ to the best fit values we
found by fitting Eq. (4.1) with $\mathscr{B}=1$ to the data in [39] (file
‘C3P1_DetA_full.dat’), namely: $p_{0}=28$ keV-1, $p_{1}=-0.8$ keV-2,
$p_{2}=19776$ keV-1 and $p_{3}=0.0423$ keV. In the signal-plus-background
hypothesis, we model the DM signal contribution in terms of spin-independent
interactions. Specifically, we set the DM-nucleon scattering cross section,
the proton to neutron coupling ratio and the DM particle mass to two benchmark
set of values: 1) the best fit values found in a fit of the latest DAMA
results [40], namely $\sigma_{\rm SI}=2.67\times 10^{-38}$ cm2,
$r=c_{p}/c_{n}=-0.76$ and $m_{\chi}=11.17$ GeV, respectively [41]; and 2)
$\sigma_{\rm SI}=4\times 10^{-42}$ cm2, $r=c_{p}/c_{n}=1$ and $m_{\chi}=3$
GeV, respectively. The first benchmark allows us to test the relative impact
of energy and timing information on the discovery of a signal resembling the
one reported by DAMA. The second one corresponds to the 90% C.L. upper bound
on the DM-nucleon scattering cross section reported by the DarkSide-50
collaboration for $m_{\chi}=3$ GeV [42]. Notice that $m_{\chi}=3$ GeV is the
smallest DM mass that can induce values of $v_{\rm min}$ below $200$ km s-1
for recoil energies below $\simeq 200$ eV, and that below $200$ km s-1 the
modulation amplitude changes sign in the case of spin-independent
interactions.

Let us now denote by $s_{i}$ the expected number of DM signal events in the
$i$-th bin when the model parameters are set at one of the above benchmark
points. The expected number of signal events in the signal region $\Delta S$
is hence $\lambda_{s}=\sum_{i}s_{i}$. Here, we adopt the notation of [43] and
introduce the signal strength parameter, $\mu$. For $\mu=0$ there is no DM
signal contribution to the observed nuclear recoil events, while for $\mu=1$
the DM contribution to the $i$-th bin is $s_{i}$. Furthermore, we denote by
$b_{i}$ the expected number of background events in the $i$-th bin and
introduce the parameter $\theta$ (see Eq. (4.6)) to model the uncertainties in
the background amplitude. Both in the 1D analysis and in the 2D analysis, we
test the background-only hypothesis against the signal-plus-background
hypothesis by means of the test statistic [43],

$\displaystyle
q_{0}=\left\\{\begin{array}[]{ll}-2\ln\lambda&\quad\hat{\mu}\geq 0\\\
0&\quad\hat{\mu}<0\end{array}\right.$ (4.4)

where

$\displaystyle\lambda=\frac{\mathscr{L}(0,\hat{\hat{\theta}})}{\mathscr{L}(\hat{\mu},\hat{\theta})}\,,$
(4.5)

and $\mathscr{L}(\mu,\theta)$ is the likelihood function,

$\displaystyle\mathscr{L}(\mu,\theta)=\prod_{i=1}^{N}\frac{(\mu s_{i}+\theta
b_{i})^{n_{i}}}{n_{i}!}e^{-(\mu s_{i}+\theta b_{i})}\,.$ (4.6)

Here, $\hat{\hat{\theta}}$ is the value of $\theta$ that maximises
$\mathscr{L}$ when $\mu=0$, while $(\hat{\mu},\hat{\theta})$ is the maximum
likelihood estimator. We sample the total number of events (background plus
signal) in the $i$-th bin, $n_{i}$, from a Poisson distribution of mean
$\overline{\mu}s_{i}+\overline{\theta}b_{i}$, where $\overline{\mu}$ and
$\overline{\theta}$ are the hypothesised values for $\mu$ and $\theta$,
respectively. For example, when sampling $n_{i}$ under the background-only
hypothesis, $\overline{\mu}=0$ and $\overline{\theta}=1$. When sampling
$n_{i}$ under the signal-plus-background hypothesis, $\overline{\mu}=1$ and
$\overline{\theta}=1$. By repeatedly sampling $n_{i}$ under the background-
only hypothesis, we obtain the probability density function (pdf)
$f(q_{0}|0)$. Sampling $n_{i}$ under the signal-plus-background hypothesis we
obtain the pdf $f(q_{0}|1)$. From $f(q_{0}|0)$ and $f(q_{0}|1)$, we obtain the
$p$-value for signal discovery,

$\displaystyle p=\int_{q_{0}^{\rm med}}^{\infty}{\rm d}q_{0}\,f(q_{0}|0)\,,$
(4.7)

where $q_{0}^{\rm med}$ is the median of $f(q_{0}|1)$. In order to test the
importance of the nuisance parameter $\theta$, we also evaluate the $p$-value
in Eq. (4.7) for $\theta$ set to the constant value
$\theta=\overline{\theta}=1$, so that $\hat{\hat{\theta}}=\hat{\theta}=1$.

###  4.2 Model selection

Focusing on DM model selection, we compare a “null-hypothesis” where DM
couples to nuclei through the familiar spin-independent interaction, with an
alternative hypothesis where DM couples to nuclei via a magnetic dipole
moment. Within the null-hypothesis, we model the DM contribution to the
observed nuclear recoils by assuming $m_{\chi}=3$ GeV, $r=c_{p}/c_{n}=1$ and
$\sigma_{\rm SI}=4\times 10^{-42}$ cm2, as in one of the realisations of the
signal-plus-background hypothesis considered in Sec. 4.1. This set of
parameters predicts 2.2 counts in an experimental setup consisting of an Al2O3
detector with an exposure of $23$ g$\times$year, when considering the interval
$[0.012-5.4]$ keV as energy window. Within the alternative hypothesis, we
assume $m_{\chi}=3$ GeV and couplings set to $\sigma_{\rm MD}\equiv
4\alpha^{2}\lambda_{\chi}^{2}=4.72\times 10^{-41}$ cm2. This reference cross
section value predicts 2.2 counts for the above experimental setup. Within
both hypotheses, we model the experimental background as in Sec. 4.1. Finally,
we compare the two hypotheses by using the log-likelihood ratio test
statistic,

$\displaystyle q=-2\ln\widetilde{\lambda}$ (4.8)

where

$\displaystyle\widetilde{\lambda}=\frac{\mathscr{L}^{(0)}(\hat{\mu},\hat{\theta})}{\mathscr{L}^{(a)}(\hat{\mu},\hat{\theta})}$
(4.9)

and

$\displaystyle\mathscr{L}^{(0)}(\mu,\theta)$
$\displaystyle=\prod_{i=1}^{N}\frac{(\mu s^{0}_{i}+\theta
b_{i})^{n_{i}}}{n_{i}!}e^{-(\mu s^{0}_{i}+\theta b_{i})}\,,$
$\displaystyle\mathscr{L}^{(a)}(\mu,\theta)$
$\displaystyle=\prod_{i=1}^{N}\frac{(\mu s^{a}_{i}+\theta
b_{i})^{n_{i}}}{n_{i}!}e^{-(\mu s^{a}_{i}+\theta b_{i})}\,.$ (4.10)

Here, $s^{0}_{i}$ ($s^{a}_{i}$) is the number of expected signal events in the
$i$-th bin under the null-hypothesis (alternative hypothesis). Importantly,
$n_{i}$ denotes the same dataset in both lines of Eq. (4.10). Analogously to
the case of DM signal discovery, we sample $n_{i}$ from a Poisson distribution
of mean $\overline{\mu}\,\overline{s}_{i}+\overline{\theta}b_{i}$. Sampling
$n_{i}$ under the null hypothesis, i.e. $\overline{\mu}=1$,
$\overline{\theta}=1$ and $\overline{s}_{i}=s_{i}^{0}$, we obtain the pdf
$f_{0}(q|1)$. Similarly, we obtain the pdf $f_{a}(q|1)$ when sampling $n_{i}$
under the alternative hypothesis, i.e. $\overline{\mu}=1$,
$\overline{\theta}=1$ and $\overline{s}_{i}=s_{i}^{a}$. Here, we quantify the
ability of a low-threshold experiment to reject the null in favour of the
alternative hypothesis in terms of the $p$-value,

$\displaystyle p=\int_{q^{\rm med}}^{\infty}{\rm d}q\,f_{0}(q|1)\,,$ (4.11)

where $q^{\rm med}$ is the median of the pdf $f_{a}(q|1)$.

Figure 6: Predicted energy (left panel) and energy$\times$time (right panel)
distributions of DM-induced nuclear recoil events for our benchmark 1
scenario: $\sigma_{\rm SI}=2.67\times 10^{-38}$ cm2, $r=c_{p}/c_{n}=-0.76$ and
$m_{\chi}=11.17$ GeV. Figure 7: Probability density functions $f(q_{0}|0)$ and
$f(q_{0}|1)$ (left panel) and associated cumulative probability density
functions (right panel) for a simulation assuming no timing information,
$\lambda_{b}=10^{3}$ and our benchmark scenario 1. The $p$-value for DM
discovery associated with this specific example is $p=0.046$.

##  5 Results

By applying the statistical methods of Sec. 4 to the theoretical framework of
Sec. 2, we now report the result of our $p$-value calculations for DM signal
discovery, Sec. 5.1, and model selection, Sec. 5.2. Here, we assume that the
number of successes in observing a value of $q$ ($q_{0}$) above $q^{\rm med}$
($q^{\rm med}_{0}$) obeys a binomial distribution of variance
$\sigma^{2}=p\mathscr{N}(1-p)$, where $\mathscr{N}$ is the number of times we
Monte Carlo generate the sample $n_{i}$. Consequently, the upper bound on the
population error of our $p$-value estimates is given by
$1/\sqrt{4\mathscr{N}}$. Since we choose $p=0.1$ (i.e. 90% C.L.) as a
reference $p$-value, we set $\mathscr{N}$ to $10^{3}$. In the tables presented
below, $p\simeq 0$ indicates that the result of our calculations is comparable
or smaller than the associated error, $1/\sqrt{4\mathscr{N}}\simeq 0.016$.

###  5.1 Signal discovery

We start by testing the background-only hypothesis against the signal-plus-
background hypothesis. We model the experimental background via Eq. (4.1) and
the DM contribution to the observed nuclear recoils in terms of spin-
independent interactions. Here, we set the DM-nucleon scattering cross
section, the proton to neutron coupling ratio and the DM particle mass to:

  1. 1.

$\sigma_{\rm SI}=2.67\times 10^{-38}$ cm2, $r=c_{p}/c_{n}=-0.76$ and
$m_{\chi}=11.17$ GeV, respectively [41]. Furthermore, $\Delta S=$ [105 eV, 16
keV].

  2. 2.

$\sigma_{\rm SI}=4\times 10^{-42}$ cm2, $r_{p}/r_{n}=1$ and $m_{\chi}=3$ GeV,
respectively [42]. In addition, $\Delta S=$ [12 eV, 5.4 keV].

Our motivations for focusing on these two benchmarks have been discussed in
Sec. 3. The first one corresponds to the best fit to the latest DAMA results,
the second one to the 90% C.L. upper limit on the DM-nucleon scattering cross
section reported by the DarkSide-50 experiment. In both cases, we assume Al2O3
as target material and a minimum exposure of 23 g$\times$158 days. In the 1D
analyses nuclear recoils are distributed in 10 energy bins, while in the 2D
analyses they are distributed in 102 energy$\times$time bins. We test the
results using both linear and logarithmic binning size. We use the test
statistic in Eq. (4.4).

####  5.1.1 Benchmark 1

Focusing on the first benchmark scenario, we now present our results on the
ability of a low-threshold experiment to reject the background-only
hypothesis. For this scenario, Fig. 6 shows the predicted energy (left panel)
and energy$\times$time (right panel) distributions of DM-induced nuclear
recoil events. As expected from our analysis of $A_{1}(v_{\rm min})$ and
$t_{\rm max}(v_{\rm min})$ in Sec. 3, the 2D distribution peaks around
December/January at small energies and in June at high energies.

Using the information contained in Fig. 6, we first test the ability of a low-
threshold experiment to reject the background-only hypothesis when the
background parameter $\theta$ is set to the constant value
$\theta=\overline{\theta}=1$. We will treat $\theta$ as a nuisance parameter
later on. As an illustrative example of the geometric interpretation of the
$p$-value for signal discovery in Eq. (4.11), Fig. 7 shows the pdfs
$f(q_{0}|0)$ and $f(q_{0}|1)$ and the associated cumulative pdfs for a
simulation assuming no timing information, $\lambda_{b}=10^{3}$,
$\lambda_{s}\simeq 141$ and $\Delta S=$ [105 eV, 16 keV] as signal region. As
reported in the legend, the $p$-value associated with this specific example is
$p=0.046$.

mass | time | $\lambda_{b}$ | $\lambda_{s}$ | N inter | p-value 1D | p-value 2D
---|---|---|---|---|---|---
23 g | 1 yr | $10^{3}$ | $\simeq 141$ | $10^{3}$ | $\simeq 0$ | $\simeq 0$
′′ | ′′ | $10^{4}$ | ′′ | ′′ | 0.046 | 0.045
′′ | ′′ | $10^{5}$ | ′′ | ′′ | 0.284 | 0.247
2 x 23g | ′′ | ′′ | $283$ | ′′ | 0.161 | 0.163
3 x 23g | ′′ | ′′ | $424$ | ′′ | 0.055 | 0.077
Table 2: $p$-values for DM signal discovery in the case of a DAMA-like signal (benchmark scenario 1), known background and no nuisance parameters. The number of bins is set to 10 both for energy and time. The first and second column correspond to detector mass and time of data taking, respectively. The fifth column gives the number of independent test statistic evaluations in our Monte Carlo simulations. The $p$-values that we obtained in 1D (2D) analyses are reported in the sixth (seventh) column. mass | time | $\lambda_{b}$ | $\lambda_{s}$ | N inter | p-value 1D | p-value 2D
---|---|---|---|---|---|---
23 g | 158 days | $500$ | $\simeq 69$ | $10^{3}$ | 0.086 | 0.083
′′ | 1 yr | $10^{3}$ | $\simeq 159$ | ′′ | $\simeq 0$ | $\simeq 0$
′′ | ′′ | $10^{4}$ | ′′ | ′′ | 0.21 | 0.212
′′ | ′′ | $10^{5}$ | ′′ | ′′ | 0.398 | 0.415
23 g | 2 yr | $2\cdot 10^{3}$ | $\simeq 319$ | ′′ | $\simeq 0$ | $\simeq 0$
′′ | ′′ | $2\cdot 10^{4}$ | ′′ | ′′ | 0.127 | 0.14
′′ | ′′ | $2\cdot 10^{5}$ | ′′ | ′′ | 0.34 | 0.362
Table 3: Same as Tab. 2, but now assuming an unknown background amplitude as a
nuisance parameter.

Tab. 2 extends the results in Fig. 7 and shows the $p$-values we find when
$\theta$ is set to the constant value $\theta=\overline{\theta}=1$ for
different combinations of $\lambda_{b}$, the expected number of background
events, and experimental exposures. More specifically, we consider values of
$\lambda_{b}$ ranging from $10^{3}$ to $10^{5}$ and experimental exposures
varying from 23 g$\times$year to $3\times$23 g$\times$year. Importantly, in
Tab. 2 we report $p$-values for signal discovery that we obtain when timing
information is taken into account (2D analysis) and when timing information is
not available (1D analysis). As one can see from Tab. 2, timing information
has a marginal impact on the ability of a low-threshold experiment to reject
the background-only hypothesis when the background events are modelled as in
Eq. (4.1) and the observed DM signal resembles the one reported by DAMA.

We now present our results on the ability of a low-threshold experiment to
reject the background-only hypothesis when the background parameter $\theta$
is considered as a nuisance parameter. Tab. 3 shows the $p$-values for signal
discovery that we find in simulations where $\lambda_{b}$ varies from $500$ to
$2\times 10^{5}$, the experimental exposure varies from 23 g$\times$158 days
(corresponding to the exposure reached in the CRESST-III analysis [17]) to 23
g$\times$2 years and $\theta$ is treated as a nuisance parameter. As in the
case of Tab. 2, Tab. 3 shows the results we find when timing information on
the observed nuclear recoils is taken into account, and when the latter is not
available. Comparing Tab. 2 with 3, one can see that in all cases the
$p$-values we find when $\theta$ is a nuisance parameters are larger than when
the experimental background is assumed to be known. This result is expected,
as the Neyman-Pearson lemma (see [43] and references therein) states that the
likelihood ratio in Eq. (4.4) for $\theta=\overline{\theta}=1$ is the test
statistic with the highest significance power, i.e. the highest probability of
rejecting the background-only hypothesis if the background-plus-signal
hypothesis is true. Tab. 3 also shows that timing information has a marginal
impact on the ability of a low-threshold experiment to reject the background-
only hypothesis in favour of benchmark 1 in the case of an unknown background
amplitude.

mass | time | $\lambda_{s}$ | $\lambda_{b}$ | N iter | p-value 1D | p-value 2d
---|---|---|---|---|---|---
23 g | 158 days | $\simeq 1$ | 500 | $10^{3}$ | 0.436 | 0.482
230 g | ′′ | $\simeq 22$ | $10^{4}$ | ′′ | 0.323 | 0.346
′′ | 2 yr | $\simeq 45$ | $2\cdot 10^{4}$ | ′′ | 0.223 | 0.263
′′ | 5 yr | $\simeq 112$ | $5\cdot 10^{4}$ | ′′ | 0.151 | 0.151
1 kg | ′′ | $\simeq 491$ | ′′ | ′′ | $\simeq 0$ | $\simeq 0$
′′ | ′′ | ′′ | $5\cdot 10^{5}$ | ′′ | $0.065$ | $0.094$
Table 4: Same as Tab. 3, but now modelling the DM-nucleon interaction as in
our benchmark 2 scenario.

####  5.1.2 Benchmark 2

Now we turn our attention to the second benchmark scenario considered in this
investigation, where $\sigma_{\rm SI}=4\times 10^{-42}$ cm2, $r=c_{p}/c_{n}=1$
and $m_{\chi}=3$ GeV. Focusing on this second scenario, below we present our
results on the ability of a low-threshold experiment to reject the background-
only hypothesis. We only consider the case of a background of unknown
amplitude and treat $\theta$ as a nuisance parameter. Tab. 4 shows the
$p$-values for signal discovery that we find comparing the background-only
hypothesis with a DM signal model corresponding to our benchmark 2. In our
analysis of benchmark 2, we let the number of background events,
$\lambda_{b}$, vary from $500$ to $5\times 10^{5}$ and consider experimental
exposures varying from 23 g$\times$158 days to 1 kg$\times$5 years.

Our results for benchmark 2 are summarised in Tab. 4. In the 1D analysis, the
$p$-values we obtain range from about 0.44, for $\lambda_{b}=500$ and an
exposure of 23 g$\times$158 days, to about 0.06, for $\lambda_{b}=5\times
10^{5}$ and an exposure of 1 kg$\times$5 years. Consistently with benchmark 2
being about two orders of magnitude below the current CRESST-III exclusion
limit on $\sigma_{\rm SI}$ for $m_{\chi}=3$ GeV, we find that the CRESST
experiment must operate with an exposure about 100 larger than the the current
one, i.e. 23 g$\times$158 days, in order to achieve a 90% C.L. discovery of a
DM model like benchmark 2. Remarkably, in all cases investigated here we find
comparable $p$-values in the 1D analysis and in the 2D analysis, where timing
information is taken into account. We therefore conclude that timing
information has a marginal impact on the ability of a low-threshold experiment
to reject the background-only hypothesis in favour of benchmark 2. These
results have been verified to be stable both using the linear and the
logarithmic distributed binning size.

###  5.2 Model selection

We conclude this section by presenting our results on the ability of a low-
threshold experiment to reject spin-independent in favour of magnetic dipole
DM-nucleon interactions. We refer to the former case as “null-hypothesis” and
to the latter as “alternative hypothesis”. As anticipated, we model the null-
hypothesis by assuming $m_{\chi}=3$ GeV, $r_{p}/r_{n}=1$ and $\sigma_{\rm
SI}=4\times 10^{-42}$ cm2. This set of parameters predicts 2.2 counts in the
$[0.012-5.4]$ keV energy window for an Al2O3 detector with an exposure of $23$
g$\times$year. We model the alternative hypothesis by assuming $m_{\chi}=3$
GeV and a cross section $\sigma_{\rm MD}\equiv
4\alpha\lambda_{\chi}^{2}=4.72\times 10^{-41}$ cm2. Consistently, this
reference cross section value produces 2.2 counts in the same Al2O3 detector
and energy window.

Tab. 5 shows the $p$-values we find when comparing null and alternative
hypotheses. In this investigation, we consider experimental exposures ranging
from 23 g$\times$158 days to 1 kg$\times$5 years, and $\lambda_{b}$ varying
from $10^{3}$ to $2\times 10^{4}$. We report both $p$-values for the 1D
analysis and $p$-values for the 2D analysis, where timing information is taken
into account. In all cases, we treat the background parameter $\theta$ as a
nuisance parameter. The information in Tab. 5 is two-fold. On the one hand we
find that a 2D-analysis based on the energy and time distribution of the event
rate produces more significant results with respect to the 1D analysis, at
equal exposure. While previous results do not depend on the linear or
logarithmic binning method, the different significance of the 2D and 1D
analyses in this case strictly depends on the better accuracy of the
logarithmic method. On the other hand, our results show that a 90% C.L.
rejection of spin-independent interactions in favour of a DM magnetic dipole
coupling is feasible at the CRESST experiment with an exposure of 460
g$\times$year and $2\times 10^{3}$ background events.

mass | time | $\lambda_{b}$ | $\lambda_{s}$ | N inter | p-value 1D | p-value 2D
---|---|---|---|---|---|---
23 g | 1 yr | $10^{3}$ | $\simeq 2.2$ | $10^{3}$ | 0.486 | 0.442
230 g | 2 yr | $2\cdot 10^{3}$ | $\simeq 44$ | ′′ | 0.107 | 0.065
′′ | 5 yr | $5\cdot 10^{3}$ | $\simeq 110$ | ′′ | 0.024 | 0.005 ($\simeq 0$)
1 kg | ′′ | $2\cdot 10^{4}$ | $\simeq 478$ | ′′ | $\simeq 0$ | $\simeq 0$
Table 5: $p$-values for rejecting spin-independent interactions with
$m_{\chi}=3$ GeV and $\sigma_{\rm SI}=4\times 10^{-42}$ cm2 in favour of a
magnetic dipole coupling producing the same number of events in $\Delta S=$
[12 eV, 5.4 keV]. We report results with and without timing information,
obtained with logarithmic binning. The notation is the same as in previous
tables.

##  6 Summary and conclusions

We assessed the prospects for DM signal discovery and model selection via
timing information in a low-threshold experiment. In all calculations, we
assumed a time-independent experimental background with an energy spectrum
resembling the one of the low-energy events observed in [17]. We focused on
Al2O3 as a detector material, as it is composed of two relatively light
elements (oxygen and aluminium) and showed very good performances, with energy
thresholds as low as 19 eV [1].

As a preparatory step, we critically reviewed the timing information that is
available to low-threshold multi-target experiments, focusing on interaction-
dependent features, such as the shape of the $t_{\rm max}(v_{\rm min})$ curve,
which can depend on the target material. We also explored the significance of
higher-order harmonics in a Fourier series expansion of the predicted nuclear
recoil rate. We discussed these aspects within the general non-relativistic
effective theory of DM-nucleon interactions, and for specific interaction
models, such as the DM magnetic dipole coupling.

Focusing on DM signal discovery, we compared a background-only hypothesis with
an alternative, background-plus-signal hypothesis. Here, we modelled the
interaction between DM and nuclei in terms of spin-independent interactions
and considered two benchmark scenarios separately. The first scenario
corresponds to the best fit mass and couplings to the latest DAMA/LIBRA
results. The second one corresponds to the 90% C.L upper limit on the spin-
independent DM-nucleon scattering cross section reported by DarkSide-50 for a
DM particle mass of 3 GeV. Focusing on DM model selection, we compared a
“null-hypothesis” where DM couples to nuclei through the familiar spin-
independent interaction (and parameters as in the second scenario described
above) with an alternative hypothesis where DM couples to nuclei via a
magnetic dipole moment. We presented our results on the ability of a low-
threshold experiment to discover a signal above a background resembling the
excess found in the CRESST-III energy spectrum, or reject the spin-independent
interaction in favour of a magnetic dipole coupling in terms of $p$-values,
that we list in Tabs. 2, 3, 4 and 5. Importantly, we performed our $p$-value
calculations under two different assumptions: 1) Taking timing information
into account. 2) Assuming that the latter is not available.

Remarkably, we found that while timing information has a marginal impact on
the potential for DM signal discover, it can be beneficial for model selection
in a low-threshold experiment. This applies to all cases studied here and
reported in Tabs. 2, 3, 4 and 5. On the one hand, our results indicate that
for this class of experiments it might not be worth tackling the challenge of
operating cryogenic detectors under stable conditions for long periods, until
the goal remains observing the first dark matter signal. On the other hand, if
cryogenic experiments undertake such an effort, in case of detection of a
positive dark matter signal, the timing information could provide more
reliable results on the dark matter nature. For example, we found that the
$p$-value for rejecting spin-independent interactions in favour of a magnetic
dipole coupling is of about 0.1 when the exposure is 460 g$\times$year, while
about 0.06 if timing information is taken into accout. Our simulations are
limited to at most 5 years of data taking, specific benchmarks and to the
assumption of a background constant in time; it should be considered as a
starting point towards further, refined analyses of the impact of timing
information on the performance of future low-threshold experiments.

In addition, our results also show that for the set of model parameters
considered here a 90% C.L. rejection of spin-independent interactions in
favour of a magnetic dipole coupling is feasible with an upgrade of the CRESST
experiment [2], where timing information is not still available; see Tab. 5,
second row. In order to consolidate our findings, we compared spin-independent
and magnetic dipole interactions under the assumption of a detector efficiency
below 100%. Specifically, in our likelihood analysis we multiplied the
expected number of signal events in each bin, $s_{i}$, by a detector
efficiency sampled from a uniform distribution in the 50% - 70% range. For an
exposure of 1 kg$\times$2 years and a background level of $5\times 10^{3}$, we
found a $p$-value for model selection of 0.072. As expected, the experimental
sensitivity diminishes when detector efficiency is taken into account. At the
same time, we found that a 90% C.L. rejection of spin-independent interactions
in favour of a magnetic dipole coupling is feasible at a CRESST-III
experiment, even when the detector efficiency is below 100%.

## Acknowledgments

This work has been performed within a double doctoral degree agreement between
Chalmers University and Technology and the Gran Sasso Science Institute. The
authors thank, Federica Petricca, Paolo Gorla, Karoline Schäffner and
Francesco Vissani for stimulating and supporting the implementation of the
agreement between the two institutes, and Florian Reindl and the CRESST
Collaboration for the many valuable discussions and suggestions which inspired
the research contained in this work. RC acknowledges support from an
individual research grant from the Swedish Research Council, dnr. 2018-05029.
The results presented here made use of the computer programmes Wolfram
Mathematica [44] and DMFormFactor [36].

##  Appendix A DM response functions

$\displaystyle R_{M}^{\tau\tau^{\prime}}\left(v_{T}^{\perp 2},{q^{2}\over
m_{N}^{2}}\right)$ $\displaystyle=$ $\displaystyle
c_{1}^{\tau}c_{1}^{\tau^{\prime}}+{J_{\chi}(J_{\chi}+1)\over
3}\left[{q^{2}\over m_{N}^{2}}v_{T}^{\perp
2}c_{5}^{\tau}c_{5}^{\tau^{\prime}}+v_{T}^{\perp
2}c_{8}^{\tau}c_{8}^{\tau^{\prime}}+{q^{2}\over
m_{N}^{2}}c_{11}^{\tau}c_{11}^{\tau^{\prime}}\right]$ $\displaystyle
R_{\Phi^{\prime\prime}}^{\tau\tau^{\prime}}\left(v_{T}^{\perp 2},{q^{2}\over
m_{N}^{2}}\right)$ $\displaystyle=$ $\displaystyle{q^{2}\over
4m_{N}^{2}}c_{3}^{\tau}c_{3}^{\tau^{\prime}}+{J_{\chi}(J_{\chi}+1)\over
12}\left(c_{12}^{\tau}-{q^{2}\over
m_{N}^{2}}c_{15}^{\tau}\right)\left(c_{12}^{\tau^{\prime}}-{q^{2}\over
m_{N}^{2}}c_{15}^{\tau^{\prime}}\right)$ $\displaystyle
R_{\Phi^{\prime\prime}M}^{\tau\tau^{\prime}}\left(v_{T}^{\perp 2},{q^{2}\over
m_{N}^{2}}\right)$ $\displaystyle=$ $\displaystyle
c_{3}^{\tau}c_{1}^{\tau^{\prime}}+{J_{\chi}(J_{\chi}+1)\over
3}\left(c_{12}^{\tau}-{q^{2}\over
m_{N}^{2}}c_{15}^{\tau}\right)c_{11}^{\tau^{\prime}}$ $\displaystyle
R_{\tilde{\Phi}^{\prime}}^{\tau\tau^{\prime}}\left(v_{T}^{\perp 2},{q^{2}\over
m_{N}^{2}}\right)$ $\displaystyle=$ $\displaystyle{J_{\chi}(J_{\chi}+1)\over
12}\left[c_{12}^{\tau}c_{12}^{\tau^{\prime}}+{q^{2}\over
m_{N}^{2}}c_{13}^{\tau}c_{13}^{\tau^{\prime}}\right]$ $\displaystyle
R_{\Sigma^{\prime\prime}}^{\tau\tau^{\prime}}\left(v_{T}^{\perp 2},{q^{2}\over
m_{N}^{2}}\right)$ $\displaystyle=$ $\displaystyle{q^{2}\over
4m_{N}^{2}}c_{10}^{\tau}c_{10}^{\tau^{\prime}}+{J_{\chi}(J_{\chi}+1)\over
12}\left[c_{4}^{\tau}c_{4}^{\tau^{\prime}}+{q^{2}\over
m_{N}^{2}}(c_{4}^{\tau}c_{6}^{\tau^{\prime}}+c_{6}^{\tau}c_{4}^{\tau^{\prime}})+{q^{4}\over
m_{N}^{4}}c_{6}^{\tau}c_{6}^{\tau^{\prime}}\right.$ $\displaystyle+$
$\displaystyle\left.v_{T}^{\perp
2}c_{12}^{\tau}c_{12}^{\tau^{\prime}}+{q^{2}\over m_{N}^{2}}v_{T}^{\perp
2}c_{13}^{\tau}c_{13}^{\tau^{\prime}}\right]$ $\displaystyle
R_{\Sigma^{\prime}}^{\tau\tau^{\prime}}\left(v_{T}^{\perp 2},{q^{2}\over
m_{N}^{2}}\right)$ $\displaystyle=$ $\displaystyle{1\over 8}\left[{q^{2}\over
m_{N}^{2}}v_{T}^{\perp 2}c_{3}^{\tau}c_{3}^{\tau^{\prime}}+v_{T}^{\perp
2}c_{7}^{\tau}c_{7}^{\tau^{\prime}}\right]{J_{\chi}(J_{\chi}+1)\over
12}\left[c_{4}^{\tau}c_{4}^{\tau^{\prime}}+{q^{2}\over
m_{N}^{2}}c_{9}^{\tau}c_{9}^{\tau^{\prime}}\right.$ $\displaystyle+$
$\displaystyle\left.{v_{T}^{\perp 2}\over 2}\left(c_{12}^{\tau}-{q^{2}\over
m_{N}^{2}}c_{15}^{\tau}\right)\left(c_{12}^{\tau^{\prime}}-{q^{2}\over
m_{N}^{2}}c_{15}^{\tau\prime}\right)+{q^{2}\over 2m_{N}^{2}}v_{T}^{\perp
2}c_{14}^{\tau}c_{14}^{\tau^{\prime}}\right]$ $\displaystyle
R_{\Delta}^{\tau\tau^{\prime}}\left(v_{T}^{\perp 2},{q^{2}\over
m_{N}^{2}}\right)$ $\displaystyle=$ $\displaystyle{J_{\chi}(J_{\chi}+1)\over
3}\left[{q^{2}\over
m_{N}^{2}}c_{5}^{\tau}c_{5}^{\tau^{\prime}}+c_{8}^{\tau}c_{8}^{\tau^{\prime}}\right]$
$\displaystyle R_{\Delta\Sigma^{\prime}}^{\tau\tau^{\prime}}\left(v_{T}^{\perp
2},{q^{2}\over m_{N}^{2}}\right)$ $\displaystyle=$
$\displaystyle{J_{\chi}(J_{\chi}+1)\over
3}\left[c_{5}^{\tau}c_{4}^{\tau^{\prime}}-c_{8}^{\tau}c_{9}^{\tau^{\prime}}\right],$

where $J_{\chi}$ is the DM particle spin. In all numerical application in the
work, for DM we assume $J_{\chi}=1/2$.

## References

  * [1] CRESST collaboration, G. Angloher et al., _Results on MeV-scale dark matter from a gram-scale cryogenic calorimeter operated above ground_ , _Eur. Phys. J. C_ 77 (2017) 637 [1707.06749].
  * [2] CRESST collaboration, G. Angloher et al., _Probing low WIMP masses with the next generation of CRESST detector_ , 1503.08065.
  * [3] G. Bertone and D. Hooper, _History of dark matter_ , _Rev. Mod. Phys._ 90 (2018) 045002 [1605.04909].
  * [4] A. Drukier and L. Stodolsky, _Principles and Applications of a Neutral Current Detector for Neutrino Physics and Astronomy_ , _Phys. Rev._ D30 (1984) 2295.
  * [5] M. W. Goodman and E. Witten, _Detectability of Certain Dark Matter Candidates_ , _Phys. Rev._ D31 (1985) 3059.
  * [6] T. Marrodán Undagoitia and L. Rauch, _Dark matter direct-detection experiments_ , _J. Phys._ G43 (2016) 013001 [1509.08767].
  * [7] Particle Data Group collaboration, P. A. Zyla et al., _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * [8] R. Trotta, R. Ruiz de Austri and L. Roszkowski, _Prospects for direct dark matter detection in the Constrained MSSM_ , _New Astron. Rev._ 51 (2007) 316 [astro-ph/0609126].
  * [9] R. Trotta, R. Ruiz de Austri and L. Roszkowski, _Direct dark matter detection around the corner? Prospects in the constrained MSSM_ , _J. Phys. Conf. Ser._ 60 (2007) 259.
  * [10] D. G. Cerdeno, T. Kobayashi and C. Munoz, _Prospects for the direct detection of neutralino dark matter in orbifold scenarios_ , _JHEP_ 01 (2008) 009 [0709.0858].
  * [11] D. S. Alves, M. Lisanti and J. G. Wacker, _Poker face of inelastic dark matter: Prospects at upcoming direct detection experiments_ , _Phys. Rev. D_ 82 (2010) 031901 [1005.5421].
  * [12] R. Catena, _Prospects for direct detection of dark matter in an effective theory approach_ , _JCAP_ 07 (2014) 055 [1406.0524].
  * [13] E. Del Nobile, G. B. Gelmini and S. J. Witte, _Prospects for detection of target-dependent annual modulation in direct dark matter searches_ , _JCAP_ 02 (2016) 009 [1512.03961].
  * [14] P. Sandick, K. Sinha and F. Teng, _Simplified Dark Matter Models with Charged Mediators: Prospects for Direct Detection_ , _JHEP_ 10 (2016) 018 [1608.00642].
  * [15] B. J. Kavanagh, F. S. Queiroz, W. Rodejohann and C. E. Yaguna, _Prospects for determining the particle/antiparticle nature of WIMP dark matter with direct detection experiments_ , _JHEP_ 10 (2017) 059 [1706.07819].
  * [16] S. J. Witte, V. Gluscevic and S. D. McDermott, _Prospects for Distinguishing Dark Matter Models Using Annual Modulation_ , _JCAP_ 02 (2017) 044 [1612.07808].
  * [17] CRESST collaboration, A. Abdelhameed et al., _First results from the CRESST-III low-mass dark matter program_ , _Phys. Rev. D_ 100 (2019) 102002 [1904.00498].
  * [18] E. Del Nobile, G. B. Gelmini and S. J. Witte, _Target dependence of the annual modulation in direct dark matter searches_ , _Phys. Rev. D_ 91 (2015) 121302 [1504.06772].
  * [19] M. S. Alenazi and P. Gondolo, _Phase-space distribution of unbound dark matter near the Sun_ , _Phys. Rev. D_ 74 (2006) 083518 [astro-ph/0608390].
  * [20] S. K. Lee, M. Lisanti, A. H. G. Peter and B. R. Safdi, _Effect of Gravitational Focusing on Annual Modulation in Dark-Matter Direct-Detection Experiments_ , _Phys. Rev. Lett._ 112 (2014) 011301 [1308.1953].
  * [21] C. McCabe, _The Earth’s velocity for direct detection experiments_ , _JCAP_ 1402 (2014) 027 [1312.1355].
  * [22] S. K. Lee, M. Lisanti and B. R. Safdi, _Dark-Matter Harmonics Beyond Annual Modulation_ , _JCAP_ 11 (2013) 033 [1307.5323].
  * [23] J. Fan, M. Reece and L.-T. Wang, _Non-relativistic effective theory of dark matter direct detection_ , _JCAP_ 1011 (2010) 042 [1008.1591].
  * [24] A. L. Fitzpatrick, W. Haxton, E. Katz, N. Lubbers and Y. Xu, _The Effective Field Theory of Dark Matter Direct Detection_ , _JCAP_ 1302 (2013) 004 [1203.3542].
  * [25] P. Gondolo, S. Kang, S. Scopel and G. Tomar, _The effective theory of nuclear scattering for a WIMP of arbitrary spin_ , 2008.05120.
  * [26] D. Gazda, R. Catena and C. Forssén, _Ab initio nuclear response functions for dark matter searches_ , _Phys. Rev._ D95 (2017) 103011 [1612.09165].
  * [27] M. Hoferichter, P. Klos and A. Schwenk, _Chiral power counting of one- and two-body currents in direct detection of dark matter_ , _Phys. Lett._ B746 (2015) 410 [1503.04811].
  * [28] F. Bishara, J. Brod, B. Grinstein and J. Zupan, _Chiral Effective Theory of Dark Matter Direct Detection_ , _JCAP_ 1702 (2017) 009 [1611.00368].
  * [29] R. Catena, J. Conrad and M. B. Krauss, _Compatibility of a dark matter discovery at XENONnT or LZ with the WIMP thermal production mechanism_ , _Phys. Rev._ D97 (2018) 103002 [1712.07969].
  * [30] S. Baum, R. Catena, J. Conrad, K. Freese and M. B. Krauss, _Determining Dark Matter properties with a XENONnT/LZ signal and LHC-Run3 mono-jet searches_ , _Phys. Rev._ D97 (2018) 083002 [1709.06051].
  * [31] S. Baum, R. Catena and M. B. Krauss, _Impact of a XENONnT Signal on LHC Dijet Searches_ , 1812.01594.
  * [32] S. Baum, R. Catena and M. B. Krauss, _Constraints on Simplified Models for Dark Matter from LHC Dijet Searches_ , 1812.01585.
  * [33] B. J. Kavanagh, P. Panci and R. Ziegler, _Faint Light from Dark Matter: Classifying and Constraining Dark Matter-Photon Effective Operators_ , _JHEP_ 04 (2019) 089 [1810.00033].
  * [34] M. I. Gresham and K. M. Zurek, _On the Effect of Nuclear Response Functions in Dark Matter Direct Detection_ , 1401.3739.
  * [35] R. Catena and B. Schwabe, _Form factors for dark matter capture by the Sun in effective theories_ , _JCAP_ 1504 (2015) 042 [1501.03729].
  * [36] N. Anand, A. L. Fitzpatrick and W. C. Haxton, _Weakly interacting massive particle-nucleus elastic scattering response_ , _Phys. Rev._ C89 (2014) 065501 [1308.6288].
  * [37] K. Freese, M. Lisanti and C. Savage, _Annual Modulation of Dark Matter: A Review_ , _Rev. Mod. Phys._ 85 (2013) 1561 [1209.3339].
  * [38] D. Schmiedmayer, _Calculation of dark-matter exclusions-limits using a maximum likelihood approach_ , Master’s thesis, Wien, 2019.
  * [39] CRESST collaboration, A. Abdelhameed et al., _Description of CRESST-III Data_ , 1905.07335.
  * [40] R. Bernabei et al., _First model independent results from DAMA/LIBRA-phase2_ , _Nucl. Phys. Atom. Energy_ 19 (2018) 307 [1805.10486].
  * [41] S. Kang, S. Scopel, G. Tomar and J.-H. Yoon, _DAMA/LIBRA-phase2 in WIMP effective models_ , _JCAP_ 07 (2018) 016 [1804.07528].
  * [42] DarkSide collaboration, P. Agnes et al., _Low-Mass Dark Matter Search with the DarkSide-50 Experiment_ , _Phys. Rev. Lett._ 121 (2018) 081307 [1802.06994].
  * [43] G. Cowan, K. Cranmer, E. Gross and O. Vitells, _Asymptotic formulae for likelihood-based tests of new physics_ , _Eur. Phys. J._ C71 (2011) 1554 [1007.1727].
  * [44] Wolfram Research, Inc., _Mathematica v12.0_ , 2019.