Analysis of cryogenic calorimeters with light and heat read-out for double beta decay searches

The suppression of spurious events in the region of interest for neutrinoless double beta decay will play a major role in next generation experiments. The background of detectors based on the technology of cryogenic calorimeters is expected to be dominated by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document}α particles, that could be disentangled from double beta decay signals by exploiting the difference in the emission of the scintillation light. CUPID-0, an array of enriched Zn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{82}$$\end{document}82Se scintillating calorimeters, is the first large mass demonstrator of this technology. The detector started data-taking in 2017 at the Laboratori Nazionali del Gran Sasso with the aim of proving that dual read-out of light and heat allows for an efficient suppression of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document}α background. In this paper we describe the software tools we developed for the analysis of scintillating calorimeters and we demonstrate that this technology allows to reach an unprecedented background for cryogenic calorimeters.


Introduction
As of today, we do not know any process in nature that violates the total number of leptons L or the number of baryons B, even if the Standard Model of Particle Physics does not predict the conservation of such quantities. On the contrary, the Standard Model predicts, also non-perturbatively, the conservation of a simple combinations of these numbers: B-L [1]. A violation of this quantity would be a clear hint of physics beyond the Standard Model and this is one of the reasons motivating the endeavor to search for a never-observed process: the neutrino-less double beta decay (0νDBD) [2,3]. This process is a hypothesized nuclear transition in which a nucleus decays with no neutrino emission: The importance of 0νDBD resides also in the fact that it can occur only if neutrinos coincide with their anti-particles, so its detection would allow to establish the ultimate nature of this elusive particle. Finally, the measurement of the 0νDBD half-life T 0ν 1/2 would provide some insight into the absolute mass of neutrinos [4].

Scintillating cryogenic calorimeters
The analysis techniques described in this paper apply to experiments using the technology of cryogenic calorimeters (historically also called bolometers). A cryogenic calorimeter is made by a temperature sensor coupled to a crystal, which acts as energy absorber. The interactions in the crystal release an amount of energy that gives rise to a sizable temperature variation ( T ∝ E/C), provided that the crystal thermal capacity C is low enough. To this aim, the crystals are cooled at cryogenic temperatures (about 10 mK). The main advantages of this technique, originally proposed by Fiorini and Niinikoski [5], are the energy resolution (as good as 0.1%) and an efficiency on 0νDBD larger than 80%. Furthermore, the crystals can be grown with high intrinsic radio-purity starting from most of the emitters of interest for 0νDBD.
The first tonne-scale experiment based on cryogenic calorimeters is the Cryogenic Underground Observatory for Rare Events (CUORE [6]), now in data-taking at Laboratori Nazionali del Gran Sasso (LNGS) in Italy. The analysis of the first months of data (corresponding to an exposure of 86. 3 kg year) proved that the detector can reach the target energy resolution and background, and allowed to place a 90% C.L. lower limit of T 0ν 1/2 ( 130 Te) >1.3×10 25 years alone, and of T 0ν 1/2 ( 130 Te) >1.5×10 25 years combined with its ancestors Cuoricino and CUORE-0 [7][8][9].
Today, the CUPID (CUORE Upgrade with Particle IDentification [10,11]) interest group is defining the strategy for a future upgrade of CUORE that will allow to increase the sensitivity on the half-life of 0νDBD above 10 27 years [12][13][14][15].
The main challenge for the CUPID project will be the development of a background-free detector at the tonne-scale level. The first important milestone is the abatement of the dominant source of background of CUORE, i.e. α particles produced by the materials constituting the detector structure [16]. It was proved [17] that the α interactions can be rejected by coupling each calorimeter with a second detector, specialized in the measurement of the scintillation light emitted by the interactions in the crystal. Unfortunately, TeO 2 does not scintillate at cryogenic temperatures [18]. For this reason, the LUCIFER [19] and LUMINEU [20] collaborations focused on the development of a new class of scintil-lating crystals based on 0νDBD emitters characterized by a high Q-value. Indeed, choosing 0νDBD candidates with high Q-value, such as 82 Se or 100 Mo, provides a natural reduction of the background contribution from environmental γ 's, that drops above the 2.6 MeV γ -line of 208 Tl. An extensive R&D activity allowed to characterize the properties of several compounds grown with such emitters, like ZnSe [21,22], ZnMoO 4 [23][24][25][26][27][28][29] or Li 2 MoO 4 [30][31][32][33][34]. These R&Ds demonstrated that the simultaneous read-out of light and heat in scintillating calorimeters enables a very effective suppression of the α background. To prove the potential of this technology on a medium-scale experiment, we designed and constructed the CUPID-0 detector [35], now in datataking at LNGS.

The CUPID-0 detector
The 0νDBD emitter chosen by the CUPID-0 collaboration is 82 Se. This isotope features a Q-value (2997.9 ± 0.3 keV [36]) well above the 2615 keV end-point of the natural radioactivity, and its half-life for the 2νββ mode is long enough (T 2ν 1/2 = (9.2 ± 0.7)×10 19 years [37]) to prevent background from pile-up in a tonne-scale experiment.
The Se powder was enriched from its natural abundance to 95% [38], and embedded in 24 Zn 82 Se crystals (plus 2 natural ZnSe crystals) [39]. The total mass of the 24 Zn 82 Se crystals amounts to 9.65 kg (5.13 kg of 82 Se), while the two natural crystals have a total mass of about 850 g (42 g of 82 Se). The ZnSe crystals are surrounded by a VIKUITI multi-layer reflecting foil produced by 3M, and arranged in 5 towers using a NOSV copper 1 structure and PTFE supports. Each ZnSe is interleaved with two light detectors (LD). These LDs consist of disk-shaped Ge crystals (170 μm thick and 4.4 cm in diameter) similar to those described in Ref. [40]. One of the best ways to obtain high-performance LDs at 10 mK consists of operating also the LDs themselves as cryogenic calorimeters: photons impinging on the LD increase its temperature and are recorded as thermal pulses.
To convert the energy deposits in ZnSe crystals and LD in readable voltage signals, each crystal was equipped with a Neutron Transmutation Doped (NTD) Ge thermistor [41] using a semi-automatic gluing system. In addition, a Si Joule heater was attached to both detectors to inject a reference pulse, which allows to correct for thermal drifts [42,43].
The CUPID-0 detector is hosted in the same 3 He/ 4 He dilution refrigerator that was used for the CUORE-0 experiment, after some major upgrades to the electronics and to the vibration reduction system. The reader can find in Refs. [35,44] a more extensive description of the cryogenic facility, elec-tronics and data-acquisition, as well as more details about the detector construction, operation and optimization.

Data collection and production
An interaction in the ZnSe crystal results in an amplified signal with amplitude ranging from tens to hundreds of mV, a rise time of 10 ms and a decay time of 40 ms. Due to the smaller detector sizes, the LD signals are usually faster, with rise-time of a few ms and decay-time of about 8 ms. Each channel can be biased, amplified and filtered using a dedicated read-out chain, which allows to optimize the amplification gain and the cut-off frequency of the anti-aliasing filter [45][46][47][48][49][50][51][52]. Due to the slow time-development of the recorded pulses and the low detector rate (2 mHz in physics runs), the data are digitized with sampling frequencies of 1 kHz for ZnSe and 2 kHz for LD, and the continuous data stream is transferred to disks for the off-line analysis. The data collection is made with a DAQ software package ("Apollo" [53,54]) that in the past was used for CUORE-0 and it is now being used by CUORE.
The trigger is software generated, and allows to use different algorithms according to the experimental needs. The data presented in this paper are processed with two triggers. For the ZnSe calorimeters we use a trigger algorithm with a channel-dependent configuration that fires when the signal derivative stays above threshold for a certain amount of time. For the LD we use simultaneously the derivative trigger and a second (off-line) trigger that forces the acquisition of the LD waveforms every time the ZnSe trigger fires. The implementation of the second trigger was motivated by the fact that most of the energy produced by an interaction is dissipated as heat in the ZnSe, while only a few % escapes the ZnSe crystal in the form of scintillation light: a ∼MeV deposit in the ZnSe crystal corresponds to about 10-100 keV in the LD (depending on the crystal as well as on the nature of the interacting particle). To prevent the loss of small (or noisy) light signals, when a signal is detected in the ZnSe we also associate to it the corresponding waveforms in the LD. In the following, we will use only this trigger for the LD. The derivative trigger is still run on the light detectors for future analyses (for example, to study events that are not produced by scintillation of ZnSe).
The complete data-stream of all channels recorded by the DAQ, as well as the trigger positions, are saved in NTuples based on the ROOT software framework. At the first stage, we convert the continuous data into acquisition windows of 5 s for the ZnSe crystals (4 s after the trigger and 1 s before) and 1 s for the LD (800 ms after the trigger and 200 ms before). The pre-trigger window is used to compute the baseline value, and thus the detector gain before the interaction occurred [55,56].
Other informations to be accessed during the off-line analysis, such as the geometrical configuration of the array, the correspondences between ZnSe and LD, the run type (physics, calibration, test...), possible time intervals that have to be rejected because of known problems (earthquakes, electronics problems, major underground activities) are stored in a PostgreSQL database.
Each physics run lasts about 2 days, and it is followed by a stop of a couple of hours to allow the liquid helium refill of the cryostat and the subsequent stabilization of the detectors. Approximately every month, we perform a calibration of 4 days with 232 Th sources. Since the most energetic γ line produced by 232 Th (2.6 MeV) is below the 82 Se Q-value, we exploit also other sources to characterize the energy region of the 0νDBD. To study the energy dependency of the shape parameters in the region of interest (Sect. 5), we use an Am:Be neutron source, emitting a broad distribution of γ rays up to several MeV. This calibration was made every time we modified the working parameters of the detectors, in order to prevent possible changes in the shapedependency on the energy. In the first year of data taking, we performed three Am:Be neutron source calibrations: one during the detector commissioning, one between the physics runs presented in this paper, and one immediately after (as a cross-check). Furthermore, we validated the 232 Th calibration with a 56 Co source, producing γ peaks well above the 82 Se Q-value (Sect. 7).
The collection of the initial plus final calibration and all the physics runs in between forms a DataSet. With the exception of the first two DataSets, devoted to the detector optimization, the percentage of live-time for physics analysis (thus excluding calibrations) exceeds 80%.

Heat pulses reconstruction
The conversion of the continuous data stream into NTuples containing all the quantities of interest is performed with a C++ based analysis framework ("DIANA") originally developed for Cuoricino. In this section we summarize the processing stages that allow to derive the parameters of interest of the heat pulses. Most of these analysis techniques are very similar to those developed by the Cuoricino, CUORE-0 and CUORE collaborations and are extensively described in Refs. [57][58][59][60].
The heat and light pulses are processed with a matched filter algorithm to suppress the signal frequencies mostly affected by noise and improve the reconstruction of the pulse amplitude [61,62]. This software filter requires (for each channel) a template for the detector response and the noise power spectrum, shown in Fig. 1. The average noise power spectrum is constructed by averaging hundreds of waveforms acquired during the entire DataSet with a random trigger. A The signal template is obtained by averaging hundreds of high amplitude events collected during the 232 Th calibrations, and aligned by the pulse maxima. In such a way, the random noise contributions are suppressed (see Fig. 1 left). In Sect. 8 we describe how the production of the signal template was improved to match the needs of scintillating crystals. After the matched filter we extract also some parameters related to the pulse shape: rise-time (time difference between the 90 and the 10% of the leading edge), decay-time (time difference between the 30 and 90% of the trailing edge), slope of the baseline before the pulse, delay of the position of the maximum of the filtered pulse with respect to the maximum of the template, and two shape parameters called Test Value Left (TVL) and Test Value Right (TVR), that correspond to the χ 2 value between the filtered signal template and the filtered pulse computed on the left and right side of the signal maximum, respectively: where y i is the pulse, A and i M its amplitude and maximum position, s i the ideal signal pulse scaled to unitary amplitude and aligned to y i , ω L (ω R ) the left (right) width at half maximum of s i . The signal amplitude computed with the matched filter is corrected for temperature instabilities by exploiting the periodic reference pulse injected with the Si heater [57]. After the correction we expect a residual instability negligible with respect to the noise fluctuations of the detector. The response of the ZnSe detectors is then equalized by energy-calibrating  In the last stage of the analysis, we compute time coincidences between ZnSe crystals. Rejection of coincidences between crystals plays a major role in the background suppression as, from GEANT-4 simulations, we expect 81.0 ± 0.2% of the 0νDBD events to be fully contained in a single crystal [63]. The coincidence window is optimized by selecting events produced by the 232 Th source in which two ZnSe crystals trigger with a total energy of 2615 keV, and is set to 20 ms. Given the counting rate of the detectors during the physics runs, we compute the probability of random coincidences among ZnSe crystals to be 1.7 × 10 −6 .
Finally, we estimate and remove the energy dependency of the shape parameters on both their absolute values and resolutions (see Fig. 2 left), that otherwise would limit our capability to define robust cuts on the pulse shape. To correct the energy dependency on a wide energy-range, we exploit the periodical calibration with the Am:Be neutron source, producing γ interactions (and thus particles with the same shape of the 0νDBD signal) up to several MeV. First of all, for each channel, the energy spectrum is divided in slices, for each of them the median and the MAD (median average deviation) of the considered shape parameter are evaluated. Then, we interpolated both the median and MAD points with a polynomial functions obtaining, in such a way, their trends in the whole energy range (up to 4 MeV), including the region where we would expect the 0νDBD signal. Finally, we use the parameters extracted from the fits to correct for the energy dependency in all the physics and calibration runs. Each shape parameter is scaled in such a way to be centered around zero (Fig. 2).
To monitor the effect of this normalization, we developed a software tool which checks that the scaled parameters do not depend on the energy, on the channel and on the measurement time. This analysis led to discard some of the shape parameters (such as the TVR) that are not stable enough to perform uniform cuts over a dataset, as they show a time

Selection of heat pulses and efficiency evaluation
We perform a first selection of thermal pulses by exploiting the shape parameters listed in Sect. 5. At this stage, we do not use the information provided by the light detectors, as the purpose is to reject spurious events, such as those barely affected by pile-up or electronics noise.
First, we exclude all the time-intervals that are marked as bad in the database because of known problems (see Sect. 4). The effect of this selection is a reduction in live-time by 1%. Moreover, we require the pulses to be triggered only by a single ZnSe, as expected from 0νDBD events.
As explained before, we exclude from the analysis the shape parameters that are not robust enough because of fluctuations in time, and we use only the decay-time, rise-time, baseline slope, delay and TVL.
To investigate the effects of cuts on the shape parameters, we study the γ peak of 65 Zn, a product of the activation of the Zn contained in the crystals, that decays via electron-capture with a half-life of 224 days and a Q-value of 1351.9 keV. This signature acts as a signal sample, while the side-bands close to the γ peak are chosen as background samples. The odd events are used to optimize the cut while the even ones to compute the selection efficiency.
We scan the distribution of each scaled shape parameter by cutting at different integer values. In Fig. 3, we report the efficiency on the γ peak of 65 Zn ( S ) and on its side-bands ( B K G ) as a function of the value at which we cut the scaled decay-time.
This plot shows that S is larger than B K G , proving that the choice of the signal/background samples was reasonable. The reason why they do not differ dramatically, is events that pass the pulse-shape cuts and the anti-coincidence cut. Bottom: events rejected by the pulse-shape cuts. We fit both the plots simultaneously with an unbinned extended maximum likelihood fit with two components: a Gaussian function and an exponential background using the RooFit analysis framework that the background sample contains also a large fraction of events due to the 2ν double beta decay that, as expected, is not affected by the shape cuts. When the cut becomes wide enough, both the efficiencies reach a plateau. To set the proper cut value, i.e. keep the highest efficiency on signal while suppressing the spurious events, we compute the ratio r = S / √ B K G (Fig. 3) and choose the cut in which r reaches the plateau.
As explained before, we evaluate the total efficiency of these cuts on the even events belonging to the γ peak of 65 Zn (Fig. 4). Even if we keep the same shape-cuts during the entire analysis, we compute the efficiency separately on each DataSet to account for possible time-variations of the shape parameters.
Weighting the efficiencies by the DataSet exposure, we obtain an average efficiency of 95 ± 2%, with a maximum variation of 6% across all DataSets. This value is cross-checked using events in which two crystals triggered that, given the negligible amount of random coincidences, can be considered as an almost pure sample of signal-like events. We obtain an efficiency compatible with the one evaluated on 65 Zn and constant from 300 to 2615 keV.
The energy region chosen for the analysis of the background is a 400 keV interval centered around the 82 Se Qvalue (2800-3200 keV). At higher energies, indeed, we expect the background to decrease, as the contributions from 214 Bi and 208 Tl (the dominant background sources) drop The red circle indicates the value extrapolated from the 232 Th calibration in Ref. [63] above 3200 keV. Therefore, further enlarging the analysis window would result in a lower background. The lowest bound of the interval was chosen to have a symmetric region around the Q-value and, at the same time, to avoid contributions from the 2615 keV photon or from the tail of the 2ν double beta decay.
Applying these cuts on the pulse-shape parameters, and requiring that each pulse is triggered by a single ZnSe, we obtain a background index of (3.6 ± 0.5)×10 −2 counts/(keV kg year).

Validation of the 232 Th calibration
Before introducing the information on the light detectors, it is worth observing that the energy calibration of the heat channels was cross-checked in a dedicated measurement. Indeed, since the Q-value of 82 Se exceeds the largest γ ray produced by the 232 Th source, we usually extrapolate the calibration function and the energy resolution in the region of interest. This procedure was used in the analysis of the 0νDBD reported in Ref. [63], in which the extrapolation at high energies resulted in an uncertainty of 3 keV on the Q-value, and an energy resolution of (23.0 ± 0.6) keV FWHM. In this paper we validate the 232 Th calibration by using a 17 dayslong measurement with a 56 Co source (T 1/2 ∼ 77.2 days) emitting γ rays up to ∼3.5 MeV. We apply to these runs the calibration coefficients derived for each ZnSe from the 232 Th calibrations, as done in a "standard" DataSet. We then fit the most prominent γ peaks using a double-gaussian model [63] and study the difference between the obtained position and their nominal energy, as well as their energy resolution (Fig.  5).
The residuals show a dependency on the energy that can be modeled with a parabolic function, resulting in a uncer-tainty on the position 0νDBD peak of 1.3 keV. This value is negligible compared to the energy resolution in that region, and proves that the choice of an uncertainty of 3 keV in the analysis of the 0νDBD reported in Ref. [63] was very conservative. The energy dependency of the energy resolution is modeled with a linear function. In the region of interest we obtain a FWHM of 22.5 ± 1.2 keV, fully consistent with the value extracted from the 232 Th calibration [63].

Reconstruction of light pulses
The first step for a correct reconstruction of the light pulses consists of generating a dedicated signal template for the matched filter. In the past, the signal response of the LD was made by averaging many pulses with good signal-to-noise ratio, obtained for example using a 55 Fe X-ray source permanently exposed to the detector. Nevertheless, this is not the best approach, as X-rays, α particles and electrons are characterized by a different time development of the light pulses, and constructing the ideal detector response on a class of events that is not similar to the one of 0νDBD can spoil the evaluation of the light shape parameters. This has a particular importance for CUPID-0 that, as explained in Sect. 9, takes advantage from the shape of the light pulses for particle identification. For this reason, we developed a new algorithm that selects only events with a shape similar to the one of 0νDBD. Exploiting the γ energy calibration made with the 232 Th source, we select heat pulses with energy of 1. 8-2.64 MeV. This energy interval is wide enough to provide a large sample of pulses without introducing non-linearities in the pulse shape. Moreover, we require the emitted light to be compatible with the one produced by scintillation of electrons to discard events produced by scintillation of α particles, or events with no associated light emission (electronics noise, interactions in the NTD Ge sensor). Finally, we reject spurious events, such as those affected by random pile-up, or those in which a second pulse was detected in the same acquisition window. The effect of the selection is shown in Fig. 6.
The events that pass all the selection cuts (a few hundreds for each ZnSe) are finally averaged to suppress the random noise contributions. For each ZnSe, we construct 3 signal templates: the template of the ZnSe itself, and the templates of the light recorded by the top/bottom LD. We stress that now, in contrast to the past [32], each LD has two different signal templates, corresponding to the light emitted by the top/bottom ZnSe (see Fig. 7).
The new structure of the signal templates demands for a new version of the matched filter algorithm with respect to the one described in Ref. [57]. From the new matched filter we extract again all the shape parameters of the light pulses. In the analysis presented in this paper we considered only the  Fig. 7 Left: a typical template of a LD response (black line) overlapped to a single light pulse acquired by the same detector (magenta line). The signal template was evaluated averaging hundreds of light pulses emitted by the bottom ZnSe crystal in order to suppress the random noise fluctuation. Right: a typical average noise power spectrum of a LD. The microphonic noise picks and the roll-off due to the anti-aliasing active Bessel filter are clearly visible parameters of the light detectors placed on top of the ZnSe crystals as the SiO coated face is more sensitive with respect to the other one. We observe that the parameter that provides the best particle identification is the TVR (Sect. 6) that in the following will be called Shape Parameter, or S P, for simplicity. A preliminary study indeed, allowed to infer that the background rejection obtained with S P over-performs the one obtained with the light yield alone.
Thanks to the new algorithm used for the light signals analysis, the TVR of the light pulses does not show channeldependent behaviour. Furthermore, since the LDs are operated at a slightly higher temperature with respect to the ZnSe detectors their working points result very stable over the time. Therefore, the normalization procedure is not needed for this parameter because it turns out to be very stable and repro-ducible both over the channels and time, as shown in the next section.
Finally, we improve the evaluation of the amplitude of the light pulses that, given the worse signal-to-noise ratio with respect to heat pulses, could be affected by larger uncertainties. This problem is corrected by measuring the matchedfiltered amplitude of the light pulses at a fixed time-delay with respect to the ZnSe scintillation that triggered the event (see Ref. [64]). The main difference with respect to the algorithm described in that paper is the calculation of the time delay. In Ref. [64], the time-delay was computed as the median of the time intervals between the heat pulses and the corresponding light pulses. In CUPID-0 we compute the time delay using the filtered signal templates of heat and light. This algorithm does not give a more precise evaluation of the delay, but it is more simple and fast to implement as it requires only the templates of the signals.
In contrast to the ZnSe channels, it is not possible to energy-calibrate the amplitude of the LD using the 232 Th strings placed in the external part of the cryostat. The energy of the scintillation light produced by interactions in the ZnSe crystals usually ranges from a few keV to tens of keV; γ 's of this energy are too weak to penetrate the external shield of the refrigerator. In the past, this problem was overcome by depositing an X-ray source on a support permanently exposed to the surface of the LD. In CUPID-0, we decided to avoid the presence of sources to be conservative from the point of view of the radioactivity, considering also that the absolute energy scale of the light pulses is not an important information, as long as the light emitted by different particles permits their discrimination.

Alpha background rejection
The simultaneous read-out of the heat and light emitted by the ZnSe allows to reduce the background in the energy region of interest without affecting the signal efficiency. After the selection of "good" thermal pulses (Sect. 6), we add the information provided by the LD. We make a first selection by requiring each pulse to be associated to a non-zero light emission, to discard events that interacted in the temperature sensor, or electronics spikes and other spurious events not rejected by the (not aggressive) pulse-shape cuts. Then, we study the shape of the light pulse S P as a function of the heat released in the ZnSe crystal (Fig. 8).
From a qualitative point of view, it is clear that the population of α events, that would produce a background of about 2×10 −2 counts/(keV kg year) in the region of interest, can be clearly distinguished and rejected. However, the absence of peaks close to the 82 Se Q-value, as well as the energydependency of the S P, prevents a simple estimation of the To compute the signal efficiency, we select a pure sample of β/γ events in the ZnSe detectors: the ones that come from the electromagnetic showers produced by muons interacting in the materials that surround the detector. These events are produced in cascades, resulting in simultaneous triggers in several ZnSe crystals. For this reason, we select events in which at least five detectors triggered, obtaining a sample of 113 events. The sample is further selected by imposing a reasonable value for the detected light (larger than the noise fluctuations of the LD and smaller than the maximum light emitted by α particles). Indeed, a muon could cross the light detector and ionize it before producing the γ cascade, and the effect of the ionization would be an unpredictable value of the S P. Finally, we remove events in which the light detectors feature more than one pulse in the same acquisition window that, again, could lead to a wrong calculation of the S P. Figure 9 shows the distribution of the S P of the selected events which follows, as expected, the distribution of the electrons, and thus can be considered a good sample for the signal. Looking at this distribution we set the cut S P < 6, as this is the smallest value that yields a 100% efficiency on the signal. This cut allows to reduce the background in the analysis region to (1.5 ± 0.3)×10 −2 counts/(keV kg year).
Finally, we study the probability of mis-identifying an α interaction by selecting events with energy between 4 and 8 MeV. The distribution of the S P of these events can be modeled using a Gaussian function with mean value μ =13.33 ± 0.01 and σ = 1.38 ± 0.01. The probability for events that follow this distribution to occur below S P = 6 (the selected cut) is 5×10 −8 , proving that the probability of a misidentification is negligible. Even if it provides a satisfactory description of the data, the choice of modelling this distribution with a Gaussian function is not supported by physics considerations and could therefore lead to an underestimation of the background events in the region of interest, due to the presence of some outliers. Nevertheless, the number of events far from the cluster of alpha events is very small, proving that the large majority of the background events can be efficiently rejected.
Summarizing, the combination of light and heat allows to suppress the α background by almost a factor three without affecting the signal efficiency. Moreover, the alpha rejection capability already matches the requirements of next generation experiments, such as CUPID, in which the background must be close to zero at the tonne-scale level.

Improving the time veto
The identification of α particles down to low energies can help in reducing also the β/γ background. Indeed, one of the Table 1 Summary of the background index (counts/keV/kg/y) and signal efficiency averaged on the DataSets exposure, measured in the region 2800-3200 keV with a ZnSe exposure of 3.44 kg·y (1.34 × 10 25 emitters·y). Uncertainties are reported at 68% C.L.. First row: events that pass the cuts on the heat described in Sect. 6. Second row: the events are further selected requiring that the shape parameter of the light is consistent with interactions of electrons (α rejection) as described in Sect. 9. Third row: we added a time-veto of 3 half-lives after the detection of an α particle with energy compatible with the Q-value of 212 Bi. Fourth row: we added a time-veto of 3 half-lives after the detection of an α particle with energy larger than 2 MeV. Last row: we report the total efficiency, including the data selection efficiency computed as in fourth row, the trigger efficiency, and the electrons containment efficiency of (81.0 ± 0.2) % (see Ref [63] most worrisome background sources in the region of interest is 208 Tl, an isotope belonging to the 232 Th chain that decays via β/γ with a Q-value of about 5 MeV. Nevertheless, the β/γ interactions produced by 208 Tl can be efficiently rejected by exploiting the time coincidence with its parent, 212 Bi. This isotope decays to 208 Tl with the emission of an α particle (Qvalue of ∼6207 keV), and 208 Tl subsequently decays with a half-life of about 3.05 min. Therefore, the β/γ background from 208 Tl can be suppressed by vetoing the detectors for a few minutes after the occurrence of an α-like event with an energy corresponding to the 212 Bi Q-value. This technique was already exploited in the past with satisfying results [22,26]. To show its effect in CUPID-0, we report in Fig. 10 the high energy region of the β/γ spectrum. This spectrum is obtained applying cuts on the pulse-shape, on the number of triggering ZnSe crystals (Sect. 6), and the α particles rejection (Sect. 9) and, as explained in the previous sections, results in (1.5 ± 0.3)×10 −2 counts/(keV kg year).
From a Monte Carlo simulation taking as input the crystals contaminations (the reader can find in Ref. [35] the α spectrum of the detector and the activities of the main peaks), we expect a large fraction of this background ((1.10±0.2)×10 −2 counts/(keV kg year)) to be dominated by 208 Tl.
As shown in Fig. 10, the time-veto is very effective in the abatement of high-energy β/γ events. By applying a veto of 3 half-lives (3×3.05 min) after the detection of an α particle at the Q-value of 212 Bi, the background reaches a value of 5.1 +2.4 −2.0 ×10 −3 counts/(keV kg year) with a dead-time of 1%.
The innovative idea in CUPID-0 consists in enlarging the window for the identification of an α produced by 212 Bi down to much lower energies, by exploiting the excellent discrimination capability between α's and electrons. When the 212 Bi decay occurs inside the crystal, indeed, it releases the whole decay energy (α + nuclear recoil) inside the detector, producing a characteristic peak at the Q-value of the transition.
On the contrary, when the decay occurs on the crystal surface, or on the surface of the materials surrounding the crystal, the α particle can loose a variable fraction of its initial energy, resulting in a low-energy deposit inside the detector. By exploiting the possibility of disentangling α particles from electrons through the read-out of the scintillation light, we can tag also 212 Bi interactions that do not produce a peak at the Q-value of the decay. For this purpose, we select the possible 212 Bi parents by choosing events with energy larger than 2 MeV and light S P between 7 and 25 (Fig. 8).
In Fig. 10 we compare the background obtained with a veto that exploits only the 212 Bi peak at the Q-value, with the background obtained with a veto that exploits all the α's down to low energies. In the latter case, the background in the region of interest becomes 3.6 +1.9 −1.4 ×10 −3 counts/(keV kg year) with a dead-time of 2.6%.

Summary and conclusions
In this paper we presented the analysis methods to exploit the simultaneous read-out of light and heat to suppress the background of cryogenic calorimeters. We showed a new method to create the signal templates for the matched filter, that allows to simplify the data processing and, at the same time, to select pulses as similar as possible to the expected signal. We presented a technique to discriminate against α particles and a new method to estimate the efficiency on the signal and on the background rejection. Finally, we discussed a new time-veto for the suppression of the 208 Tl background, that can deal with both internal and surface contaminations. We summarize the main results of the analysis in Table 1.
Thanks to the analysis tools presented in this paper we were able to prove that the simultaneous read-out of light and heat allows to reach a background of 3.6 +1.9 −1.4 ×10 −3 counts/(keV kg year). The achievement of this background level, the lowest among detectors based on cryogenic calorimeters, sets a key milestone for next generation experiments.