{'fraction_non_alphanumeric': 0.0400302909190349, 'fraction_numerical': 0.021392961568396474, 'mean_word_length': 4.1797777293528, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Prospects for Dark Matter signal discovery and model selection via timing information in a low-threshold experiment

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 (SI) interaction in favour of a magnetic dipole (MD) 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 SI-interactions in favour of a MD-coupling is about 0.11 when the experimental exposure is 460 g x 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 x year. At the same time, our results show that a 90% C.L. rejection of SI-interactions in favour of a MD-coupling is within reach of an upgrade of the CRESST experiment [2].


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×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.

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, is the minimum velocity a DM particle must have to deposit an energy E R in the detector, m T and µ T are the target nucleus and DM-nucleus reduced mass, respectively, m χ is the DM particle mass and ρ χ is the local DM density. In multitarget detectors, like CRESST, the contribution to Eq. (2.1) of each element in the detector material is weighted by the corresponding mass fraction, ξ T . The differential rate of DMnucleus scattering events also depends on the differential cross section for DM-target nucleus scattering, dσ T /dE R , and on the local DM velocity distribution expressed as a function of the DM velocity in the detector reference frame, f (v, t). The velocity distribution f (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 (v, t) by taking the gravitational focusing of the Sun into account, which implies [19], wheref is the DM velocity distribution in the Galactic reference frame (where, by construction, the mean velocity of DM particles is zero), v (11, 232, 7) km s −1 is the relative velocity of the Sun with respect to the Galactic centre in Galactic coordinates, v is the DM particle velocity in the detector rest frame and V ⊕ (t) is the relative velocity of the Earth with respect to the Sun. Consistently, we denote by v s = v + V ⊕ the DM particle velocity in the Sun rest frame. With this notation, the solar frame velocity v ∞ a DM particle must have at infinity to move with velocity v s with respect to the Sun at the detector position is given by [19], where r s is Earth-Sun distance,r s is the time dependent unit vector that points from the Sun to the Earth, v 2 ∞ = v 2 − 2GM /r s , M is the mass of the Sun, and G is Newton's constant. Explicit expressions for V ⊕ andr s in Galactic coordinates can be found in [20] and [21,22], respectively. Forf , we assume a Gaussian distribution truncated at the escape velocity v esc ,f Table 1. Interaction operators defining the non-relativistic effective theory of spin 1/2 DM-nucleon interactions [23,24].

4)
where µ N is the DMnucleon reduced mass, is the transverse relative velocity and 1 χN is the identity in the DM-nucleon spin space.
whereṽ = v +v ∞ . Here, we assume the one-dimensional velocity dispersion σ v = 164 km s −1 and a local escape velocity of v esc , = 544 km s −1 . The normalisation factor, N esc , in Eq. (2.4) is given by, (2.5) The differential cross section for DM-nucleus scattering, dσ T /dE R encodes all particle and nuclear physics inputs required to evaluate Eq. (2.1). Within the non relativistic effective theory of DM-nucleon interactions, dσ T /dE R can be expressed as follows, 2m T E R is the momentum transfer and m N is the nucleon mass. At most eight DM and nuclear response functions, R τ τ k and W τ τ k , k = M, Σ , Σ , Φ , Φ M,Φ , ∆, ∆Σ , respectively, can appear in Eq. (2.6) [24]. The R τ τ k functions depend on q 2 /m 2 N , v ⊥2 T and on the coupling constants for DM-nucleon interactions, c τ j . 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 χ interaction types that are invariant under Galilean transformations and spatial rotations for a DM particle of spin J χ [25]. We list the 14 independent interaction operators arising for spin 1/2 DM in Tab. 1. The eight nuclear response functions W τ τ k 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 τ and τ run from 0 to 1: 0 corresponds to "isoscalar" interactions and 1 to "isovector" interactions [24]. In terms of, e.g. c τ 1 , the neutron to proton coupling ratio can be expressed as r = c p /c n = (c 0 1 − c 1 1 )/(c 0 1 + c 1 1 ). Finally, k = 0 for k = M, Σ , Σ , and k = 1 otherwise.
The way the cross section dσ T /dE R depends on v determines how the nuclear recoil rate depends on t.  [24] and here associated with the coupling constants c τ 1 , c τ 7 , c τ 8 and c τ 11 are representative of these four classes. The corresponding DM-nucleus scattering cross section reads as, where we assumed that the DM particle has spin J χ = 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 DMnucleon interactions and to [27,28] for a discussion on its limitations and extensions. Here, we just mention that the O 1 interaction is the familiar spin-independent interaction [24] and that the interactions O 7 , O 8 and 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 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 L = 1 2 λ χ χσ µν χF µν , where χ is the DM spinor and F µν the electromagnetic field strength tensor. Here, the coupling constant λ χ has dimension GeV −1 . From this Lagrangian one finds the non-relativistic amplitude for DM-nucleon scattering [33], 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 twocomponent 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/|q| 2 with the momentum transfer. Consequently, if one promotes the coupling constants c τ j 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, where α is the fine structure constant, whileμ N = g N /2, withμ p = 2.8 andμ 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 pp M = W 00 M + 2W 01 M + W 11 M 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 CaWO 4 , Li 2 MoO 4 and Al 2 O 3 targets. Here, we primarily focus on detectors employing Al 2 O 3 crystals. This choice is motivated by the fact that Al 2 O 3 is made of two light elements, aluminium and oxygen, and Al 2 O 3 detectors can operate with energy thresholds as low as 19 eV [1]. Both properties make Al 2 O 3 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 τ τ Σ , W τ τ ∆ and W τ τ ∆Σ functions are different from zero. Consequently, Al 2 O 3 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 Al 2 O 3 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].

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 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 max (v min ) curve introduced in [18], where t 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 min and t max vs v 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.

Annual modulation amplitude
Since the vectors V ⊕ andr 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 where ω = 2π/year and for t 0 we choose the value of t that maximises the vector v +V ⊕ . The expansion coefficients A 0 , A n and B n depend on the nuclear recoil energy or, equivalently, on v 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 (v, t) changes slowly over velocity variations of the order of |V ⊕ |, then A 0 A 1 A n≥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 where we emphasised the dependence of the differential rate dR/dE R on v min and the time of the year. We now characterise the interactions in Tab. 1 in terms of the predicted A 1 (v 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 min . We set the DM mass to m χ = 10 GeV and one coupling constant at the time to 1/m 2 V , 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 O 1 and O 11 (top panels), the modulation amplitude changes sign for v 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 min ∼ 200 km s −1 . For the interactions O 7 and O 8 (bottom panels) there is no inversion of phase, i.e. A 1 > 0 for every v min . Fig. 1 also shows that for the O 1 , O 8 and O 11 interactions, the heavier the target the larger A 1 . This is expected, as the corresponding scattering cross section depends on W τ τ M , and, therefore, scales with the number of nucleons squared in the small momentum transfer limit. In contrast, for the spin-dependent O 7 interaction, aluminium exhibits the largest modulation amplitude. Fig. 2 shows A 1 as a function of v min for the same interactions and targets as Fig. 1, but now for m χ = 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×day.

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 min . In the second one, A 1 > 0 for every v min . In this subsection we investigate the validity of the single cosine approximation, i.e. the hierarchy A 0 A 1 A n≥2 , for the two families   Here, we consider two benchmark scenarios: 1) A DM particle with interactions of 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 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 Al 2 O 3 crystals. The spike around 200 km s −1 in the |A 1 |/A 0 curve for O 1 reflects the change of sign of A 1 for this interaction. Similarly, the smooth |A 1 |/A 0 curve associated with 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 O 1 interaction, and at 415-420 km s −1 and 690-700 km s −1 for the O 7 interaction. Consequently, we conclude that the single cosine approximation is valid at all v min for interactions where A 1 does not change sign, like O 7 . We also find that for interactions like O 1 the single cosine approximation is valid at all v 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 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

Target dependence of the time of maximum rate
We now characterise the interactions in Tab. 1 in terms of the predicted t max (v 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, In terms of the η(v min , t) andη(v min , t) functions, the expected rate of nuclear recoils can be written as    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 max (v 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 min and, in the case of aluminium, for different DM particle masses. In the case of calcium, the curve t max (v min ) is approximately flat, as the term proportional to theη function in Eq. (3.5) dominates over the term proportional to η. This is due to the fact that 16 O and 40 Ca have spin zero, which implies W Σ = 0. As a result, the surviving term in Eq. (2.9) is proportional to W pp M for both the η and theη contributions and the coefficient ofη is proportional to 1/v 2 min , which explains why theη contribution dominates over the η one at small v min .. In contrast, the curve t max (v min ) has a minimum above 100 km s −1 in the case of aluminium. This arises from the interplay of the η andη terms in Eq. (3.5). However, this minimum is pronounced only for m χ > 10 GeV. Therefore, for magnetic dipole DM-nucleus interactions, aluminium, calcium and oxygen are characterised by the same, approximately flat t max (v min ) curve for small m χ values. On the other hand, for the familiar spin-independent interaction, the t max (v min ) curve span a period of about five months. Consequently, it is reasonable to expect that timing information on the observed nu- clear 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 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 dσ T /dE R .

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, ∆S, that we divide into N bins. ∆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×time intervals, when we assume that timing information on the nuclear recoils is also accessible (2D analysis).

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, λ b , in a pre-defined signal region ∆S is an input. We set λ b to the desired value by varying the overall background normalisation, 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 CaWO 4 , whereas here we consider Al 2 O 3 as target material. Specifically, for the experimental background we assume the energy spectrum [38], 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 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 σ SI = 2.67 × 10 −38 cm 2 , r = c p /c n = −0.76 and m χ = 11.17 GeV, respectively [41]; and 2) σ SI = 4 × 10 −42 cm 2 , r = c p /c n = 1 and m χ = 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 χ = 3 GeV [42]. Notice that m χ = 3 GeV is the smallest DM mass that can induce values of v min below 200 km s −1 for recoil energies below 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 ∆S is hence λ s = i s i . Here, we adopt the notation of [43] and introduce the signal strength parameter, µ. For µ = 0 there is no DM signal contribution to the observed nuclear recoil events, while for µ = 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 θ (see Eq. (4.4)) 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], and L (µ, θ) is the likelihood function, Here,θ is the value of θ that maximises L when µ = 0, while (μ,θ) 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 µs i + θb i , where µ and θ are the hypothesised values for µ and θ, respectively. For example, when sampling n i under the background-only hypothesis, µ = 0 and θ = 1. When sampling n i under the signal-plus-background hypothesis, µ = 1 and θ = 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, where q med 0 is the median of f (q 0 |1). In order to test the importance of the nuisance parameter θ, we also evaluate the p-value in Eq. (4.5) for θ set to the constant value θ = θ = 1, so that θ =θ = 1.

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 χ = 3 GeV, r = c p /c n = 1 and σ SI = 4 × 10 −42 cm 2 , 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 Al 2 O 3 detector with an exposure of 23 g×year, when considering the interval [0.012−5.4] keV as energy window. Within the alternative hypothesis, we assume m χ = 3 GeV and couplings set to σ MD ≡ 4α 2 λ 2 χ = 4.72 × 10 −41 cm 2 . 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, and 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.8).
Analogously to the case of DM signal discovery, we sample n i from a Poisson distribution of mean µ s i + θb i . Sampling n i under the null hypothesis, i.e. µ = 1, θ = 1 and s i = s 0 i , 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. µ = 1, θ = 1 and s i = s a i . 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, where q med is the median of the pdf f a (q|1).

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 med (q med 0 ) obeys a binomial distribution of variance σ 2 = pN (1 − p), where 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/ √ 4N . Since we choose p = 0.1 (i.e. 90% C.L.) as a reference p-value, we set N to 10 3 . In the tables presented below, p 0 indicates that the result of our calculations is comparable or smaller than the associated error, 1/ √ 4N 0.016.

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 DMnucleon scattering cross section, the proton to neutron coupling ratio and the DM particle mass to: 1. σ SI = 2.67 × 10 −38 cm 2 , r = c p /c n = −0.76 and m χ = 11.17 GeV, respectively [41]. Furthermore, ∆S = [105 eV, 16 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 Al 2 O 3 as target material and a minimum exposure of 23 g×158 days. In the 1D analyses nuclear recoils are distributed in 10 energy bins, while in the 2D analyses they are distributed in 10 2 energy×time bins. We test the results using both linear and logarithmic binning size. We use the test statistic in Eq. (4.2).

Benchmark 1
Focusing on the first benchmark scenario, we now present our results on the ability of a lowthreshold experiment to reject the background-only hypothesis. For this scenario, Fig. 6 shows the predicted energy (left panel) and energy×time (right panel) distributions of DM-induced nuclear recoil events. As expected from our analysis of A 1 (v min ) and t max (v 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 θ is set to the constant value θ = θ = 1. We will treat θ as a nuisance parameter later on. As an illustrative example of the geometric interpretation of the p-value for signal discovery in Eq. (4.9), 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, λ b = 10 3 , λ s 141 and ∆S = [105 eV, 16 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.
as signal region. As reported in the legend, the p-value associated with this specific example is p = 0.046. Tab. 2 extends the results in Fig. 7  constant value θ = θ = 1 for different combinations of λ b , the expected number of background events, and experimental exposures. More specifically, we consider values of λ b ranging from 10 3 to 10 5 and experimental exposures varying from 23 g×year to 3×23 g×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 θ is considered as a nuisance parameter. Tab. 3 shows the p-values for signal discovery that we find in simulations where λ b varies from 500 to 2×10 5 , the experimental exposure varies from 23 g×158 days (corresponding to the exposure reached in the CRESST-III analysis [17]) to 23 g×2 years and θ 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 θ 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.2) for θ = θ = 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.

Benchmark 2
Now we turn our attention to the second benchmark scenario considered in this investigation, where σ SI = 4 × 10 −42 cm 2 , r = c p /c n = 1 and m χ = 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 θ 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, λ b , vary from 500 to 5 × 10 5 and consider experimental exposures varying from 23 g×158 days to 1 kg×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 λ b = 500 and an exposure of 23 g×158 days, to about 0.06, for λ b = 5 × 10 5 and an exposure of 1 kg×5 years. Consistently with benchmark 2 being about two orders of magnitude below the current CRESST-III exclusion limit on σ SI for m χ = 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×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.

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 χ = 3 GeV, r p /r n = 1 and σ SI = 4×10 −42 cm 2 . This set of parameters predicts 2.2 counts in the [0.012−5.4] keV energy window for an Al 2 O 3 detector with an exposure of 23 g×year. We model the alternative hypothesis by assuming m χ = 3 GeV and a cross section σ MD ≡ 4αλ 2 χ = 4.72 × 10 −41 cm 2 . Consistently, this reference cross section value produces 2.2 counts in the same Al 2 O 3 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×158 days to 1 kg×5 years, and λ b varying from 10 3 to 2 × 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 θ 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×year and 2 × 10 3 background events.  Table 5. p-values for rejecting spin-independent interactions with m χ = 3 GeV and σ SI = 4 × 10 −42 cm 2 in favour of a magnetic dipole coupling producing the same number of events in ∆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.

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 Al 2 O 3 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 max (v 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×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 spinindependent 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×2 years and a background level of 5 × 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%.