Scintillating bolometers: a key for determining WIMP parameters

In the last decade direct detection Dark Matter (DM) experiments have increased enormously their sensitivity and ton-scale setups have been proposed, especially using germanium and xenon targets with double readout and background discrimination capabilities. In light of this situation, we study the prospects for determining the parameters of Weakly Interacting Massive Particle (WIMP) DM (mass, spin-dependent (SD) and spin-independent (SI) cross section off nucleons) by combining the results of such experiments in the case of a hypothetical detection. In general, the degeneracy between the SD and SI components of the scattering cross section can only be removed using targets with different sensitivities to these components. Scintillating bolometers, with particle discrimination capability, very good energy resolution and threshold and a wide choice of target materials, are an excellent tool for a multitarget complementary DM search. We investigate how the simultaneous use of scintillating targets with different SD-SI sensitivities and/or light isotopes (as the case of CaF2 and NaI) significantly improves the determination of the WIMP parameters. In order to make the analysis more realistic we include the effect of uncertainties in the halo model and in the spin-dependent nuclear structure functions, as well as the effect of a thermal quenching different from 1.


Introduction
Weakly Interacting Massive Particles (WIMPs) can be directly detected through their scattering off target nuclei of a detector. 1 In the last decades, numerous experiments, using different targets and detection techniques, have been searching for WIMPs or are currently taking data. Some of them have searched for distinctive signals, such as an annual modulation in the detection rate: DAMA 2 and DAMA/LIBRA, 3,4 using NaI scintillators, have reported a highly significant signal (9.3σ) and CoGeNT 5, 6 claimed a less significant evidence (2.2σ) in the first three years of its data, gathered with a Ge semiconductor. Moreover, CoGeNT, 7 CRESST 8 (using CaWO 4 scintillating bolometers) and CDMS II (with data from its Si detectors) 9 have reported excesses of events at low energies that could be compatible with a signal produced by light WIMPs with a mass of the order of 10 GeV. On the other hand, XENON10, 10 XENON100, 11 LUX 12 (also based on Xe), the abovementioned CDMS II, 13 EDELWEISS 14,15 (with Ge), KIMS 16 (with CsI), PICASSO 17 (with C 4 F 10 ), SIMPLE 18 (with C 2 ClF 5 ) and COUPP 19 (with CF 3 I) have obtained negative results setting more stringent upper bounds on the WIMP-nucleon cross sections. Currently the strongest limits are obtained by the LUX collaboration, excluding spin-independent WIMP-nucleon elastic scattering cross sections larger than 7.6×10 −46 cm 2 for a WIMP mass of 33 GeV, and the SuperCDMS collaboration for low mass WIMPs. 20,21 In the next years new experiments and upgraded versions of the existing ones are going to explore even smaller cross sections, closing in on DM searches.
The final goal of all these experiments is to determine the nature of DM, measuring some of its properties (namely its mass and interaction cross section with ordinary matter). Signals from different targets are needed, since they can provide complementary information which can lead to a better determination of the DM parameters. 22,23 In a previous paper 24 we analysed the complementarity of a Ge and a Xe experiment with energy thresholds and resolutions already achieved by CDMS and XENON100 experiments, respectively, and with background levels expected for their corresponding extensions (SuperCDMS 25 and XENON1T 26 ). For different WIMP scenarios, we assumed hypothetical detections with an exposure of 300 kg×yr in both experiments and we concluded that the combination of data from Xe and Ge-based detectors might not lead to a good reconstruction of all the WIMP parameters, since there is a degeneracy in the SI and SD parts of the scattering WIMP-nucleus cross section, and both targets have very similar SI over SD sensitivity (see also Ref. 27 for a recent study on the non-complementarity of Xe and Ar). We showed that incorporating targets with different sensitivities to SI and SD interactions could significantly improve the reconstruction. We considered the case of some of the most promising scintillating bolometric targets: CaWO 4 (currently used by CRESST), Al 2 O 3 and LiF (studied by ROSEBUD, 28 that could be considered in the future as additional targets in EURECA, 29 a European collaboration that plans to search for WIMPs with a 1-ton cryogenic hybrid detector).
We observed that the inclusion of CaWO 4 (being mainly sensitive to SI couplings) only leads to a total complementary result for a WIMP of 50 GeV in a small region of the plane (σ SI , σ SD ) in which the expected events in Ge and Xe are mainly due to SD interactions. On the other hand, Al 2 O 3 and LiF (being more sensitive to SD interactions) achieve complementarity with germanium and xenon in regions of the parameter space where the rate in the latter is dominated by SI couplings. We also determined the exposures and background levels required by the bolometers to be complementary to Ge-and Xe-based experiments.
In this paper we follow the same strategy and reanalyze the role of Ge-and Xe-based experiments in light of the improved (or potential) energy thresholds in CDMS and LUX 3 . We also study the complementarity with two additional targets: CaF 2 and NaI. The first one has already been used as scintillating bolometer, 31,32 whereas the construction of a bolometer based on NaI (which is a hygroscopic and fragile material) is an ongoing R&D project of the Zaragoza group. 33 We include in our analysis not only the effect of the previously considered uncertainties in the halo parameters and SD structure functions, but also the possible influence of the thermal quenching between nuclear and electron recoils in the complementarity of these targets.
The structure of this article is as follows: Sec. 2 is a short summary of the methodology we follow in reconstructing the WIMP parameters from the (simulated) data in direct detection experiments. In Sec. 2.1 we address the most relevant uncertainties in the analysis, in particular the astrophysical ones (due to our imperfect knowledge of the DM halo of the Milky Way), those related to the SD Structure Functions (SDSF) parametrizing the spin content of the nucleons in the target and, finally, the effect of changing the thermal quenching q. In Sec. 2.2 we present the results for some selected benchmarks when considering only Ge and Xe experiments, finding that the combination of data from these two targets contributes to a better measurement of the WIMP parameters, but a degeneracy in the SD and SI independent cross section usually remains. In Sec. 3 we describe the characteristics of the scintillating targets under study (i.e. CaF 2 and NaI ). In Sec. 4 we show how their inclusion can lead to a better determination of the DM mass and scattering cross section, breaking in some cases the SI-SD degeneracy. Finally, conclusions are presented in Sec. 5.

Reconstructing WIMP parameters from signals in direct detection experiments
In the standard analysis framework for WIMP direct detection 34,35 (see also Refs. 36 and 37 for recent reviews) the WIMP-nucleus scattering cross section is separated into a SI and a SD contribution, with f p and a p being the corresponding effective couplings to protons and f n and a n to neutrons. In order to reduce the number of parameters that characterize the expected event rate, we assume here that the SI coupling is isospin-invariant (f p = f n ) and we take a specific relation between a p and a n (namely a p /a n = −1). Under these assumptions, the generic WIMP is completely determined by its mass m χ , the SI contribution to the WIMP-nucleon cross section σ SI and the SD component σ SD . Thus, the total number of WIMP recoil events in a given energy window can be expressed as where, for each isotope, f is its mass fraction in the detector, S p and S n are the expectation values of the total spin operators for protons and neutrons respectively and the coefficients C SI and C SD can be written as follows, (2) ǫ is the experimental exposure, ρ 0 is the local WIMP density, f (v) is the WIMP speed distribution in the Earth reference frame normalized to unity, µ n is the WIMP-nucleon reduced mass, E R is the nucleus recoil energy, F 2 SI (F 2 SD ) is the SI (SD) nuclear form factor, A is the nucleus mass number, and J its nuclear spin.
We focus the analysis on two benchmark cases. For each of them (and for each target independently) we calculate the signal that such WIMPs would produce in that specific detector by computing the number of recoil events, {λ i } a , expected for target a in the i-th bin of N evenly-spaced energy bins contained in the energy window for WIMP search of each experiment. These expected events, {λ i } a , represent our experimental data, D, and we are interested in estimating how well such simulated measurements can be used to reconstruct the WIMP parameters.
In order to do so, we perform scans over the parameter space (m χ = 1 − 10 5 GeV, σ SI = 10 −12 −10 −6 pb, and σ SD = 10 −8 −1 pb). For every point in the scan we compute the number of recoil events N a i in the i-th energy bin for every target a and then compute the likelihood comparing N a i with the prediction of the benchmark model in the same energy bin for the same target, assuming that data from each experiment follow independent Poissonian distributions. We present the results as 68% and 99% confidence regions in the profile likelihood (PL). The nuclear and astrophysical uncertainties are considered as nuisance parameters. The scans are performed with MultiNest 3.0 38 interfaced with our own code for the computation of the number of recoil events and of the likelihood. Logarithmic flat priors are assumed for the three variables. We refer to Ref. 24 for a detailed description of how the scans are performed.
If only one target is considered, the reconstruction of WIMP parameters is affected by degeneracies, since the number of events detected can be explained by different combinations of (m χ , σ SI , σ SD ). Such degeneracy can be broken by including more targets in the analysis: in Ref. 24 we defined "complementarity" as the situation in which a certain set of experiments manages to determine m χ , σ SI and σ SD with a certain finite accuracy, or, equivalently, when 68% confidence level of the 2-dimensional contours are closed simultaneously in the three planes (m χ , σ SI ), (m χ , σ SD ) and (σ SI , σ SD ).
The following two WIMP benchmarks will be considered in the remaining sections: • VL-SI: m χ =20 GeV, σ SI =10 −9 pb, σ SD =10 −5 pb, corresponding to a very light WIMP for which the SI contribution dominates the detection rate in Ge and Xe, • L-SD: m χ =50 GeV, σ SI =10 −10 pb, σ SD =1.5×10 −4 pb, a light WIMP for which the SD contribution dominates in Ge and Xe.

Including uncertainties
The expected DM signal depends on parameters affected by large uncertainties. In the following, we will take into account uncertainties in the velocity distribution of DM in the Milky Way halo, the spin dependent form functions for the target nuclei and the performance of the detector. We considered a velocity distribution function that differs from the standard halo model by the presence of a high-velocity tail. Such a model, adopted from Ref. 39, is well motivated by N-body simulations and the velocity distribution can be written as follows, where N k = v 3 0 e −y 2 e ye 0 y 2 (e −(y 2 −y 2 e )/k − 1) k dy, y e = v esc /v 0 and k is the parameter that quantifies the deviation from the standard halo model, recovering it in the limiting case of k=0. 39 This expression for the velocity distribution depends on three parameters: v esc , v 0 and k. In order to account for our ignorance on the true velocity distribution of the DM in the halo of our Galaxy we leave such parameters free to vary within the following ranges: v esc = [478, 610] km/s, v 0 = [170, 290] km/s, and k = [0.5, 3.5]. We also scan over the local WIMP density ρ 0 , in the range between 0.2 and 0.6 GeV/cm 3 . All these parameters are subject to a uniformly flat prior distribution.
Regarding the WIMP interaction with the nucleus, it has been shown in Ref. 40 that uncertainties in the SI and SD form factors of the target nuclei play a very different role. In the case of SI interactions, differences in the form factor can be safely neglected (in the present paper we have used the Helm form factor). On the other hand, for SD interactions the expressions of the form factors are more dependent on the nuclear model. These differences can significantly affect the expected WIMP rate and, thus, the reconstruction of the WIMP parameters (especially when SD interactions play a relevant role).
To take into account such uncertainties, the SDSFs are parametrized as follows, 40 where u is an adimensional quantity proportional to the square of the momentum transfer, u = (qb) 2 /2, in terms of the oscillator size parameter b = A 1/6 . Note that for the case a p /a n = −1 the only contribution comes from the S 11 SDSF. Tab. 1 shows the ranges in which the three parameters N , α and β have been allowed to vary for each nucleus in order to reproduce the various determinations of the form factors available in the literature (see Ref. 40 for a detailed explanation). Results for the SDSFs of the isotopes relevant in this work in the case a p /a n = −1 are displayed in Fig. 1 (light blue area) together with the most relevant nuclear calculations.
Finally, important systematics can also arise from the detection technique itself. Among these, we consider the effect of the thermal quenching factor, q, that measures the relative efficiency in the conversion into measurable thermal signal of the nuclear recoils energy deposition with respect to that corresponding to electron recoils, since the detectors are calibrated with gamma sources and the measured spectra are given in electron-equivalent energy. This factor is typically assumed equal to one for bolometers but small deviations (of about 10-15%) have been measured in different detectors (see for example Ref. 47

Results for Ge and Xe
As in our previous paper, 24 we start the analysis studying the complementarity of two experiments, based, respectively on Ge and Xe. Such elements are employed by the collaborations producing the currently most stringent limits on WIMP properties and are contemplated in projects planning to extend the search to the ton scale (e.g. EURECA, SuperCDMS and XENON1T) or even to the multi-ton scale (LZ 48 and DARWIN 49 ). Consequently, these targets are expected to represent the most sensitive experiments (at least in the most general WIMP scenarios) in the near future. For our study, we have assumed a positive result (WIMP detection) in two experiments, one using a Ge-based target and the other using Xe. We consider the two detections combined when reconstructing the WIMP parameters. The same exposure (ǫ = 300 kg×yr) is assumed for both experiments, as well as zero background 4 . The energy window is set to  keV for Ge and  keV for Xe, where the lower values account for the recent or potential improvements in nuclear recoil energy thresholds of some Xe and Ge experiments. Tab. 2 shows the expected number of WIMP recoil events for the considered benchmarks over the whole energy range. Fig. 2 shows the 68% and 99% confidence level contours for the three WIMP parameters projected onto the corresponding two-dimensional plots (σ SI , m χ ), (σ SD , m χ ), and (σ SD , σ SI ) for the benchmark VL-SI. The yellow dot represents the nominal value and the circled cross is the best-fit point. As we showed in Ref. 24, the combination of data from Ge and Xe leads to a substantial reduction in the contours of the reconstructed WIMP parameters. The improved energy threshold also contributes to this. In particular, for this benchmark the mass of the WIMP can be well determined (the contours using only one target would not be closed). However, there remains a degeneracy in both cross sections, σ SI and σ SD , for which only upper limits are derived. This is due to the similar sensitivity to SI/SD interactions of Ge and Xe in this point.
Analogously, Fig. 3 displays the contour plots for the benchmark L-SD. The reconstruction of WIMP parameters is similar, although in this case σ SD is better bounded and even a lower limit is derived at 68% C.L. Nevertheless, at 99% CL the degeneracy between σ SI and σ SD still remains and only an upper bound is obtained for σ SI which is far from the nominal value.
Thus, although the combined data from Ge and Xe experiments can be used to significantly improve the determination of WIMP parameters, the degeneracy in the SI and SD components of the scattering cross section might be difficult to break. As we will argue in the following sections, incorporating data from a third target with a different sensitivity to these components can help solving this problem.

Scintillating bolometers for dark matter searches
Compared to other hybrid detectors with discrimination capability, scintillating bolometers, originally proposed in 1988, 50, 51 have the advantage of a wide target choice. This makes it possible to select intrinsically radiopure materials and combine different nuclei to maximize the explored region of the WIMP parameter space (high mass number A for large SI coupling, low A to enhance sensitivity to light WIMPs, or non-zero nuclear spin for sensitivity to SD interaction, to name just a few possibilities).
The energy threshold that has been achieved in the heat signal with cryogenic detectors is as low as ∼1 keV. However, when looking for nuclear recoils, the discrimination threshold is determined by the target light yield and the sensitivity of the optical detector. Usually this role is played by a second low-mass large-area bolometer facing the primary one. Optimizing the sensitivity and response of the optical bolometer is a very active ongoing research field (see for example Refs. 52 and 53) and lower thresholds are expected in a near future. Nevertheless, in this paper, we follow the same approach of Ref. 24 and take a reference energy threshold of 10 keV (a value already observed in some scintillating bolometers 8,54 ) for the bolometric targets under study.
In our previous work 24 we studied the complementarity of three scintillating materials: CaWO 4 , Al 2 O 3 and LiF. CaWO 4 is the current target of the CRESST experiment 8 and was used also by the ROSEBUD collaboration in the first underground DM search with light and heat discrimination. 55 It constitutes the baseline for the EURECA scintillating targets. Al 2 O 3 , used by ROSEBUD 56,57 and by CRESST in the first phase of the experiment, 58 is particularly interesting for its sensitivity to low mass WIMPs. Finally, LiF, also sensitive to light WIMPs and SD interaction, has been used by the ROSEBUD collaboration for DM searches and as neutron detector, showing that its use in a DM experiment could allow for thermal neutron monitoring. 59,60 However, the light yields achieved so far do not provide a good discrimination threshold, so further developments are needed in order to use this target in a DM experiment.
In this paper we focus on two other targets: CaF 2 and NaI. Fluorine-based scintillators are particularly attractive for DM searches because of the sensitivity of 19 F (J=1/2, 100% isotopic abundance) to SD interactions. Among them, CaF 2 presents the highest light yield 61 and has already been used in several DM searches as scintillator at room temperature. [62][63][64] It was the target material of the first scintillating bolometer ever constructed, 31 although in that experiment the light measurement was performed with a silicon photodiode, less sensitive than the semiconductor bolometers usually used in recent setups. 32 Scintillation at low temperature has been studied for pure and europium-activated targets, resulting in good scintillation at 1 K specially for doped samples, 65,66 although the radiopurity levels achieved in this case are usually worse. NaI, on the other hand, is one of the most widely used scintillators for γ spectroscopy due to its very high light yield. As mentioned above, this is the target used by DAMA/LIBRA and other proposed DM experiments looking for annual modulation. 67,68 Although NaI is usually doped with Tl for room temperature applications, the pure material is known to scintillate better at temperatures of a few Kelvin 69 (nevertheless, an increase in light yield of the Tl-doped material below 30 K has been recently reported 33,70 ). Despite its high light yield at low temperature and intrinsic interest for DM searches, this material has not been tested yet as a bolometer due to its fragility and high hygroscopicity.

Results with bolometric targets
Let us now investigate the complementarity potential of scintillating bolometer targets of CaF 2 and NaI with the Ge and Xe experiments. For both bolometric targets, we assume an energy window from 10 to 100 keV, a 5% energy resolution and, as we have done previously for Ge and Xe, a total exposure of 300 kg×yr and a zero background experiment. Tab. 3 gives the number of recoil events for each of the bolometric targets. In the case of NaI, three different quenching factors have been considered (q = 0.85, 1 and 1.15). Following the same procedure of Ref. 24, for each benchmark and target we have derived the contour plots after the combination of data from a Ge detector, a Xe detector and the corresponding bolometric target (see Figs. 4 to 7). Results are shown as blue contours, while black lines correspond to the case when only Ge and Xe are used.   In benchmark VL-SI, for which Ge and Xe exhibited a degeneracy in the (σ SI , σ SD ) plane, CaF 2 provides a good complementarity, allowing the full reconstructions of the WIMP parameter space (see Fig. 4). This is because for this BM the 97% (95%) of the signal in Ge (Xe) is due to the SI component whereas in the case of CaF 2 , a target very sensitive to SD WIMP-nucleon interactions, the 80% of the total rate is due to the SD component.
The results using NaI are represented in Fig. 5, where the three rows correspond to the three different values considered for the quenching factor. As we can observe, we are able to obtain closed contours for σ SI , but not for σ SD (see Fig. 5). The reason is that, as in the case of Ge and Xe, the signal for this target is dominated by the SI contribution (approximately 85% of the total rate). The change in the quenching factor (which can be understood as a shift in the energy window of nuclear recoils) leads to variations in the number of events due to SD and SI interactions. More importantly these do not change by the same amount, since the energy dependence of the SD and SI form factors is different. For NaI we observe that the relative contribution due to the SD term increases as the quenching factor decreases, shifting from 14% at q=1.15 to 17% at q=0.85. This implies that, for this benchmark, the complementarity with Ge and Xe is better for q = 0.85, as we can observe in Fig. 5. The effect is clearer in the 1-D profile likelihood of the SD cross section shown in Fig. 8. Notice also that, although the upper limit on σ SD is more stringent for q=0.85, the derived 1-D profile likelihood is practically flat (Fig. 8 right) and that leads to a failure in the estimation of σ SD by the best-fit point (Fig. 5).
The results for benchmark L-SD are shown in Fig. 6 for the combination of data from Xe, Ge, and CaF 2 , and in Fig. 7 for NaI. As we can see in Tables 2 and 3, in this benchmark the WIMP interactions are dominated by the SD contribution for all the targets. Consequently, the degeneracy in the (σ SI , σ SD ) plane is not completely removed, although the contours are substantially reduced with respect to the case with Ge and Xe alone. In particular, closed contours appear for σ SD around the nominal value with both CaF 2 and NaI, but only an upper bound for σ SI is obtained. In this benchmark the effect of the quenching factor is quite imperceptible (see also Fig. 9) because the relative contribution of the SD term is practically the same (approximately 7.5%) for the three values of q.

Conclusions
Following the work done in Ref. 24, where we investigated the determination of WIMP parameters (m χ , σ SI , σ SD ) from a hypothetical direct DM detection with multiple targets, in this paper we have extended the analysis to consider the effect of lower thresholds in Ge and Xe targets, as well as the complementarity potential of two new bolometric targets: CaF 2 and NaI.
We first considered the combination of data from Ge and Xe targets, for both of which we assumed a low energy threshold of 3 keV to account for recent or projected experimental improvements. We studied two benchmark scenarios, featuring a very light WIMP (m χ =20 GeV, σ SI =10 −9 pb, σ SD =10 −5 pb) in which SI contribution dominates the detection rate in both Ge and Xe, and a light WIMP (m χ =50 GeV, σ SI =10 −10 pb, σ SD =1.5×10 −4 pb) in which the SD contribution dominates. Although the combination of data from both targets allows a significant improvement in the reconstruction of DM parameters, a degeneracy in the (σ SI , σ SD ) plane usually remains in the points in the parameter space where both targets have similar SI/SD ratios.
Scintillating bolometers, with very good energy threshold and resolution and particle discrimination capability, provide a wide choice of absorber materials that allows to select interesting targets form the point of view of its complementarity with other experiments. In Ref. 24 we studied how certain bolometric targets (CaWO 4 , Al 2 O 3 and LiF) could provide complementary information to data from Ge or Xe based experiments. In this work we have extended the analysis to other two scintillating targets (CaF 2 and NaI), and considered also the effect of an uncertainty in the thermal quenching factor of ±15%. Both targets are sensitive to the SD component of the WIMP-nucleus interaction (particularly CaF 2 thanks to the presence of 19 F). We have shown how the inclusion of one of these targets together with Ge and Xe can help breaking the degeneracy in the (σ SI , σ SD ) plane. In particular, in the points of the parameter space for which the rate in Ge and Xe is dominated by the SI contribution and the rate in CaF 2 is mostly SD, the three DM parameters can be reconstructed. In other examples, although the degeneracy cannot completely removed, at least one of the components of the WIMP-nucleus scattering cross section can be determined.
We have also shown how a small uncertainty in the thermal quenching factor can modify noticeably the parameter reconstruction.