Low Energy Analysis Techniques for CUORE

CUORE is a tonne-scale cryogenic detector operating at the Laboratori Nazionali del Gran Sasso (LNGS) that uses tellurium dioxide bolometers to search for neutrinoless double-beta decay of $^{130}$Te. CUORE is also suitable to search for low energy rare events such as solar axions or WIMP scattering, thanks to its ultra-low background and large target mass. However, to conduct such sensitive searches requires improving the energy threshold to 10 keV. In this paper, we describe the analysis techniques developed for the low energy analysis of CUORE-like detectors, using the data acquired from November 2013 to March 2015 by CUORE-0, a single-tower prototype designed to validate the assembly procedure and new cleaning techniques of CUORE. We explain the energy threshold optimization, continuous monitoring of the trigger efficiency, data and event selection, and energy calibration at low energies in detail. We also present the low energy background spectrum of CUORE-0 below 60keV. Finally, we report the sensitivity of CUORE to WIMP annual modulation using the CUORE-0 energy threshold and background, as well as an estimate of the uncertainty on the nuclear quenching factor from nuclear recoils in CUORE-0.

Abstract CUORE is a tonne-scale cryogenic detector operating at the Laboratori Nazionali del Gran Sasso (LNGS) that uses tellurium dioxide bolometers to search for neutrinoless double-beta decay of 130 Te. CUORE is also suitable to search for low energy rare events such as solar axions or WIMP scattering, thanks to its ultra-low background and large target mass. However, to conduct such sensitive searches requires improving the energy threshold to 10 keV. In this paper, we describe the analysis techniques developed for the low energy analysis of CUORE-like detectors, using the data acquired from November 2013 to March 2015 by CUORE-0, a single-tower prototype designed to validate the assembly procedure and new cleaning techniques of CUORE. We explain the energy threshold optimization, continuous monitoring of the trigger efficiency, data and event selection, and energy calibration at low energies in detail. We also present the low energy background spectrum of CUORE-0 below 60 keV. Finally, we report the sensitivity of CUORE to WIMP annual modulation using the CUORE-0 energy threshold and background, as well as an estimate of the uncertainty on the nuclear quenching factor from nuclear recoils in CUORE-0.  [1,2]. In 0νβ β decay, two neutrons in an atomic nucleus simultaneously decay to two protons and two electrons, without emitting any electron antineutrinos. The experimental signature of 0νβ β decay is a sharp peak at the tail end of the two-neutrino double-beta decay (2νβ β ) summed energy spectrum. The bolometric technique of CUORE offers an excellent energy resolution of ∼ 5 keV FWHM at the Q-value of 130 Te, 2527.5 keV [3][4][5], which suppresses the 2νβ β decay background leaking into the 0νβ β decay signal region of interest (ROI) [6].
While CUORE will be one of the leading 0νβ β decay experiments during its scheduled 5 years of data-taking, it will also benefit from the ultra-low background and large target mass to search for lower energy rare events, such as the direct detection of Weakly Interacting Massive Particle (WIMP) dark matter or solar axions [9]. WIMP direct detection is possible with terrestrial detectors by measuring nuclear recoils produced as WIMPs scatter off nuclei in the target material [10]. The resulting energy spectrum falls quasi-exponentially as a function of energy and extends to only a few tens of keV for typical WIMPs with masses of O(100 GeV/c 2 ). For WIMPs in the galactic halo, an annual modulation of the event rate is expected due to the Earth's motion relative to the dark matter halo of the Milky Way [11,12], with event rates highest in June, when the Earth's relative velocity with respect to the halo is maximal, and lowest in December. Alternatively, solar axions can be detected by the inverse Primakoff effect in the Coulomb field of the crystal, with a signal from the M1 transition of 57 Fe expected at 14.4 keV [13]. A critical requirement for both these searches is the achievement of an energy threshold of < 10 keV and sufficient rejection of low energy noise and/or spurious events, in addition to a detailed understanding of the low energy backgrounds and adequate detector stability.
Since Cuoricino, which proved that the CUORE detector technology is well suited for searching for 0νβ β decay [7], the CUORE collaboration has worked to lower the energy thresholds to perform low energy rare event searches. We developed a new low energy software trigger, the "optimal trigger" (OT) [14], based on filtering the continuous data stream before the application of the trigger condition. We implemented this algorithm in test measurements of the bolometric performance of a small number of CUORE crystals in a dedicated setup (CUORE Crystal Validation Runs, CCVR) [15]. The results were encouraging; we were able to identify events with energies as low as 3 keV with a trigger efficiency of > 90% in three bolometers out of four [16]. In CUORE-0, we improved the OT technique and developed new software and hardware tools for the low energy analysis. These tools include the continuous monitoring of the trigger efficiency, the development of low energy event selection criteria, and a low energy calibration.
In this paper, we describe our low energy analysis techniques in detail, present the energy thresholds and spectrum from the CUORE-0 experiment and report the sensitivity of CUORE to WIMP-induced annual modulation assuming the same energy thresholds and background. We also evaluate the uncertainty on the nuclear quenching factor of TeO 2 obtained from CUORE-0 data. Specifically, Section 2 describes the experimental setup and data production of CUORE-0. Section 3 explains the OT algorithm and trigger efficiency evaluation, and section 4 details the data and event selection criteria. Section 5 presents the low energy spectrum as well as the determination of the analysis threshold and the evaluation of the low energy calibration uncertainty. Section 6 outlines the analysis developed for a WIMP search, including the nuclear quenching factor of TeO 2 estimation, as well as a study of the sensitivity of CUORE to WIMP-induced annual modulation. Finally, we present the summary in Sec. 7.
2 The CUORE-0 experiment CUORE-0 comprised 52 TeO 2 crystals with a total active mass of 38.4 kg. The crystals were arranged in a single tower, with 13 planes of four 5 × 5 × 5 cm 3 crystals held securely inside a copper frame by polytetrafluoroethylene supports. The detector was hosted in the same cryostat used for Cuoricino at a base temperature of ∼ 10 mK, and used the same shielding and electronics. The detector design, construction, and operation are detailed in [17]. Each crystal was instrumented with a neutron-transmutation-doped (NTD) germanium thermistor [18] to read the thermal signal, and a silicon resistor [19], used as a Joule heater to inject reference pulses of constant energy every 300 s. The reference pulses were mainly used to correct the thermal gain against the drift in temperature, but they also played a fundamental role in determining the OT efficiency, as explained in Sec. 3. The thermal readout of each bolometer thermistor was in the form of a voltage waveform continuously acquired with a sampling frequency of 125 S/s. Taking advantage of this relatively low acquisition rate, we could record the continuous data stream without hardware trigger. This allowed us to reprocess the raw data with different software trigger algorithms and to optimize the energy thresholds offline.
We collected data in one-day-long runs. Approximately once per month, we calibrated the detector by irradiating it for about 3 days using thoriated tungsten wires inserted between the outer vacuum chamber of the cryostat and the external lead shielding. The basic analysis unit is the dataset, which is composed of initial calibration runs, approximately 3 weeks of physics runs, and final calibration runs. For the low energy analysis, we also performed a dedicated measurement before each final calibration using low energy pulses generated by the Joule heater, with energies ranging from 0 to 50 keV.
A comprehensive description of the standard CUORE-0 data processing procedure for 0νβ β decay and 2νβ β decay can be found in [20]. The following summarizes the major steps of the data processing that are common to both low energy and high energy (0νβ β decay and 2νβ β decay) analyses. After application of the software trigger, we store events in 5 s windows and evaluate the pulse amplitude using the optimal filter (OF) [21]. The OF weights each frequency component by the expected signal-to-noise ratio, calculated as the ratio between the average pulse (obtained from the 2615 keV γ rays in the calibration data), and the average noise power spectra (NPS). To calculate the NPS we average baselines recorded in windows without a signal, that are acquired simultaneously on all bolometers every 200 s and selected for the NPS after some quality checks (basically we require that no pulses or pulse tails are present in the window). The drift in signal gain due to temperature fluctuations in the bolometer is corrected by performing a linear regression between the detector baseline voltage, a proxy for the detector temperature, and the amplitude of the reference pulser events. Finally, the voltage readout is converted to energy using the numerous γ ray peaks between 511 and 2615 keV from the daughter nuclei of 232 Th in the calibra-tion data. The mapping from pulse amplitude to energy is described as a second-order polynomial with zero intercept to take into account possible nonlinearities, as those originated from the pulse shape dependence on energy.

Optimal trigger optimization and efficiency evaluation
The energy ROI of the standard data processing for 0νβ β decay analysis is in the MeV range, and we use a simple trigger algorithm which flags an event when the slope of the waveform exceeds a given threshold for a certain amount of time [17]. This results in energy thresholds ranging from 30 to 120 keV, depending on the bolometer, while the use of the OT algorithm is critical for lowering the threshold below 30 keV.
The OT algorithm works as follows. The data buffer is divided into slices that are continuously filtered in the frequency domain with the OF described in Sec. 2. The filtered waveforms are less noisy than the original waveforms, and baseline fluctuations are reduced. This allows us to trigger on the filtered trace in the time domain with a simple threshold as low as < 10 keV. Furthermore, the filter is sensitive to the shape of the expected signal, suppressing trigger on spurious noise-induced pulses.
The algorithm used in CUORE-0 is improved relative to that described in [14] in order to achieve a higher trigger efficiency. In particular, we have removed the veto around high energy pulses that prevented the algorithm from retriggering the symmetric side lobes generated by the OF, which was identified in [14] as the main source of trigger inefficiency. The new algorithm recognizes the side lobes of a high energy pulse as OF artifacts and does not flag them.
We set the trigger threshold independently for every bolometer in every dataset (hereafter "BoDs") based on the noise level. First, we calculate an OT trigger level at θ = 3σ OF , where σ OF is the baseline resolution after applying the OF. The energy-dependent trigger efficiency ε(E) is modeled by the Gaussian cumulative density function which is 50% for E = θ . At this energy we reject 99.86% of baseline noise. The trigger threshold, E 99% trig , is set to the value at which 99% efficiency is reached. The validity of the efficiency calculation is checked at the end of each dataset with dedicated measurements injecting low energy pulses [22]. The bottom plot of Fig. 1 shows an example of trigger efficiency (dashed line) as a function of energy, along with data obtained from corresponding low energy pulser measurements (circles). Vertical gray lines indicate the energies where OT trigger efficiency reaches 50% (θ ) and 99%  (E 99% trig ). In this case, E 99% trig threshold is set at 7.5 keV.
The top plot shows the difference between the modeled trigger efficiency and the data. The CUORE-0 E 99% trig thresholds range from 4 to 12 keV for most BoDs, slightly above those of the CCVR measurement due to a larger noise contribution, as explained in Section 4.3. In Table 1 we summarize the different energy thresholds considered in this work together with the range of values obtained for the CUORE-0 BoDs.

Data selection and energy threshold determination
In this section we detail three steps to select legitimate low energy events. First, we only choose data whose low energy response and stability are verified. Second, we identify and remove non-legitimate events that pass the trigger requirement, such as electronic noise, tower vibrations, pile-ups or particle interactions in a thermistor, where the last appears as a narrow pulse with fast decay time. Last, we remove events that occur simultaneously in more than one bolometer since the probability that WIMP or solar axion interactions occur in more than one bolometer within the coincidence window is essentially zero. In this way we reject muons passing through the tower and radioactive decays that deposit energy in several bolometers (α-decays in the surfaces of the crystals, Compton scattering, cascade γ-rays...).

Dataset-bolometer selection criteria
While CUORE-0 ran from March 2013 to March 2015, we only use 11 datasets from the second data-taking campaign, lasting from November 2013 to March 2015, because of its more stable cryogenic conditions. We exclude some runs (a total of ∼ 5 data-taking days) with an abnormally higher (∼ 10 times) low energy event rate, which we attribute to cryostat instability following a helium refill. To preserve the data quality, we reject the time intervals for each bolometer that exhibit degraded bolometric performance due to large baseline excursions or elevated noise levels, as described in [20]. The total exposure after these preliminary quality checks is 23.15 kg · yr of TeO 2 .
To ensure a stable energy calibration and sufficient resolution at low energies, we additionally require that the residual gain variation of the energy pulser of every BoDs after temperature stabilization do not vary more than ±(σ DS + 1 keV) from the mean over the entire data-taking period, where σ DS is the uncertainty in the pulser position associated with the dataset and the pulser energy ranges from 13 to 27 keV. We also discard BoDs with fewer than 11 events in the region where we evaluate the pulse shape parameter event selection efficiency (35 − 50 keV, see Sec. 4.2). Finally, we exclude run-bolometer pairs with baseline RMS values that are greater than 2σ above the median, where median and σ are calculated for each bolometer over all the datasets. After all the data selection, we use 490 BoDs out of 539 with a total TeO 2 exposure of 20.02 kg · yr.

Event selection criteria
For the standard 0νβ β decay analysis we use a set of six pulse shape parameters to select physical events in TeO 2 based on the pulse shape characteristics [20]. These parameters, however, lose rejection power at low energy due to the worse signal-to-noise ratio. Instead, the OT algorithm provides us with a powerful shape parameter, OT χ 2 , to select legitimate signal candidate events at low energy [14]. We define OT χ 2 as the reduced χ 2 computed between the triggered event and a cubic spline of the filtered average pulse obtained from the calibration γ rays at 2615 keV. This parameter is sensitive to the shape of the expected signal, suppressing pulses with shape deviating from the nominal one.  Fig. 2 (color online) A typical distribution of OT χ 2 as a function of energy for one BoDs. Physical events due to the particle interactions in the TeO 2 crystals are distributed in a nearly-horizontal band around OT χ 2 ∼ 1. Non-physical events such as electronic noise and tower vibrations, as well as particle interactions in the thermistors, follow an oblique distribution. The green solid line corresponds to the 90th percentile of the OT χ 2 distribution (OT χ 2 90% ) calculated in the region [35 − 50] keV and OT χ 2 < 10 using physics data. The magenta (orange) dashed line corresponds to the 90th (50th) percentile calculated using calibration data in the region [100 − 500] keV and OT χ 2 < 10, assuming linear dependence on energy. Red and blue dashed vertical lines represent the trigger threshold E 99% trig and the analysis threshold E thres , respectively. Fig. 2 shows a typical OT χ 2 distribution as a function of energy for the triggered events which pass the selection criteria described in Sec. 4.1. Between the OT trigger level and E 99% trig there can be a leakage of baseline noise. Physical events due to particle interactions in the TeO 2 crystals scatter around OT χ 2 ∼ 1 forming an almost horizontal distribution, while spurious events due to electronic noise or particle interactions in the thermistors follow an oblique distribution with OT χ 2 values as high as ∼ 10 4 at 200 keV. Pile-up events lie between the two bands. We select legitimate physical events with a requirement on the OT χ 2 parameter, and we evaluate the selection efficiency by counting the number of events before and after the cut in a region free of noise. As is evident in Fig. 2, OT χ 2 has a slight dependence on energy; this dependence is more or less pronounced depending on the bolometer, but mostly is imperceptible below 100 keV. Assuming no energy dependence at low energy in the range between 10 and 50 keV, we compute the selection efficiency in the region with OT χ 2 < 10 and energy in the range 35 − 50 keV, where the statistics are higher and the noise leakage is negligible. We choose the values of the selection to achieve 90% efficiency and calculate it as the 90th percentile of the OT χ 2 distribution (OT χ 2 90% , green solid line in Fig. 2). The selection efficiency is computed independently for every BoDs, and the uncertainty, evaluated taking into account the statistical fluctuation in counting for each BoDs, ranges from 0.3 to 1%. We have investigated the OT χ 2 dependence on energy using calibration data, as the low statistics above 60 keV make the results significantly uncertain in background data. The behaviour is well described by a linear fit up to 500 keV and the OT χ 2 90% value at 35 − 50 keV agrees with that calculated in background data for bolometers featuring low rate during calibration, like the one on Fig. 2, where the magenta dashed line corresponds to 90% efficiency and the orange dashed line to 50% efficiency. However, in general the larger pile-up probability during calibration runs shifts the OT χ 2 90% selection upwards with respect to the value in background runs, so in the following we use the value calculated in background and assume no energy dependence down to threshold. In order to estimate the uncertainty related to the choice of an energy independent cut efficiency, we have calculated the counting rate below 35 keV after the cut assuming the same energy dependence measured in calibration, being the difference with respect to the energy-independent selection lower than the statistical error.
The validity of the selection efficiency computation relies on the assumption that the region 35 − 50 keV is free of noise and the shape of the OT χ 2 distribution does not change at lower energies. To verify this hypothesis, we compare the OT χ 2 90% selection with the 50th percentile selection (OT χ 2 50% ), assuming their selection efficiencies are 90% and 50%, respectively. If there exists a significant noise contribution with the OT χ 2 90% selection, the efficiency-corrected spectra would differ, as noise rejection is stronger for the OT χ 2 50% selection. The residual spectrum is shown in red in Fig. 3. The selection efficiency corrected rate difference between the two spectra is compatible with zero down to ∼ 25 keV. Below ∼ 25 keV the rate difference increases, suggesting the presence of noise in the data. In fact, the noise contribution for most bolometers overlaps the physicalevents band at energies directly above E 99% trig . In order to avoid noise contribution in the spectrum we set the most stringent analysis energy threshold E thres >E 99% trig , independently for every BoDs, as described in Sec. 4.3.

Energy threshold determination
The CUORE-0 cryostat at LNGS was more noisy than the R&D cryostat in which the CCVR bolometer performance tests were performed. It means that in CUORE-0, at E = E 99% trig we mainly trigger noise events. The pulse shape parameter OT χ 2 presented in Sec. 4.2 provides powerful discrimination between physics and spurious events, but the two populations overlap as the energy decreases.
In order to avoid a leakage of spurious events in the data we set an analysis threshold (E thres ) at the minimum energy where the populations are well separated for each BoDs. Specifically, we perform a Kolmogorov-Smirnov (KS) test to quantify the similarity between the OT χ 2 populations in different energy slices with respect to the OT χ 2 distribution at 35 − 50 keV and OT χ 2 <10, the same pure signal sample region used to calculate the OT χ 2 selection efficiency. Starting from E 99% trig (the vertical red dashed line in Fig. 2), we compare the distribution of OT χ 2 in 4 keV-width energy windows with the reference distribution, for OT χ 2 < 10. We set E thres at the lower edge of the energy range that provides KS probability larger than 0.1. This method is not valid when the reference sample is contaminated with noise. We ensure the validity of the reference by requiring OT χ 2 90% < 6. The value 6 is obtained from the OT χ 2 90% distribution of all BoDs as the 2σ above the median. We discard 37 BoDs that do not fulfill this requirement from a total of 490. Once the KS threshold was fixed, the technique described worked for all 490 BoDs without manual adjustments, making it suitable for an O(1000) bolometers experiment.
Final analysis thresholds range from 8 to 35 keV, with only 16 BoDs having a threshold lower than 10 keV, thus not being representative of the whole data-taking. Therefore, we set 10 keV as the minimum CUORE-0 energy threshold. The exposure ranges from 1.6 kg · yr at 10 keV up to 18.56 kg · yr at 35 keV (see inset of Fig. 5). We verify that the noise acceptance is negligible with the same procedure used in Sec. 4.2. As shown as the blue band in Fig. 3, the 90%-50% residual with the analysis thresholds is compatible with zero down to 10 keV.

Anti-coincidence requirement
The last event selection criterion for the low energy rare event searches is anti-coincidence; i.e., no signal events are triggered in other bolometers in a certain temporal window. We use a coincidence window of ±100 ms, 20 times wider than that used for the standard 0νβ β decay analysis [20], due to the larger difference in characteristic rise time between low and high energy events. To evaluate the event loss due to random coincidences between physical events and unrelated events on another bolometer (anticoincidence selection efficiency) we use the 1461 keV γ ray peak in the single crystal energy spectrum. This peak, coming from 40 K EC, does not belong to any cascade, so the only true coincident event is the ∼ 3 keV X-ray from the Ar de-excitation, which is below our threshold. Counting the number of events in the 1461 keV peak of the single crystal spectrum before and after the selection we find the anti-coincidence selection efficiency to be 99.2±0.3%. We combine this efficiency with the 90% event selection efficiency on OT χ 2 to obtain the total detection efficiency. The uncertainty on the efficiency is BoDs dependent.

Energy calibration
During the calibration runs, the 232 Th sources are outside the cryostat, so the γ rays pass through a 1.4-cm-thick ancient Roman lead shield before reaching the detector. Consequently, the peaks in the low energy region of the spectrum are highly attenuated, and only those between 511 and 2615 keV are clearly visible and used to calibrate the energy response of each bolometer. As stated in Sec. 2, in the standard 0νβ β decay analysis we use a second-order polynomial with zero intercept to fit the reconstructed peak positions in the calibration spectrum to their nominal energies. To improve the energy resolution at the Q-value of 0νβ β decay, a combination of four energy estimators with slightly different calibration coefficients are used in the final analysis presented in [8]. For the low energy analysis, we use a single energy estimator to avoid complexity.
The calibration uncertainty in the 0νβ β decay ROI is 0.05 ± 0.05 (stat.) ± 0.09 (syst.) keV [20], but this value is energy dependent. To validate the extrapolation to energies below 100 keV, we use the characteristic Te X-rays that can follow a γ ray interaction, which can escape the crystal and be detected in another crystal. These events can be selected by requiring a coincidence in an adjacent crystal. The most intense X-rays are eight K-shell peaks ranging from 26 to 32 keV (see Table 2). Due to the energy resolution, these peaks are noticeable as a main and a secondary peak around 27 and 31 keV, in both the calibration and background spectra.
To determine their reconstructed energies, we fit the region from 22 to 34 keV with an eight-Gaussian line shape plus a linear background, where all of the Gaussians are constrained to have the same width. The relative intensities and positions of each Gaussian are fixed with respect to the main K α1 peak using nuclear data from [23]. In order to take into account any possible discrepancy in the relative intensities of the peaks arising from systematic effects in the detector, we determine these intensities with a Monte Carlo (MC) simulation based on the Geant4 package [24] (version 4.9.6.p03, see [25] for details) that includes the bolometerdependent energy thresholds and the analysis coincidence window. Fig. 4 displays the fit results of both calibration and physics data. The most intense K α1 peaks for calibration and physics data are measured to be 27.60 ± 0.05 keV and 27.65 ± 0.13 keV, respectively. The corresponding residuals with respect to the nominal energy are 0.13 ± 0.05 keV and 0.18 ± 0.13 keV, respectively. The latter indicates the systematic upward shift on the energy scale in the physics data, and we take into account its impact on the WIMP sensitivity reported in Sec. 6.
The difference in peak positions between the nuclear data and the simulation is found to be less than 0.01 keV. Nevertheless, the amplitude ratio K β 1 /K α1 in the MC simulation is 0.27 instead of 0.17 from Table 2 due to the strong change in X-ray attenuation length between 27 and 31 keV. This effect is also appreciable in the CUORE-0 data; as shown in Fig. 4 (upper panel) the 31 keV peak is underestimated by the fit function. Leaving the relative intensity K β 1 /K α1 as a floating parameter improves the goodness of the fit and reproduces the relative intensities estimated by MC, but the position of the K α1 peak does not change within the uncertainty. Fig. 5 shows the low energy spectrum of CUORE-0, using the selected BoDs with the event selection criteria described in Sec. 4 and the detection efficiency. The background rate above 50 keV is 0.05 counts/(kg · keV · day), consistent with the results obtained in the standard 0νβ β decay analysis. Below 50 keV, the background rate increases substantially to 1.7 counts/(kg · keV · day) at 10 keV; it is, however, two times lower than the background rate measured with the four best Cuoricino bolometers for which thresholds below 10 keV were attained [16].  The most noticeable feature in the spectrum is a peaklike structure around 31 and 37 keV. Given that this structure is observed in all the bolometers and it was also present in the Cuoricino background, its origin is likely physical. The current background prediction based on our MC simu-lation [25] does not fully account for the low energy spectrum including this peak-like structure, and further investigation is on-going under the hypothesis that the contamination is due to materials facing the detectors (e.g. copper shielding...). We expect to have better insights on this peaklike structure with CUORE data, where inner bolometers mostly face other bolometers and not the copper shielding. Through the comparison between the innermost and outermost bolometers, we may be able to attribute the origin of this peak-like structure to a certain process.

CUORE sensitivity to WIMP annual modulation
In this section we present the sensitivity of CUORE to the annual modulation in the detection rate induced by dark matter in the galactic halo. We restrict our analysis to WIMPs interacting with the target nuclei in the detector via elastic scattering off nuclei; for this reason, we present a study of the nuclear quenching factor of TeO 2 obtained in CUORE-0.

TeO 2 nuclear quenching factor
One of the prerequisites to perform a WIMP dark matter search is a good understanding on the low energy response of the detector for both nuclear recoils (NRs), produced by WIMPs or background neutrons, and electronic recoils (ERs), produced by electromagnetic backgrounds. The nuclear quenching factor is defined as the ratio of the measured signal generated by a NR to that generated by an ER depositing the same energy in the detector, and depends on the energy and recoiling nucleus. Given that any energy conversion in the TeO 2 bolometers finally produces signal through phonons, the nuclear quenching factor in the bolometers is expected to be close to one. The nuclear quenching factors of several recoiling nuclei in TeO 2 have been measured previously using the daughter nuclei of the α decays from 224 Ra, 220 Rn, 216 Po, 212 Po, and 212 Bi in the energy range between 100 and 170 keV. The result was found to be 1.025 ± 0.01 (stat) ± 0.02 (syst) [26]. We exploit the same technique and estimate the nuclear quenching factors using the daughter nuclei of some α emitters at energies around 100 keV.
Specifically, we measure the recoiling energy of the daughter nuclei following α decays of 210 Po, 222 Rn, 224 Ra, and 218 Po present in the CUORE-0 crystal surfaces where either the α particle or the daughter nucleus escapes and is detected in an adjacent crystal. We tag these events by requiring coincidence in two crystals with a total energy within some tens of keV of the Q-value of the decay. Then we fit the spectrum of the recoiling nuclei with an asymmetrical Gaussian function with a smooth power-law tail relative to the mean to obtain the peak position. Table 3 sum-marizes the obtained nuclear quenching factors for the selected recoiling nuclei. While the nuclear quenching factor obtained from 218 Po ( 214 Pb and 220 Rn) is close to unity, we notice that the one obtained from 206 Pb exhibits significant deviation from unity. Acknowledging that these recoiling nuclei are surface events and energy losses might happen at the surface, we use unity as the nuclear quenching factor to set the energy scale of WIMPs in the following analysis, and conservatively estimate an uncertainty of nuclear quenching factor as 7% using the largest deviation from unity observed by 206 Pb. Its impact is integrated in Fig. 6 to report the WIMP sensitivity of CUORE.

WIMP sensitivity of CUORE
The sensitivity of CUORE-0 to annual modulation of WIMPs is limited by its relatively small exposure. However, we can use results of CUORE-0 to estimate the CUORE sensitivity assuming the same background rate and energy thresholds. This is a conservative hypothesis since a significant background reduction is expected to be achieved in CUORE thanks to its close-packed detector configuration and the careful selection of radio-pure detector materials [25].
TeO 2 is an interesting DM target, as it combines a heavy nucleus (tellurium), which provides a large scattering amplitude (assuming coherent interaction, that scales as A 2 ) and a light one (oxygen) to enhance sensitivity in the low-mass WIMP region. To calculate the expected WIMP rate in the detectors, we follow the commonly used analysis framework for WIMP direct detection [27,28]. We consider only the spin-independent (SI) contribution since the spin-dependent contribution is comparatively reduced in TeO 2 ; the main isotopes with non-zero nuclear spin are 125 Te and 123 Te, with isotopic abundance of 7.1% and 0.9%, respectively. We assume coherent isospin-invariant coupling and the Helm model [29] for the nuclear form factors. Under these assumptions, the generic WIMP is completely determined by its mass m W and SI WIMP-nucleon cross section σ SI . For the velocity distribution of dark matter, we use the standard halo model (SHM) [30] commonly adopted for comparisons of direct detection experiments. Consequently, the annually modulating WIMP recoil rate due to the motion of the Earth around the Sun can be approximated using a constant term S 0 plus a cosine-modulated term S m , as given by where ω = 2π/yr and t 0 is around June 1.
To obtain the sensitivity to annual modulation of WIMPs, we scan over the WIMP parameter space (m W , σ SI ) looking for the region at which a WIMP interaction would produce an annual modulation in the detection rate over the measured background at a certain confidence level (C.L.). For each (m W , σ SI ) we generate 100 toy MC simulations and for each MC spectrum, we perform a maximum likelihood (ML) analysis for both the annual-modulation (AM) and the absence of modulation (null) hypotheses. We quote the significance of the modulation as the log-likelihood ratio of the best fits χ 2 =2log(L AM /L null ). The likelihood L AM is calculated using the probability density function (PDF) where M det is the target mass, ε BoDs (E,t) is the BoDsdependent detection efficiency, and φ b is the background PDF, which we model with a Chebychev polynomial with coefficients a i and for which we do not consider any temporal dependence. The likelihood L null is calculated from φ b alone. We choose the ROI to perform the analysis as 10 − 28 keV, which excludes the peak-like structures above 30 keV shown in Fig. 5. Given that the differential rate of WIMPs quasi-exponentially falls as a function of energy, most of the signal is contained at the low energy and the expected contribution to the WIMP sensitivity from 30−60 keV is negligible compared to that from 10−28 keV. We consider a target mass of 742 kg and the scheduled 5 years of data-taking with 75% duty cycle, accounting for the calibration time. Based on the CUORE-0 energy threshold, we use 10 keV but we also show the sensitivity that could be attained under the more optimistic hypothesis that we reach a 3 keV threshold as demonstrated in the CCVR experiment with a linear extrapolation of the CUORE-0 background to lower energies. Fig. 6 shows CUORE sensitivity requiring a 90% C.L. in 90% of the toy-MC experiments. The results are consistent with those obtained with a pure statistical calculation following [36]. This figure assumes a WIMP local density ρ χ = 0.3 GeV/c 2 , local circular velocity v 0 = 220 km/s,  Fig. 6 (color online) 90% sensitivities on the SI elastic WIMP-nucleon cross section as a function of WIMP mass of CUORE, assuming 5 years of data-taking with 75% of duty cycle and 10 keV threshold (red), as well as 3 keV threshold (blue). Uncertainty on the energy scale dominated by the nuclear quenching factor is taken into account. DAMA/LIBRA positive signal reported in [31] is shown as yellow/dark yellow/orange islands. The results from CRESST-II (dashed green) [32], CDMS Lite (dashed red) [33], XMASS (dashed violet) [34], and LUX (black solid) [35] are also shown. galactic escape velocity v esc = 650 km/s, and orbital velocity of the Earth around the Sun v orb = 29.8 km/s. Uncertainty on the energy scale dominated by the nuclear quenching factor is taken into account. For comparison, we also show the 5σ , 3σ and 90% C.L. regions resulting from a ML analysis reported in [31] on the DAMA/LIBRA annual modulation positive result [37,38] using the same parameters for the SHM. Thanks to the 741 kg of target mass of CUORE, we expect to achieve the sensitivity required to fully explore the parameter region implied by the DAMA/LIBRA positive annual modulation signal with 5 years of data-taking. Other recent experimental results from CRESST-II, CDMS Lite, XMASS and LUX [32][33][34][35] are also shown. The results from CRESST-II, CDMS, and LUX were obtained using v esc = 544 km/s. The impact of using v esc = 544 km/s instead of v esc = 650 km/s for CUORE sensitivity is found to be less than 10 −5 pb at 6 GeV. Also for the other experiments only a minor impact of the escape velocity on the exclusion limit is expected.

Summary
We have presented the analysis techniques developed for low energy rare event searches with CUORE and their validation using the data acquired with the CUORE-0 experiment. We have optimized the software trigger developed in previous CUORE prototypes, removing an intrinsic dead time that prevented the algorithm from reaching 100% efficiency, and designed a protocol to periodically monitor the efficiency by injecting low energy reference pulses at the end of every dataset. With the new trigger, we have reduced the CUORE-0 trigger thresholds from several tens of keV to values between 4 and 12 keV.
We have also demonstrated that a pulse shape analysis can efficiently select legitimate physics events in TeO 2 bolometers against spurious ones at energies below 100 keV. In addition, we have developed a technique, scalable to an experiment with one thousand bolometers, to independently establish the analysis threshold of each bolometer in each dataset. In CUORE-0 the analysis thresholds range between 8 and 35 keV. Using characteristic X-rays from Te, we have found the energy scale shift to be 0.18 ± 0.13 keV upward at ∼ 27 keV.
After the data and event selection, the CUORE-0 background rate ranges from 1.7 counts/(kg · keV · day) at 10 keV to 0.05 counts/(kg · keV · day) at 50 keV, two times less than that attained with the best Cuoricino bolometers. Nevertheless, the low energy spectrum requires further investigation including the explanation of the peak-like structures between 30 and 40 keV. We use the nuclear quenching factors of TeO 2 obtained by tagging the recoiling daughter nuclei of α decays in the CUORE-0 data to estimate the uncertainty of the WIMP energy scale. We incorporate it to report the CUORE sensitivity to WIMP annual modulation.
CUORE will search for low energy rare events such as solar axions, WIMP dark matter in the galactic halo, or coherent scattering of galactic supernova neutrinos using the analysis techniques presented in this paper. In particular, we expect to reach a sensitivity to annual modulation of WIMPs sufficient to fully explore the parameter region indicated by the positive annual modulation signal of the DAMA/LIBRA experiment with 5 years of CUORE data-taking.