1 Introduction

Two-neutrino (\(2\nu \beta \beta \)) [1] and neutrinoless (\(0\nu \beta \beta \))  [2] double beta decay are among the rarest decay processes studied. While the former is allowed by the Standard Model and has been experimentally detected in a number of isotopes [3], the latter has never been observed; its discovery would imply that lepton number is not conserved and that neutrinos are in fact Majorana particles [4, 5].

The CUORE experiment (Cryogenic Underground Observatory for Rare Events) [6,7,8,9], which is currently running at Laboratori Nazionali del Gran Sasso (LNGS), is designed to perform a high-sensitivity search for \(0\nu \beta \beta \) decay of \(^{130}\hbox {Te}\) to the ground state of \(^{130}\hbox {Xe}\)[10]. The active isotope is contained in \(\hbox {TeO}_{2}\) crystals, which are operated as thermal detectors in a cryostat capable of reaching temperatures below 10 mK. At this temperature, the crystal heat capacity becomes very small and consequently a release of energy within a crystal results in a detectable increase of its temperature. The sought-after experimental signature of \(0\nu \beta \beta \) decay is a monochromatic peak in the summed energy spectrum of the final state electrons at 2527.518 ± 0.013 keV [11,12,13], which is the transition energy of the decay. To maximize the sensitivity of the search, the radioactive background at the transition energy must be kept as low as possible.

The first tower assembled in the CUORE assembly line was operated as a standalone experiment, named CUORE-0 [14], from 2013 to 2015. CUORE-0 was designed to validate several key aspects of CUORE, including detector construction, data acquisition and the analysis framework. In addition to this, CUORE-0 provided a sensitive probe of several rare decays, including \(0\nu \beta \beta \)[15, 16] and \(2\nu \beta \beta \)[17] decay of \({}^{130}\hbox {Te}\) to the ground state of \({}^{130}\hbox {Xe}\), and the \(\beta ^{+}EC\) decay of \({}^{120}\hbox {Te}\)[18]. In this work we focus on a search for double-beta decay of \({}^{130}\hbox {}\) to the first \(0^+\) excited state of \({}^{130}\hbox {Xe}\) (\(\beta \beta _{0^+_1}\)) with CUORE-0. As shown in Fig. 1 this decay emits two electrons, which share a maximum energy of 734 keV, followed by a gamma cascade to the ground state of \({}^{130}\hbox {Xe}\). The most probable de-excitation pattern, which has a branching ratio of 86.0%, involves the emission of two gamma rays with energies of 1257.5 (3) keV and 536.1 keV. Two more patterns are also possible, namely the emission of three gamma rays with energies of 536.1 keV, 586.0 keV and 671.3 keV (branching ratio of 12.2%) and the emission of two gamma rays with energies of 671.3 keV and 1122.1 keV (branching ratio of 1.8%). These gamma lines result in multi-detector coincidence signatures which we exploit in our analysis to achieve very powerful background rejection. As is the case for the decay to the ground state, the summed energy spectrum of the emitted electrons is distinctly different in the \(2\nu \beta \beta _{0^+_1}\) vs. the \(0\nu \beta \beta _{0^+_1}\) case. The former is a continuous spectrum (0–734 keV), whereas the latter is a monochromatic peak centered at 734 keV.

Prior to the current work, the most stringent constraints on these decays came from the CUORICINO  [19] experiment, a predecessor to CUORE-0, which reported the following limits on the decay half-lives [20]:

$$\begin{aligned}&T_{0^+_1}^{0\nu }>9.4\cdot 10^{23}\,\hbox {yr}\,, \\&T_{0^+_1}^{2\nu }>1.3\cdot 10^{23}\,\hbox {yr}. \end{aligned}$$

2 Experiment

The CUORE-0 tower, just like all the CUORE towers, contains 52 \(^{\text {nat}}\hbox {TeO}_2\) crystals (i.e., \(^{130}\hbox {Te}\) is present at its natural abundance of 34.2% [21]). The crystals are arranged in a copper frame into 13 floors, with each floor containing four crystals [14]. The mass of each crystal is 750 g, for a total detector mass of \(\hbox {TeO}_2\) of 39 kg or 10.8 kg of \(^{130}\hbox {Te}\).

The crystals are operated as thermal detectors at a working temperature of \(\sim 10\,\hbox {mK}\) to minimize the heat capacity. In CUORE-0, this temperature is obtained using the same dilution refrigerator infrastructure as CUORICINO  [22]. The detector thermal link to the fridge is provided by the polytetrafluoroethylene supports which hold each crystal in the copper frames and by the golden wires used to carry the electrical signal.

The temperature of each crystal is continuously sensed by a neutron transmutation doped (NTD) thermistor [23] glued to the crystal surface. The thermistor converts the thermal signal to a voltage output which is digitized at an acquisition rate of 125 samples/s. Signal pulses corresponding to thermal events are identified through a software trigger, and a 5 s long window – 1 s before and 4 s after the trigger – is selected for further analysis. The initial second (pre-trigger) is used to establish the baseline temperature of the crystal just prior to the event. The pulse amplitude is used to determine the deposited energy. In order to monitor and correct for changes in detector gain due to temperature drifts, a silicon resistor (heater) [24] is coupled to each crystal and is used to generate reference thermal signals every 300 s [25, 26]. A time-coincidence analysis can be performed to search for events that involve multiple crystals simultaneously. To account for the time response of the detector we use a coincidence window of \(\pm 5\,\hbox {ms}\). As the measured event rate is approximately 1 mHz/crystal, the probability of accidental (i.e., causally unrelated) coincidences is extremely small (\(\simeq 10^{-5}\)).

In order to reduce background due to environmental radioactivity, the tower and cryogenic infrastructure are surrounded by several layers of shielding, including an internal low-background Roman lead layer and an external anti-radon box. The details of the CUORE-0 detector design, operation and performance are described in [14, 16]. The dilution refrigerator, shielding, and other cryostat components are those from the CUORICINO experiment [22, 27].

Fig. 1
figure 1

Decay scheme of \(^{130}\hbox {Te}\) showing the energy levels and the branching ratios for the \(\gamma \) rays [10]

3 Analysis

In this work we exploit the multi-detector coincidence patterns expected to accompany \(\beta \beta _{0^+_1}\) decays to maximize our sensitivity. The emitted electrons and gamma rays can interact in multiple crystals, producing a variety of experimental scenarios or signatures depending on the number of detectors involved in the process (event multiplicity). We first identify all the possible signatures that can be detected and rank them by their expected sensitivity to \(0\nu \beta \beta _{0^+_1}\) or \(2\nu \beta \beta _{0^+_1}\) decays. Our determination of the most significant signatures to be used in the final analysis makes use of both real and simulated data.

3.1 Signature identification

To simplify the analysis we restrict ourselves to signatures in which the energy deposition from electrons is limited to the crystal where the decay took place. We further require that each individual de-excitation gamma is completely absorbed in a single crystal, thus discarding events where these gamma rays scatter but subsequently escape the detector. Scenarios where gamma rays are absorbed in the same crystal where they originated from are considered valid. Since the maximum number of gamma rays emitted in the decay is three (Fig.  1), these choices imply that the maximum number of crystals involved in an event is four (one crystal contains the two electrons). We considered adding a constraint on the proximity of the involved crystals but found that it does not increase our sensitivity significantly.

For every possible signature satisfying these conditions, both for the \(2\nu \beta \beta _{0^+_1}\) and the \(0\nu \beta \beta _{0^+_1}\) decay modes, we identify a single monochromatic line in one of the up-to-four crystals which can be used to estimate the decay rate through the fit described in Sect. 3.5. For this reason, at least two crystals must be involved in \(2\nu \beta \beta _{0^+_1}\) decay signatures, so that at least one of them records a monochromatic peak – when only one crystal is involved the energy deposited is a continuous distribution spread over a 734 keV-wide range. The analysis therefore ultimately involves searching for the peak associated with each signature. Taking as an example the most likely de-excitation pattern (two gammas with energies of 1257.5 keV and 536.1 keV), the allowed scenarios resulting from this set of rules are reported in Table 1.

Table 1 Allowed scenarios for the most likely de-excitation pattern (two gammas with energies of 1257.5 keV and 536.1 keV), resulting from the rules defined in Sect. 3.1. The decay always takes place in detector A. \(\beta \beta \), which represents the energy deposited by electrons, corresponds to 734 keV for the \(0\nu \beta \beta _{0^+_1}\) decay and to any energy in the \(0\div 734\) keV range for the \(2\nu \beta \beta _{0^+_1}\) decay. An energy deposition of 1793.5 keV is recorded when both gammas deposit their energy in the same crystal. Scenarios where only one crystal is involved (Multiplicity 1) are valid only for \(0\nu \beta \beta _{0^+_1}\)

3.2 Sensitivity evaluation

Considering all the possible gamma cascade patterns and the aforementioned constraints, a total of 57 signatures remain. However, only a few of them produce a sizable contribution to the half life sensitivity to double beta decay. In the case of a peak search in the presence of nonzero background the half-life sensitivity has the following dependence on the experimental parameters [28]:

$$\begin{aligned} T^{1/2}_{0^+_1}=k\epsilon \sqrt{\frac{M\cdot t}{b\cdot \varDelta E}}, \end{aligned}$$
(1)

where \(\epsilon \) is the total detection efficiency, \(M\cdot t\) is the exposure, b is the background rate per unit energy (background index) in the energy region of interest (ROI), \(\varDelta E\) is the energy resolution near the ROI and the constant k incorporates the isotopic abundance of \({}^{130}\hbox {Te}\), the molecular mass of the source material and the required significance level. While the exposure (\(35.2\,\hbox {kg}\cdot \hbox {y}\) of \(\hbox {TeO}_{2}\)) and k are always the same, the other three parameters vary significantly depending on the signature and associated peak under consideration.

We measure the energy resolution and its energy dependence directly from CUORE-0 data. We find the resolution exhibits a linear energy dependence which we parametrize as \(\varDelta E (E)=\varDelta E(2615)\times (p_0+p_1\times E)\) [16]; here \(\varDelta E(2615)=4.9\) keV is the average FWHM of the 2615 keV line obtained during the calibration runs, \(p_0=0.49\pm 0.04\) and \(p_1=(2.22\pm 0.15)\times 10^{-4}\). To avoid biasing our ranking procedure, the background index for each signature is estimated from the CUORE-0 background model described in [17] rather than using the CUORE-0 data directly. The model involves a full reconstruction of the measured spectra with highly detailed Monte Carlo simulations based on the Geant4 package [29]. The selection conditions (cuts) associated with each analysis signature are applied to the simulated data and the background index is evaluated in the relevant energy range. Finally, the efficiency term accounts for the probability that a \(\beta \beta _{0^+_1}\) event is triggered, produces the multi-detector coincidence signature in question, and the peak is properly reconstructed at the expected amplitude. The quantities involved in the efficiency calculation are described in detail in Sect. 3.4. Some \(2\nu \beta \beta _{0^+_1}\) scenarios have similar energy deposition patterns, requiring the same event multiplicity and energy of the associated monochromatic peak, but different energy ranges in the crystal where the decay took place. If there is an overlap between these ranges, we apply an additional procedure to properly evaluate the detection efficiency. For example scenarios \(\#6\) and \(\#8\) from Table 1 both search for a time coincidence between a peak at 1257.5 keV and an event in a wide energy range (\(0\div 734\) keV for \(\#6\), \(536.1\div 1270.1\) keV for \(\#8\)). The range between 536.1 and 734 keV is included in both scenarios. If \(P_6\) and \(P_8\) are the probabilities of detecting a decay from scenario \(\#6\) and \(\#8\) respectively, the probability of detecting a decay in the overlapping energy region is \(P_6+P_8\). Any time we observe a situation like this, we modify the original scenarios to exclude the overlapping energy region, which we use instead to create an additional scenario. Referring again to the previous example, a first scenario will require events in the \(0\div 536.1\) keV range, a second in the \(536.1\div 734\) keV range and a third in the \(734\div 1270.1\) range. For this reason, some of the scenarios we considered for our final \(2\nu \beta \beta _{0^+_1}\) analysis do not cover the whole 734 keV-wide range.

We define our total sensitivity to \(0\nu \beta \beta _{0^+_1}/2\nu \beta \beta _{0^+_1}\) decay as the sum in quadrature of the sensitivities given by each signature. We consider a signature to be relevant if its contribution to the total sensitivity to the process exceeds \(1\%\). Only five such scenarios are identified for both the \(0\nu \beta \beta _{0^+_1}\) and the \(2\nu \beta \beta _{0^+_1}\) channels, and they are ranked by their sensitivity in Tables 2 and 3 respectively. The selected scenarios cover \(\sim 97\)% and \(\sim 99\%\) of the total sensitivity for the \(2\nu \beta \beta _{0^+_1}\) and \(0\nu \beta \beta _{0^+_1}\) cases respectively.

Table 2 The five most relevant scenarios that contribute to the total sensitivity to the neutrinoless decay channel by more than 1%. b is the background index (units: counts/keV/kg/y) and \(\epsilon \) the total detection efficiency. The electrons are assumed to be always fully absorbed in the detector ’A’. The numbers within parentheses show how the energy is split between electrons (734 keV) and gammas in detector ’A’. The final fit is performed at the highest energy monochromatic peak measured in each signature, marked here in bold. The observed spectra for each scenario are shown in Fig. 3. Scenario numbers are unrelated to those in Table 1
Table 3 The five most relevant scenarios that contribute to the total sensitivity to the two neutrino decay channel by more than 1%. b is the background index (units: counts/keV/kg/y) and \(\epsilon \) the total detection efficiency. The electrons are assumed to be always fully absorbed in the detector ’A’. The final fit is performed at the highest energy monochromatic peak measured in each signature, marked here in bold. The energy ranges in detector ’A’ do not always cover the full 734 keV-wide range to properly evaluate the efficiency (see Sect. 3.2). The observed spectra for each scenario are shown in Fig. 2. Scenario numbers are unrelated to those in Table 1

3.3 Data selection

We first remove time periods where the data quality are poor; the effect of this is accounted for in the exposure. We next remove events that are either poorly reconstructed by our analysis or are non-signal-like using the pulse shape methods described in [16]. We then impose cuts based on the deposited energy and on the event multiplicity (i.e., the number of crystals involved in the event). For the \(0\nu \beta \beta _{0^+_1}\) case (Table 2), the energy and multiplicity rules a candidate event must pass are listed below.

  • In scenario 1, events must involve exactly three hits in the same coincidence time window. One of the three crystals must contain a signal with energy \(E_{1}\) in the range \((734\pm 5\sigma _{734})\) keV. The notation \(\sigma _{734}\) indicates the energy resolution at 734 keV, which is estimated from the resolution function reported in Sect. 3.4. Another of the three crystals must have energy, \(E_{2}\), in the \((536\pm 5\sigma _{536})\) keV range. No requirement is imposed on the energy deposited in the third crystal.

  • In scenarios 2–5, events must have exactly two crystal hits. No requirement is set on the individual measured energies of the hits (\(E_{1}\) and \(E_{2}\)), but rather on their sum, \(E_{tot}=E_{1} + E_{2}\). Labeling as \(E_{A}\) and \(E_{B}\) the energies indicated in Table 2 that are expected to be deposited in detector A and B respectively, and as \(E_{AB}\) their sum (\(E_{A}+E_{B}\)), then \(E_{tot}\) must lie in the range \((E_{AB}\pm 5\sigma _{AB})\), where \(\sigma _{AB}=\sqrt{\sigma _{A}^2+\sigma _{B}^2}\).

We apply a similar logic to the \(2\nu \beta \beta _{0^+_1}\) decay (Table 3), but the range defined by the continuous electron spectrum replaces the \(\pm 5\sigma \) requirement.

  • In scenario 1 events must involve exactly three hits. The energy of the first hit, \(E_1\), is bound in the range \(0<E_{1}<734\); the energy of the second hit, \(E_2\), must be in the range \((536\pm 5\sigma _{536})\). No requirement is set for the hit on the third crystal.

  • For scenarios 2–5, two hits are required in the same coincidence window. Again, we don’t set a requirement on the individual hit energies, but on their sum, \(E_{tot}=E_{1} + E_{2}\). The continuous electron spectrum, contained in detector A (Table 3), defines two energy limits, \(E_{A}^{min}\) and \(E_{A}^{max}\). We require that \(E_{tot}\) satisfy the condition \(E_{A}^{min}+E_{B}< E_{tot} < E_{A}^{max} + E_{B}\).

3.4 Efficiency evaluation

The detection efficiency is the probability that a \(\beta \beta _{0^+_1}\) event is triggered and properly reconstructed, including the production of the required multi-detector coincidence signature. The detection efficiency is a product of several factors:

  • the probability that an event is triggered, which we estimate from the fraction of heater-induced events that are triggered and correctly reconstructed at the expected amplitude;

  • the probability to include only physical events whose pulse shape does not differ from the average behavior (pulse shape efficiency);

  • the probability to correctly measure the number of crystals involved in an event (i.e., the event multiplicity);

  • the fraction of \(\beta \beta _{0^+_1}\) events that deposit energy according to a particular scenario.

The trigger and pulse shape efficiencies are derived in the same way described in [16]. The multiplicity term is obtained differently depending on the number of crystals involved. The probability that an event involving a single crystal (multiplicity 1) is actually recorded as one is calculated using the \({}^{40}\hbox {K}\) line at 1461 keV [16]: since it’s the only \(\gamma \) line emitted in the decay, the only way for it to be measured in a multiplicity \(>1\) event is by chance. For events with higher multiplicity we don’t have a way to evaluate the efficiency directly, so we resort to a statistical derivation. We take into consideration the two main effects that can alter the multiplicity of one event: accidental coincidences and pile-up.

  • Accidental coincidences are events that happen simultaneously by chance and not due to causally correlated signals. This leads to an artificially increased multiplicity.

  • Pile-up refers to the situation where two (or more) signals happen randomly on the same channel within the same 5 seconds-long signal window. Such events are marked as pile-up and ignored by the coincidence calculation. This artificially reduces the multiplicity; if, for example, two signals are coincident but one of them is removed due to pile-up, the second signal is recorded as a multiplicity one event.

The magnitude of these effects can be estimated with Poisson statistics, considering the measured average event rate of 1 mHz/channel, a 5 s long signal window and a 10 ms coincidence window. The resulting efficiencies for events with multiplicities 1 through 4 are listed in Table 4. While multiplicity 4 scenarios are possible, none of them gives a significant contribution to sensitivity.

The final term in the efficiency, related to \(\beta \beta _{0^+_1}\) decay itself, is evaluated using dedicated Monte Carlo simulations. We use the same simulation software employed for the CUORE-0 background model to reproduce the effects of all decay patterns and compute the fraction of fully-contained \(\beta \beta _{0^+_1}\) events for each of the 57 scenarios. Due to the high statistics of these simulations, the error associated to this efficiency term is extremely small (\(<0.1\%\)).

Table 4 Signal detection efficiency terms. Trigger and pulse shape efficiencies are common to all scenarios. Multiplicity 1 efficiency is calculated from \({}^{40}\hbox {K}\), while multiplicities 2–4 come from a statistical calculation. The same statistical calculation applied to multiplicity 1 events yields the same result (\(99.6\%\))

3.5 Fitting technique

The final step of the analysis procedure is to obtain for each scenario the energy spectrum selected for the fit. We choose to fit the spectrum of the crystal that records the monochromatic peak with the highest energy as, due to the shape of our observed spectra, it usually has the lowest background. To simplify the fits, ranges are chosen to exclude any peak from other \(\gamma \) lines. This selection is based on the CUORE-0 background model[17] rather than the data. The final spectra for the \(0\nu \beta \beta _{0^+_1}\) decay search are shown in Fig. 2, while those for \(2\nu \beta \beta _{0^+_1}\) decay are shown in Fig. 3. The strong background reduction achieved with energy-related cuts and the excellent agreement between real data and the background model are evident in these figures.

Fig. 2
figure 2

CUORE-0 spectra for the five signatures selected for the \(0\nu \beta \beta _{0^+_1}\) decay. The blue histogram shows the data without any energy-related cut; the reconstruction of the background model is shown in red; in black are data with energy cuts. The non-shaded area represents the fit range, while the black arrow points to the location of the expected \(\beta \beta _{0^+_1}\) peak

Fig. 3
figure 3

CUORE-0 spectra for the five signatures selected for the \(2\nu \beta \beta _{0^+_1}\) decay. The blue histogram shows the data without any energy-related cut; the reconstruction of the background model is shown in red; in black are data with energy cuts. The non-shaded area represents the fit range, while the black arrow points to the location of the expected \(\beta \beta _{0^+_1}\) peak

The fit function for each signature is the following:

$$\begin{aligned} B_{const} + B_{lin}\cdot E + \dfrac{\epsilon \cdot t\cdot \varGamma _{\beta \beta }^{0^+_1} }{\sqrt{2\pi \sigma ^2}}\cdot G(E_{\beta \beta }^{0^+_1},\sigma ), \end{aligned}$$
(2)

where \(B_{const}\) and \(B_{lin}\) are parameters describing a linear background, \(\epsilon \) is the detection efficiency, t is the live time, \(\varGamma _{\beta \beta }^{0^+_1}\) is the \(\beta \beta _{0^+_1}\) decay rate, and G is a gaussian function centered at \(E_{\beta \beta }^{0^+_1}\) and with a resolution of \(\sigma \). The expected values of \(E_{\beta \beta }^{0^+_1}\) for each signature are indicated in bold in Tables 2 and 3. We perform a simultaneous unbinned extended maximum likelihood (UEML) fit on the five signatures belonging to each \(\beta \beta _{0^+_1}\) decay, using the RooFit fitting package [30]. The fit parameters are constrained as follows:

  • the decay rate \(\varGamma _{\beta \beta }^{0^+_1}\) is a common parameter for all five signatures;

  • the two background components are independent for each signature;

  • the detection efficiency is fixed at the value reported in Tables 2 and 3 for each signature;

  • the exposure is fixed at 35.2 kg\(\cdot \)y;

  • the energy resolution is fixed at the value determined from the \(\varDelta E(E)\) curve reported in Sect. 3.4;

  • the peak position (\(E_{\beta \beta }^{0^+_1}\)) is fixed at the expected value reported in Tables 2 and 3 .

Each fit has 11 free parameters: two background parameters for each of the five signatures (for a total of 10) and the decay rate \(\varGamma _{\beta \beta }^{0^+_1}\), which is common for all signatures.

4 Results

We find no evidence of a \(\beta \beta _{0^+}\) signal, either for the neutrinoless or the two neutrino decay mode. We set a \(90\%\) confidence upper limit on the decay rates for the two processes: \(\varGamma _{0^+}^{0\nu } < 8.8 \cdot 10^{-25} y^{-1}\), \(\varGamma _{0^+}^{2\nu } < 2.8 \cdot 10^{-24} y^{-1}\). In turn, these correspond to the following lower limits for the half lives:

$$\begin{aligned}&T_{0^+_1}^{0\nu }> 7.9 \cdot 10^{23}\hbox { y}, 90\%\hbox { CL},\\&T_{0^+_1}^{2\nu } > 2.4 \cdot 10^{23}\hbox { y}, 90\%\hbox { CL} \end{aligned}$$

We estimate the systematic uncertainties with a procedure identical to that applied in the analysis for \(0\nu \beta \beta \) decay to the ground state [15, 16]. We consider two components for the systematic uncertainty: one factor which scales with the decay rate (\(\sigma _{scaling}\)), and one which is independent of the decay rate (\(\sigma _{add}\)). We generate a large number of simulated spectra with a distribution taken from the best fit of each signature, but with the value of a single nuisance parameter modified by \(1\sigma \); we then fit the simulated spectra with the unmodified parameter. To probe the value of \(\sigma _{scaling}\), we repeat the analysis including in the simulated spectra a fake signal of variable strength. We regress the resulting best-fit decay rates against the simulated values to determine \(\sigma _{add}\) and \(\sigma _{scaling}\). This procedure is applied separately to get the systematic contributions from the uncertainty on efficiency, energy resolution and peak position. We also run the simulation without changing any parameter, to check for a possible fit bias. The resulting values for \(\sigma _{add}\) and \(\sigma _{scaling}\) are reported in Table 5 for \(0\nu \beta \beta _{0^+_1}\) and in Table 6 for \(2\nu \beta \beta _{0^+_1}\). In both cases, the dominant effect is a small negative bias.

Table 5 Summary of systematic uncertainties on \(\varGamma _{0^+_1}^{0\nu }\)
Table 6 Summary of systematic uncertainties on \(\varGamma _{0^+_1}^{2\nu }\)

We combine the CUORE-0 likelihood curves with those from CUORICINO [20] (Fig. 4 for \(0\nu \beta \beta _{0^+_1}\) and Fig. 5 for \(2\nu \beta \beta _{0^+_1}\)). We set limits for the decay rates taking into account both the CUORE-0 systematic effects and the combination with the CUORICINO results: \(\varGamma _{0^+_1}^{0\nu } < 4.80 \cdot 10^{-25} y^{-1}\), \(\varGamma _{0^+_1}^{2\nu } < 2.73 \cdot 10^{-24} y^{-1}\). These yield the following limits on the half lives:

$$\begin{aligned}&T_{0^+_1}^{0\nu }> 1.4 \cdot 10^{24}\hbox { y}, 90\%\hbox { CL},\\&T_{0^+_1}^{2\nu } > 2.5 \cdot 10^{23}\hbox { y}, 90\%\hbox { CL} \end{aligned}$$

Thanks to the improved background and analysis techniques, we achieve similar reach as CUORICINO with less than half of the exposure. The lower limits for the half lives of both the \(0\nu \beta \beta _{0^+_1}\) and \(2\nu \beta \beta _{0^+_1}\) decays obtained by the combination of the results from CUORE-0 and CUORICINO are the best currently available. CUORE will achieve even higher sensitivity, thanks to the improved background and to the powerful coincidence analysis made possible by its closely-packed 988 crystals.

Fig. 4
figure 4

Negative log likelihood (NLL) from CUORE-0 and the combination with CUORICINO for \(0\nu \beta \beta _{0^+_1}\) decay

Fig. 5
figure 5

Negative log likelihood (NLL) from CUORE-0 and the combination with CUORICINO for \(2\nu \beta \beta _{0^+_1}\) decay