Search of the neutrino-less double beta decay of $^{82}$Se into the excited states of $^{82}$Kr with CUPID-0

The CUPID0 experiment searches for double beta decay using cryogenic calorimeters with double (heat and light) read-out. The detector, consisting of 24 ZnSe crystals 95$\%$ enriched in $^{82}$Se and 2 natural ZnSe crystals, started data-taking in 2017 at Laboratori Nazionali del Gran Sasso. We present the search for the neutrino-less double beta decay of $^{82}$Se into the 0$_1^+$, 2$_1^+$ and 2$_2^+$ excited states of $^{82}$Kr with an exposure of 5.74 kg$\cdot$yr (2.24$\times$10$^{25}$ emitters$\cdot$yr). We found no evidence of the decays and set the most stringent limits on the widths of these processes: $\Gamma$($^{82}$Se $\rightarrow ^{82}$Kr$_{0_1^+}$)$<$8.55$\times$10$^{-24}$ yr$^{-1}$, $\Gamma$($^{82}$Se $\rightarrow ^{82}$Kr$_{2_1^+}$)$<6.25 \times10^{-24}$ yr$^{-1}$, $\Gamma$($^{82}$Se $\rightarrow ^{82}$Kr$_{2_2^+}$)$<$8.25$\times$10$^{-24}$ yr$^{-1}$ (90$\%$ credible interval


Introduction
The double beta decay is a transition among isobaric isotopes (A, Z) → (A, Z + 2) + 2e − + 2ν e . Despite being among the rarest processes in Nature, it was observed for eleven nuclei with typical half-lives of 10 18 -10 24 years [1]. In 1937, Furry hypothesized that double beta decay could occur also without the emission of neutrinos, in the form (A, Z) → (A, Z + 2) + 2e − [2]. This process, called neutrino-less double beta decay (0νDBD), is forbidden by the Standard Model of Particle Physics as it would violate the difference between the total number of baryons and leptons (B-L) [3,4]. Furthermore, 0νDBD is considered a golden channel to probe a fundamental property of neutrinos, i.e. their nature. This transition, indeed, can occur only if (in contrast to all the other known fermions) neutrinos coincide with their own anti-particles, as predicted by Majorana [5]. Finally, the measurement of the 0νDBD halflife would help in understanding the absolute mass scale of neutrinos, that today is one of the missing elements in the puzzle of Particle Physics [6].
CUPID-0 is the first medium-scale 0νDBD cryogenic experiment exploiting the dual read-out of heat and light for background suppression [7]. The detectors are operated as calorimeters [8]: each crystal acts as energy absorber converting energy deposits ∆E into temperature variations ∆T . The temperature variation ∆T is determined by the crystal thermal capacitance C (∆T ∝ ∆E/C) and, for a single-particle energy deposition of ∼1 MeV, it is possible to observe sizeable signals (hundreds of µK) only if C is of the order of 10 −9 -10 −10 J/K. Since in dielectric and diamagnetic crystals C ∝ T 3 according to the Debye law, such thermal capacitances require the crystals to be cooled down to about 10 mK. The temperature variations are converted into readable voltage signals using a Neutron Transmutation Doped (NTD) Ge thermistor [9] glued to the crystal. The resistance of this device shows a strong dependency on the temperature: R(T) = R 0 exp(T 0 /T) γ with R 0 , T 0 and γ of about 2 Ω, 4.2 K, and 0.5 respectively. Thus, biasing the thermistor with a small current allows to convert temperature variations in electrical signals with a temperature sensitivity of hundreds of mK per MeV (or hundred of µV/MeV).
The technological effort of operating tens of massive crystals at cryogenic temperatures is motivated by the advantages that this technique offers in terms of energy resolution, efficiency, and versatility in the choice of the emitter. The CUORE experiment [10,11,12] is successfully operating 988 TeO 2 calorimeters for the study of the 130 Te 0νDBD, proving the feasibility of a tonnescale experiment based on this technology. According to the CUORE background model, the dominant contribution to the region of interest stems from α particles emitted by the materials in the proximity of the detector [13].
The primary goal of the CUPID-0 experiment is proving that the dual read-out heat/light allows to reject the α interactions, reducing the background in the region of interest for 0νDBD by an order of magnitude. This would allow to set an important milestone for next-generation projects aiming at working in an almost background-free environment to increase the discovery potential [14,15,16,17]. For this purpose, each calorimeter is coupled to a light detector that enables particle identification exploiting the different light yield of different particles.
The CUPID-0 detector has been taking data since the end of March 2017 in the underground Laboratori Nazionali del Gran Sasso (LNGS) in Italy. The first data release demonstrated the potential of this technology: thanks to the strong background suppression, CUPID-0 set the most stringent limit on the half-life of the 82 Se decay to the ground state of 82 Kr [18], improving by an order of magnitude the previous limits despite the small exposure (1.83 kg·y compared to the 4.90 kg·y collected by NEMO-3 [19,20]).
In this work we search for the 82 Se decay to the 0 + 1 , 2 + 1 , 2 + 2 excited levels of its daughter, 82 Kr. These transitions were already studied with an exposure of 3.64×10 24 emitters·yr using a high purity Germanium detector operated underground at LNGS [21]. With an efficiency ranging from 0.3 to 3.2% (depending on the chosen signature), and a background of 9.6±0.5 c/keV/y, mainly ascribed to the multi-scatter of high energy γ's produced by contaminations of the set-up, it was possible to set the following 90% C.I. limits on the widths of the decays: In contrast to the measurements made with γ spectroscopy, in CUPID-0 we can distinguish the decay with two neutrinos emission from neutrino-less double beta decay. In this paper we present the results obtained with first data of CUPID-0 in the search of the neutrino-less double beta decay mode.

The CUPID-0 Detector
After an extensive R&D on scintillating crystals based on different 0νDBD emitters [22,23,24,25,26,27,28,29,30,31,32,33,34,35] the CUPID-0 collaboration decided to focus on 82 Se. The relatively long half-life of the 2-neutrino decay-mode (T 2ν 1/2 = [9.39 ± 0.17 (stat) ± 0.58 (syst)] × 10 19 yr [20]) prevents pile-up events in the region of interest despite the slow time response of cryogenic calorimeters (∼ms). The high Q-value of this isotope (2997.9±0.3 keV [36]) allows to reduce the background due to the environmental radioactivity, that drops above the 2615 keV line of 208 Tl. The rather low natural isotopic abundance of 82 Se (8.82% [37]) was increased via isotopic enrichment to 96.3% by the URE-NCO Stable Isotopes unit (Almelo, Netherlands). The obtained 82 Se was used to synthesize the ZnSe powder, that was later purified and doped using ZnSe(Al) with natural isotopic composition of Zn and Se, in order to enhance the light output.
The ZnSe powder was used to grow 24 cylindrical ZnSe crystals 95% enriched in 82 Se [38]. The detector includes also two natural ZnSe crystals, not considered for this work. Since we optimized the crystal shape in order to prevent losses of enriched material, and since we had to reduce the mass of some crystals to discard inclusions and imperfections, the ZnSe detectors feature slightly different size and mass. The total mass of the 24 Zn 82 Se crystals amounts to 9.65 kg, but two of them are not used for the analysis presented in this paper because of their poor bolometric performance. Thus, the mass considered for the analysis is 8.74 kg (3.41×10 25 nuclei of 82 Se).
The light produced by the ZnSe scintillation (a few % of the total energy released as heat) escapes the crystal and is recorded using two light detectors. The fraction of energy converted in form of light depends on the nature of the interacting particle, enabling particle identification and, ultimately, the rejection of the α background [39]. In CUPID-0 also the light detectors are operated as cryogenic calorimeters, meaning that they convert the impinging photons in temperature variations using NTD Ge thermistors [40]. Nevertheless, in this paper we do not describe the details of the light detectors, as the analysis of coincidences among detectors already provides a sufficient background suppression.
The ZnSe crystals, surrounded by a VIKUITI multilayer reflecting foil produced by 3M, and interleaved by light detectors, are assembled in 5 towers using PTFE pieces and a mechanical structure made of NOSV copper (produced by Aurubis AG). Each detector was equipped with a Si Joule resistor that periodically injects a reference pulse to correct thermal drifts [41,42].
More details concerning the detector construction and operation, the 3 He/ 4 He dilution refrigerator, the electronics and data-acquisition can be found in Ref. [7].

Expected Signatures in CUPID-0
The decay scheme of 82 Se is shown in Fig. 1. If 82 Se decays to the ground state of 82 Kr, the two emitted electrons share the entire Q-value of the transi-tion. From Monte Carlo simulations, the probability for the two electrons to release the full energy in the crystal where they were emitted, thus producing a peak at the Q-value of the 0νDBD, is 81.0 ± 0.2%.
The scenario becomes slightly more complicated if the decay occurs to an excited level of 82 Kr. In this case, the energy of the decay is split between the electrons and the γ rays emitted during the de-excitation of 82 Kr. The containment efficiency of the two electrons does not vary significantly, but the probability that a coincident γ ray with energy ranging from 698 keV (γ 2 ) to 1475 keV (γ 3 ) releases its full energy in the crystal is very low, leading to an important decrease of the detection efficiency. The γ rays produced in the deexcitation, indeed, can be fully absorbed in the crystal, or escape the crystal and be absorbed in another one, or escape the crystal and scatter in another crystal, or completely escape detection. Depending on the scenario, we expect different signatures. In addition, more decay schemes can result in the same signature, further complicating the analysis. The redundancy of states, as well as the different detection efficiency of the processes, impose a down-selection of the decay schemes.
First, we exclude from the analysis the events in which a single crystal triggers, such as those in which the γ rays escape detection. This choice is motivated by the high background produced mainly by the 2-neutrino double beta decay of 82 Se (the reader can find a plot of the physics spectrum obtained imposing that a single ZnSe triggered in Ref. [18]).
Then, we restrict the analysis to events in which only two ZnSe crystals trigger, thus rejecting interactions in three or more crystals. This choice, that excludes for example the scenario in which the electrons and the two γ's interact in three different crystals, is motivated by the very low efficiency of these signatures.
Finally, we discard the remaining signatures with efficiency lower than 0.01% that would not give a substantial contribution to the analysis.
To compute the detection efficiency, we simulate 10 7 decays in the CUPID-0 crystals, accounting also for a linear smearing of the energy resolution as a function of the energy (Sec. 5). We select events in which two detectors trigger and the energy measured by one of them (E coinc ) is compatible with the energy of a deexcitation γ ray. Then, we search for a peak in the other crystal with energy E main equal to the total energy of the two electrons E ββ , or to the sum of the energy of the electrons and another γ ray emitted in the same decay E ββ +E γi . The signatures chosen for the analysis presented in this paper and their detection efficiencies are summarized in Table 1. ; γ i are the γ rays emitted in the de-excitation to the ground state (Fig. 1); the vertical bar separates the particles releasing their full energy (E main ) in the 1 st crystal, and the particles releasing their full energy (E coinc ) in the second crystal.The detection efficiency is determined by a Monte Carlo simulation. Different decay schemes resulting in the same signature (for example, 1, 7, 11) are labelled with the same letter in the last column; the letter B indicates two states with a slightly different energy ββ, that were grouped given the resolution of the detector.
The number of decays N i corresponding to the i th signature can be written as a function of the exposure (ξ = 2.24 × 10 25 emitters·yr), the detection efficiency i (Table 1), the data selection efficiency η (Sec. 5) and the width of the corresponding decay channel (Γ ): Each signature can be modeled with a function describing the detector response to a monochromatic energy deposit (Σ), a flat background component (ρ f lat ) and, if the mean energy is close to the energy of the 40 K line (signatures B, C, D, F), a peaking background accounting for accidental coincidences due to the higher rate of this peak (Σ 40K ): . . .
In Sec. 4 we describe how data are acquired and processed. In Sec. 5 we derive a model for the detector response Σ as a function of the energy and compute the data selection efficiency. Finally, in Sec. 6 we perform the simultaneous fit of the models f A . . . f G to the data to extract the values of Γ 0 + 1 , Γ 2 + 1 and Γ 2 + 2 .

Data Collection and Processing
The temperature variation produced by 1 MeV energy deposit in a ZnSe results in a voltage signal of tens of µV, with typical rise-time of 10 ms and decay-time ranging from 15 to 60 ms, depending on the detector. The voltage signals are amplified and filtered using a Bessel 6 poles anti-aliasing filter with tunable cut-off frequency and gain (more details about the electronics and read-out can be found in Refs. [43,44,45,46,47,48,49,50,51]). The data acquisition system digitizes all the ZnSe channels with a sampling frequency of 1kHz and saves the corresponding data on disk in NTuples based on the ROOT software framework. During the measurement we run a software trigger on the acquired data and save the corresponding timestamps in the NTuples for the off-line analysis. The trigger algorithm is sensitive to the derivative of the waveforms and its configuration parameters are optimized separately for each channel [52,53].
The analysis presented in this work comprises 6 DataSets, each consisting of a collection of physics runs of about 2 days, plus an initial and final calibration with 232 Th sources to monitor the detector stability (see Table 2). In order to include the Q-value of the 0νDBD in the calibration data, we performed a run with a short living 56 Co source. The source emits γ rays up to 3.5 MeV and has been used at the end of the data taking cycle. This calibration is used also to study the energy dependency of the energy resolution (Sec. 5).
The first DataSet was devoted to the detector commissioning and optimization, and for this reason it shows the lowest fraction of live-time. This DataSet was not used for the 0νDBD analysis presented in Ref. [18] because of the poor rejection of the α background due to the variations of the working conditions of the light detectors. Concerning the analysis presented in this paper, we do not expect α particles to contribute to the background. Indeed, the searched processes produce events occurring simultaneously in two crystals. Due to the detector shielding, coinciding alpha particle interactions in multiple detectors are very unlikely and thus are not to be considered a background contribution to the analysis. Thus, we decided to include also the first DataSet to increase the statistics. The total ZnSe collected exposure (enriched crystals only) amounts to 5.74 kg·yr, corresponding to 3.05 kg·yr of 82 Se (2.24×10 25 emitters·yr). These values account for the dead-time due to detector problems (such as earthquakes or major underground activities) and also for the loss of two enriched crystals due to a non-satisfactory bolometric performance.
The collected data are processed off-line using a C++ based analysis framework originally developed by the CUORE-0 collaboration [54,55,56,57]. The continuous data stream is converted into acquisition windows of 4 s, with a pre-trigger window of 1 second to evaluate the detector instantaneous temperature before the pulse occurred.
The data are filtered with a software matched-filter algorithm [58,59] to improve the signal-to-noise ratio. The filtered amplitude is then corrected for gain instabilities using the reference pulses periodically injected through the Si resistors [60,61]. The corrected amplitude is converted into energy using the most intense peaks produced by the 232 Th source between 511 keV and 2615 keV. Each peak is modelled with a Gaussian + exponential background function. The obtained mean values are then converted into energy using a 2 nd order function.
Finally, we compute time coincidences between detectors with a coincidence window of 20 ms, earlier optimized by studying the time-distribution of real coin-cidence events collected during the 232 Th calibrations. This choice was made conservatively to ensure a 100% selection efficiency, at the cost of a possibly larger background due to accidental coincidences. Nevertheless, in the next section we show that, given the low detector rate in physics runs, the number of random coincidences is almost negligible.

Data Analysis
To infer the detector response to a monochromatic energy release Σ, we study the 2615 keV line produced by the decay of 208 Tl. As explained in Ref. [18], the simplest model describing this peak is the sum of two Gaussian functions G with two different σ and mean values: As of today, we do not know the underlying physics behind this bi-Gaussian response. Nevertheless, a multi-Gaussian description of the signal was already observed in other experiments based on cryogenic calorimeters [54,62]. To account for possible time-variations of the detector response, we fit this model to the 2615 keV peak in the sum energy spectrum of all the periodical 232 Th calibrations, obtaining f 1,2 = 0.83±0.03, µ 1 = 2613.88± 0.13 keV, σ 1 = 8.89 ± 0.12 keV, and µ 2 = 2628.37 ± 1.42 keV, σ 2 = 15.42 ± 0.49 keV.
To validate the calibration with 232 Th sources and characterize the detector response over the range of interest, we perform a run with 56 Co sources, emitting many γ rays from 511 to 3451 keV. We energycalibrated the 56 Co spectra with the coefficients derived from the 232 Th calibration, and we fit the model Σ(µ 1 , µ 2 , σ 1 , σ 2 , f 1,2 ) to the most prominent peaks by fixing the value of f 1,2 and the ratios µ 2 /µ 1 and σ 2 /σ 1 to those obtained fitting the 2615 keV line. This choice allows to fix the shape of the fit function and keep as free parameters only the mean value and resolution of the most prominent peak of the bi-gaussian function: Σ(µ 1 , µ 2 , σ 1 , σ 2 , f 1,2 ) → Σ(µ 1 , σ 1 ). In the following we will always use the bi-gaussian model, for simplicity we will drop the subscript 1 and write Σ(µ, σ).
The parameters µ and σ extracted from the fit of the most prominent 56 Co peaks are reported as a function of the energy in Fig. 2.
This study shows that the residuals of the calibration in the region of interest range from -0.5 to 2.5 keV. We check possible differences in the physics runs with respect to calibration runs by fitting the same model to the γ line of 40 K, obtaining µ = 1458.11 ± 0.41 keV, i.e, a residual of 2.7 keV that will be treated as source of systematic error. Fig. 2 shows that the resolution of the peaks scales linearly with energy: σ(E) = σ 0 +aE, with σ 0 = 2.47 ± 0.19 keV and a = (2.4 ± 0.1) × 10 −3 . σ 0 is related to the electronics noise of the detector and is negligible at higher energies, where the loss in energy resolution is dominated by the propagation of phonons in the crystal lattice. The linear dependency on energy is used to derive the correct signal model at the energy of each signature Σ(µ, σ(E)). These models were used to construct the pdfs and finally fit the selected data.
The data used for this analysis are selected by imposing a time-coincidence between two crystals in a 20 ms time-window. Given the low rate in physics runs (∼2 mHz) and the narrow time-window for coincidences, we expect a coincidence rate of 3.2×10 −7 events/s in the range from 0 to 10 MeV. Adding also the requirement on the energy of the events (Table 1) brings the rate of accidental coincidences to a negligible value ranging from <0.08 counts (signatures A, E, G) to 0.5 counts (signature B) integrated over the entire exposure of 5.74 kg·yr. As a consequence, there is no need to exploit the algorithms for background suppression developed for the search of the 0νDBD to the ground state [63]. On the contrary, we apply only a basic cut to the events, selecting windows in which a single pulse is present, to prevent a wrong estimation of the pulse amplitude due to an unpredictable response of the matched filter in presence of more pulses.
The total efficiency comprises the trigger efficiency, the energy reconstruction efficiency, and the efficiency of the quality-cuts applied to the data. The trigger efficiency, computed on the reference pulses injected by the Si resistor, is defined as the ratio of the triggered to injected pulses. The energy reconstruction efficiency is defined as the number of reference pulses reconstructed within ±3σ off the mean energy. Their combined value results in 99.50±0.02%.
The data selection efficiency is computed on the most intense γ peak of the physics spectrum, due to the decay of 65 Zn (T 1/2 ∼224 d, Q-value∼1351.9 keV). The accepted and rejected events are simultaneously fitted with an un-binned extended maximum likelihood fit (with the RooFit analysis framework), resulting in a selection efficiency of 96.0±0.4% (Fig. 3).  Fig. 3 γ peak produced by the decay of 65 Zn. Top: events in which the acquisition window contains a single pulse. Bottom: events in which the acquisition window contains more than one pulse. The spectra are fitted simultaneously with an unbinned extended maximum likelihood fit (RooFit analysis framework) with two components: the function modeling the detector response Σ(µ, σ(E)), and an exponential background.
The analysis was repeated on the γ peak produced by 40 K at 1.46 MeV, obtaining a compatible result (96.0±1.1%).
The combination of the selection efficiency with the trigger and energy reconstruction efficiencies results in a total efficiency η = 95.5 ± 0.4%.

Results
The energy spectra of signatures A, . . . , G are reconstructed starting from the selected events. Following Table 1, for each spectrum we consider a ∼400 keV range centered around E main , as this is a reasonable interval in which the background rate can be considered flat. This spectrum is further selected requiring that the events are in a 20 ms time-coincidence with events in [E coinc − 2σ, E coinc + 2σ]. The values of σ are chosen using the results shown in Fig. 2.
We report in Fig. 4 the energy spectra corresponding to signatures A (<0.08 background counts expected) and B (0.5 background counts expected), to show the two limit cases. Because of the low background rate, the spectra of the other signatures are very similar to those displayed in these plots. The fit accounts also for various sources of systematic errors. We evaluate the uncertainty related to the calibration function assuming an uniform distribution in the interval [µ-∆, µ+∆], with ∆ =2.7 keV, conservatively chosen from the study of the residuals.
The error due to the extrapolation of σ at the energies where we expect the decays, as well as the systematic error induced by the uncertainty on the selection efficiency η, are modelled by weighting the likelihood by a gaussian function with the mean fixed to σ (η) and RMS fixed to the uncertainty on σ (η). We finally integrate the likelihood by a numerical method.
The 90% credible intervals Bayesian upper limit are set using a uniform prior on the values of Γ 0 + (3)

Conclusions
In this paper we presented the first background-free search of the neutrino-less double beta decay of 82 Se to the excited states of 82 Kr with an exposure of 2.24×10 25 emitters·yr. Despite the low efficiency, we were able to set the most competitive upper limits on the decay widths of the 0 + 1 , 2 + 1 and 2 + 2 levels. The detector is still taking data at LNGS with the aim of reaching a ZnSe exposure of 10 kg·yr, that will allow to further improve this result.