Main

The standard model of particle physics is a successful paradigm for the number, properties and interactions of fundamental particles. Nevertheless, the observation of neutrino oscillations indicates the incompleteness of the standard model: they imply non-vanishing neutrino masses, requiring an extension of the standard model, and violate three accidental symmetries connected to the flavour lepton numbers Le, Lμ and Lτ, leaving the difference between the baryon and lepton number, B − L, as the only unprobed quantity. A promising process to experimentally test B − L is neutrinoless double beta (0νββ) decay, in which a nucleus of mass number A and charge Z decays by the emission of only two electrons: (A, Z) → (A, Z + 2) + 2e. We highlight that this process creates two electrons, namely two matter particles4. This decay can be mediated by various non-standard model mechanisms involving Majorana neutrino masses. A minimal extension of the standard model Lagrangian adds heavy Majorana neutrinos that mix with the known neutrinos to produce a set of light Majorana neutrinos, explaining the observed light neutrino masses5 and at the same time providing a mechanism to explain the baryon asymmetry in the universe2. At this time, experimental searches for 0νββ decay are the most sensitive means to corroborate this framework.

The 0νββ decay signature is a peak in the spectrum of summed energy of the two emitted electrons at the mass difference (Qββ) between the parent and daughter nuclei. A worldwide quest is ongoing, involving a range of nuclei such as 76Ge6,7, 136Xe8,9 and 130Te. The latter, in the form of TeO2 cryogenic calorimeters, is used by the Cryogenic Underground Observatory for Rare Events, CUORE10,11.

To fully exploit the potential of TeO2 crystals as cryogenic calorimeters, the CUORE Collaboration designed and built to our knowledge the largest dilution refrigerator ever constructed, capable of cooling approximately 1.5 t of material to a temperature of approximately 10 mK and maintaining it for years with a 90% duty cycle (1 t = 1,000 kg). In this Article, we describe the performance of CUORE over a four-year measurement campaign and the results of a new high-sensitivity 0νββ decay search with over 1 t yr of TeO2 exposure.

The CUORE experiment

CUORE is the culmination of thirty years of 0νββ decay searches with TeO2 cryogenic calorimeters12. 130Te benefits both from a high natural isotopic abundance of approximately 34%13 and a high Qββ of 2,527.5 keV14, placing the 0νββ decay region of interest above most natural γ-emitting radioactive backgrounds. The detector is an array of 988 natTeO2 cubic crystals15 (Fig. 1) of 5 × 5 × 5 cm3 size and ~750 g mass, for a total mass of 742 kg, which corresponds to 206 kg of 130Te. The array is arranged as 19 towers, each comprised of 13 floors containing four crystals. The crystals are operated as cryogenic calorimeters16 at a temperature of approximately 10 mK. To achieve this low-temperature environment, a novel cryogenic infrastructure—the CUORE cryostat—has been realized.

Fig. 1: The CUORE detector.
figure 1

a, Rendering of the six-stage cryostat, with the pulse tubes and dilution unit, the internal low-radioactivity modern and Roman lead shields, and the array of 988 TeO2 crystals (light blue). b, The detector after installation. The plastic ring was used during assembly for radon protection. c, One of the calorimeters instrumented with an NTD Ge thermistor which measures the temperature increase induced by absorbed radiation. The Si heater is used to inject pulses for thermal gain stabilization. The polytetrafluoroethylene (PTFE) supports and the gold wires instrumenting the NTD and the heater provide the thermal link between the crystal and the heat bath, that is, the Cu frames24.

In a cryogenic calorimeter, the energy deposited by impinging radiation in the absorber crystal is turned into heat, resulting in a temperature rise (Extended Data Fig. 1). Each CUORE crystal (Fig. 1c) is instrumented with a neutron-transmutation-doped germanium thermistor (NTD)17 that converts thermal pulses into electric signals and a heater18 to inject reference heat pulses for thermal gain stabilization19. Thermal signals can be induced by electrons emitted in 0νββ decays but also other background radiation, for example, γ and α particles from residual radioactive contaminants or cosmic ray muons.

CUORE is protected by several means against backgrounds that can mimic a 0νββ decay. It is located underground at the Laboratori Nazionali del Gran Sasso (LNGS) of INFN, Italy, under a rock overburden equivalent to approximately 3,600 m of water, which shields from hadronic cosmic rays and reduces the muon flux by six orders of magnitude. Environmental γ backgrounds are suppressed by a 30-cm layer of low-radioactivity lead above the detector (Fig. 1), a 6-cm-thick lateral and bottom shield of 210Pb-depleted ancient lead recovered from a Roman shipwreck20 (Extended Data Fig. 2), and a 25-cm-thick lead shield outside the cryostat. Environmental neutrons are suppressed by a 20-cm layer of polyethylene and a thin layer of boric acid outside the external lead shield. Finally, radioactive contaminants in the crystals and in the adjacent structures are minimized by careful screening of material for radio-purity and use of high-efficiency cleaning procedures and manipulation protocols21.

Cryogenic innovation and performance

Dilution refrigerator technology was originally proposed in the 1950s22 and underwent considerable development in the 1980s driven also by the application of cryogenic calorimeters for single-particle detection23. Gradually, experimental volumes of the order of tens of litres capable of hosting cold masses of up to 60 kg at 10 mK temperature24 were achieved. Ultimately, detectors were limited by the capacity, duty cycle and radio-purity of commercial or near-commercial cryogenic systems. In the context of this history, the CUORE cryostat represents a major advance in cryogenic technology, reaching an experimental volume of approximately 1 m3 and a cold mass of 1.5 t (detectors, holders, shields) at 10 mK, which corresponds to a 20-fold improvement in experimental volume and target mass compared to the previous state of the art at this temperature scale. Prior to CUORE, the ultimate temperature for comparable target masses was in the resonant-mass gravitational antenna community at 65 mK23.

The CUORE detector is hosted in a multistage cryogen-free cryostat25 (Fig. 1), equipped with five pulse tube cryocoolers that avoid pre-cooling with a liquid helium bath, thus enabling a high duty cycle. A custom-designed dilution unit with a double condensing line for redundancy provides more than 4 μW cooling power at 10 mK. The cryostat is uniquely designed to provide the necessary i) cooling power and temperature stability over a time scale of years, ii) very low radioactivity environment, and iii) extremely low-vibration conditions. As shown in Fig. 2a, b, CUORE became operational in 2017, with the initial period mostly devoted to characterization and optimization campaigns. Since 2019, the data-taking has proceeded smoothly with a duty cycle of approximately 90%. Figure 2d shows that the temperature stability achieved is at the level of 0.2% (±3σ range) over a period in excess of one year. Such a stability is important to achieve a uniform response of all detectors over time. The CUORE calorimeters are sensitive to thermal signals and feature an intrinsic thermal fluctuation limit of approximately 0.5 keV, so any process inducing heat dissipation equal to or greater than 0.5 keV degrades the energy resolution. Mechanical vibrations can be transferred to the inner components and produce heat through friction. To minimize the impact of vibrational noise, the calorimeter array is mechanically decoupled from the cryostat by a custom suspension system. Vibrations induced by the pulse tubes at the 1.4-Hz operational frequency and its harmonics are particularly relevant. In CUORE, we actively tune the pulse tube relative phases for vibration cancellation26 (Fig. 3). This solution is transferable to any cryogenic application involving signals in the same bandwidth of the pulse-tube-induced noise.

Fig. 2: Cryogenic performance.
figure 2

a, The exposure accumulated by CUORE (teal), along with the exposure used for this analysis (orange). Parts of 2017 and 2018 were dedicated to maintenance and optimization of the cryogenic set-up. b, Since then, CUORE has been operating stably with a 90% duty cycle (March 2019–October 2020). c, Examples of temperature instabilities induced by external causes. From left to right: blackout (June 2019), earthquake in Albania of magnitude 6.4, 520 km away (November 2019), regular maintenance (July 2020), and insertion of calibration sources (September 2020). d, The temperature stability of CUORE over ~1 yr of continuous operation, shown by a plot of relative temperature fluctuation versus time, and a histogram of the same data. (1 t yr = 1,000 kg yr.).

Fig. 3: Pulse tube phase optimization.
figure 3

a, Frequency spectrum of the noise for a random combination of the pulse tube phases (orange) and after the active phase tuning (teal). For reference, the frequency spectrum of physical signals is also reported. b, Integral of the power spectrum at the pulse tube frequency (1.4 Hz) and its harmonics before and after active noise cancellation.

CUORE now collects sensitive exposure with 984 out of 988 calorimeters, at a rate that is, to our knowledge, unprecedented for cryogenic particle detectors. Below, we describe the data treatment and 0νββ decay search with greater than 1 t yr of TeO2 exposure.

Data analysis and results

CUORE data are subdivided into datasets of 1–2 months of physics data, separated by a few days of calibration data collected with the detector exposed to 232Th and/or 60Co sources.

The voltage across each NTD is amplified, passed through an anti-aliasing filter, and continuously digitized with a 1-kHz sampling frequency27,28. We identify thermal pulses in the data stream and compute the pulse amplitudes, applying optimum filters that maximize the frequency-dependent signal-to-noise ratio29. To monitor and correct for possible drifts of the thermal gain of the detectors we exploit two ‘standard candles’: monoenergetic heater pulses for the calorimeters with functioning and stable heaters (>95% of the total), and events from the 2,615-keV 208Tl calibration line for the remainder. Drift-stabilized pulse amplitudes are converted to energy using the regularly acquired source calibration data30. We blind the 0νββ search via a data-salting procedure that produces an artificial peak at Qββ30. Once the full analysis procedure is finalized and frozen, we reverse the salting.

To simplify the analysis, we eliminate data from periods affected by high noise or failed data processing, which amounts to 5% of the exposure. Furthermore, calorimeters with greater than 19-keV full width at half maximum (FWHM) energy resolution at the 2,615-keV calibration line are discarded, adding 3% loss in exposure. In addition to these so-called base cuts, the following second-level cuts are then applied to suppress single background-like or low-quality events. Monte Carlo simulations show that approximately 88% of 0νββ decay events release their full energy in a single crystal31. Hence, we apply an anti-coincidence cut that excludes events depositing energy in more than one crystal. Finally, we use pulse shape discrimination to eliminate pulses that are consistent with more than one energy deposit in the pulse time window, pulses with a non-physical shape, and excessively noisy pulses that survived the previous selections (Extended Data Fig. 3). The selection efficiencies are summarized in Table 1, with details provided in Methods.

Table 1 Main parameters for the 0νββ analysis

The detector response to a monoenergetic energy deposition is an important input to the 0νββ decay search. We empirically model the response function of each calorimeter as a sum of three equal-width Gaussians and determine the function parameters from a fit to the 2,615-keV calibration line3. As a characteristic indicator of the overall energy resolution of the calorimeters we quote the exposure-weighted harmonic mean FWHM of the detectors at this calibration line, 7.78 ± 0.03 keV. All values are reported as mean ± s.d.

We quantify the scaling of energy resolution with energy and investigate energy reconstruction bias—that is, the deviation of reconstructed γ-line positions from the literature values—by fitting the calorimeter response functions to prominent γ lines in the physics data, allowing the peak means and widths to vary in the fit. The bias is allowed to scale as a quadratic function of energy, as done in our previous result32, whereas the resolution scaling has been changed to a linear function of energy, following studies showing that it was overparameterized by a quadratic scaling. The results, extrapolated to Qββ, are an exposure-weighted harmonic mean FWHM energy resolution of 7.8 ± 0.5 keV and an energy bias of less than 0.7 keV. We summarize all the relevant analysis quantities in Table 1.

Figure 4 shows the full energy spectrum along with the [2,490, 2,575] keV fit region, which contains only one background peak at 2,505.7 keV from the simultaneous absorption of two coincident γ rays from 60Co in the same crystal. We estimate that around 90% of the continuum background consists of degraded α particles from radioactive contaminants of the support structure surface, and the other approximately 10% are multi-Compton scattered 2,615-keV γ events31,33.

Fig. 4: Physics spectrum for 1,038.4 kg yr of TeO2 exposure.
figure 4

We separately show the effects of the base cuts, the anti-coincidence (AC) cut, and the pulse shape discrimination (PSD). The most prominent background peaks in the spectrum are highlighted. Inset, the region of interest after all selection cuts, with the best-fit curve (solid red), the best-fit curve with the 0νββ rate fixed to the 90% CI limit (blue), and background-only fit (black) superimposed.

We run an unbinned Bayesian fit with uniform non-negative priors on the background and 0νββ decay rates. The fit with a background-only model—that is, excluding the 0νββ component—yields a mean background rate of (1.49 ± 0.04) × 10−2 counts keV−1 kg−1 yr−1 at Qββ for a corresponding median exclusion sensitivity of \({T}_{1/2}^{0\nu } > 2.8\times {10}^{25}\,{\rm{y}}{\rm{r}}\) (90% credibility interval (CI)). The fit with the signal-plus-background model shows no evidence for 0νββ decay. We find the best fit at Γ0ν = (0.9 ± 1.4) × 10−26 yr−1 and set a limit on the process half-life of \({T}_{1/2}^{0\nu } > 2.2\times {10}^{25}\,{\rm{y}}{\rm{r}}\) (90% CI). Systematic uncertainties are included as nuisance parameters and affect both the best fit and the limit by 0.8% (Extended Data Table 1). Compared to the sensitivity, the probability of getting a stronger limit is 72%. This represents, to our knowledge, the current world-leading 0νββ sensitivity for 130Te, having improved in accordance with our increased exposure since our previous result32, and although the actual limit is weaker, it is well within the range of possible outcomes due to statistical fluctuations.

Under the common assumption of a light neutrino exchange mechanism, the 130Te half-life limit converts to a limit on the effective Majorana mass of mββ < 90–305 meV, with the spread induced by different nuclear matrix element calculations34,35,36,37,38,39,40. This limit on mββ is among the strongest in the field, proving the competitiveness of the cryogenic calorimeter technique used in CUORE. CUORE will continue to take data until it reaches its designed 130Te exposure of 1,000 kg yr.

Impact

We have demonstrated that the cryogenic calorimeter technique is scalable to tonne-scale detector masses and multi-year measurement campaigns, while maintaining low radioactive backgrounds. Next-generation calorimetric 0νββ decay searches exploiting these developments are planned. Among these, CUPID (CUORE Upgrade with Particle IDentification)41 will use the same cryogenic infrastructure as CUORE, replacing the TeO2 crystals with scintillating \({{\rm{Li}}}_{2}^{100}{{\rm{MoO}}}_{4}\) crystals and exploiting the scintillation light for greater than 100-fold active suppression of the α background42,43. In parallel, the AMoRE collaboration aims to build a large-mass calorimetric 0νββ decay experiment in Korea44. In general, the possibility to cool large detector payloads paired with the low energy thresholds achievable by cryogenic calorimeters will benefit next-generation projects at the frontier of particle physics, for example dark matter searches such as SuperCDMS45 and CRESST46, and low-energy observatories exploiting CEνNS for solar and supernova neutrino studies47 and neutrino flux monitoring of nuclear reactors48.

A serendipitous effect is that the cryogenic innovations pioneered by CUORE for 0νββ decay appear to be a solution-in-waiting for the challenges faced by the relatively young, but rapidly growing, field of quantum information technology. The need to cool increasingly large arrays of qubits to less than approximately 100 mK means there is now a commercial market for large, high-cooling-power dilution refrigerators, with some featuring technological solutions derived from CUORE. Moreover, the recent realization that ionizing radiation from natural radioactivity will be a limiting factor for the coherence time of quantum processors with an increasing number of qubits49 suggests that the one-time niche, low-radioactivity ultralow-temperature cryogenics pioneered for 0νββ decay may become mainstream in large-scale quantum computing50.

Methods

Optimum trigger and analysis threshold

The continuous data stream of CUORE is triggered with the optimum trigger, an algorithm based on the optimum filter51 characterized by a lower threshold than a more standard derivative trigger32. A lower threshold enables us not only to reconstruct the low-energy part of the spectrum, but also yields a higher efficiency in reconstructing the events in coincidence between different calorimeters, and consequently a more precise understanding of the corresponding background components52,53.

The optimum trigger transfer function of every event is matched to the ideal signal shape, obtained as the average of good-quality pulses, so that frequency components with low signal-to-noise ratio are suppressed. A trigger is fired if the filtered signal amplitude exceeds a fixed multiple of the noise root mean square (RMS), defined separately for each calorimeter and dataset. We evaluate the energy threshold by injecting fake pulses of varying amplitude, calculated by inverting the calibration function, into the data stream. We reconstruct the stabilized amplitude of the fake pulses, fit the ratio of correctly triggered events to generated events with an error function, and use the 90% quantile as a figure of merit for the optimum trigger threshold. This approach enables monitoring of the threshold during data collection, and is based on the assumption that the signal shape is not energy dependent, that is, that the average pulse obtained from high-energy γ events is also a good template for events of a few keV. The distribution of energy threshold at 90% trigger efficiency is shown in Extended Data Fig. 4.

For this work we set a common analysis threshold of 40 keV, which results in >90% trigger efficiency for the majority (97%) of the calorimeters, while at the same time allowing the removal of multi-Compton events from the region of interest through the anti-coincidence cut.

Efficiencies

The total efficiency is the product of the reconstruction, anti-coincidence, pulse shape discrimination (PSD) and containment efficiencies.

The reconstruction efficiency is the probability that a signal event is triggered, has the energy properly reconstructed, and is not rejected by the basic quality cuts requiring a stable pre-trigger voltage and only a single pulse in the signal window. It is evaluated for each calorimeter using externally flagged heater events54, which are a good approximation of signal-like events.

The anti-coincidence efficiency is the probability that a true single-crystal event correctly passes our anti-coincidence cut, instead of being wrongly vetoed owing to an accidental coincidence with an unrelated event. It is extracted as the acceptance of fully absorbed γ events at 1,460 keV from the electron capture decays of 40K, which provide a reference sample of single-crystal events.

The PSD efficiency is obtained as the average acceptance of events in the 60Co, 40K and 208Tl γ peaks that already passed the base and anti-coincidence cuts. In principle, the PSD efficiency could be different for each calorimeter, but given the limited statistics in physics data we evaluate it over all channels and over the full dataset. To account for possible variation between individual calorimeters, we compare the PSD efficiency obtained by directly summing their individual spectra with that extracted from an exposure-weighted sum of the calorimeters’ spectra. We find an average ±0.3% discrepancy between the two values and include it as a global systematic uncertainty in the 0νββ fit. This takes a Gaussian prior instead of the uniform prior used in our previous result32, which had its uncertainty come from a discrepancy between two approaches that has since been resolved.

Finally, the containment efficiency is evaluated through Geant4-based Monte Carlo simulations55 and accounts for the energy loss due to geometrical effects as well as bremsstrahlung.

Principal component analysis for PSD

In this analysis we use a new algorithm based on principal component analysis (PCA) for pulse shape discrimination. The method has been developed and documented for CUPID-Mo56, and has been adapted for use in CUORE57. This technique replaces the algorithm used in previous CUORE results, which was based on six pulse shape variables30. The PCA decomposition of signal-like events pulled from γ calibration peaks yields a leading component similar to an average pulse, which on its own captures >90% of the variance between pulses. We choose to treat the average pulse of each calorimeter in a dataset as if it were the leading PCA component, normalizing it like a PCA eigenvector. We can then project any event from the same channel onto this vector and attempt to reconstruct the 10-s waveform using only this leading component. For any waveform x and leading PCA component w with length n, we define the reconstruction error as:

$${\rm{RE}}=\sqrt{\mathop{\sum }\limits_{i=1}^{n}{({{\bf{x}}}_{i}-({\bf{x}}\cdot {\bf{w}}){{\bf{w}}}_{i})}^{2}}.$$
(1)

This reconstruction error metric measures how well an event waveform can be reconstructed using only the average pulse treated as a leading PCA component. Events that deviate from the typical expected shape of a signal waveform are poorly reconstructed and have a high reconstruction error. We normalize the reconstruction errors as a second-order polynomial function of energy on a calorimeter-dataset basis (see Extended Data Fig. 5), and cut on the normalized values by optimizing a figure of merit for signal efficiency over expected background in the Qββ region of interest. Using this PCA-based method, we obtain an overall efficiency of (96.4 ± 0.2)% compared to the (94.0 ± 0.2)% from the pulse shape analysis used in our previous results, as well as a 50% reduction in the PSD systematic uncertainty from 0.6% to 0.3%.

Statistical analysis

The high-level statistical 0νββ decay analysis consists of an unbinned Bayesian fit on the combined data developed with the BAT software package58. The model parameters are the 0νββ decay rate (Γ0ν), a linearly sloped background, and the 60Co sum peak amplitude. Γ0ν and the 60Co rate are common to all datasets, with the 60Co rate scaled by a preset dataset-dependent factor to account for its expected decay over time. The base background rate, expressed in terms of counts keV−1 kg−1 yr−1, is dataset-dependent, whereas the linear slope to the background is shared among all datasets, because any structure to the shape of the background should not vary between datasets. Γ0ν, the 60Co rate, and the background rate parameters have uniform priors that are constrained to non-negative values, whereas the linear slope to the background has a uniform prior that allows both positive and negative values.

In addition to these statistical parameters, we consider the systematic effects induced by the uncertainty on the energy bias and energy resolution59,60, the value of Qββ, the natural isotopic abundance of 130Te, and the reconstruction, anti-coincidence, PSD and containment efficiencies. We evaluate their separate effects on the 0νββ rate by adding nuisance parameters to the fit one at a time and studying both the effect on the posterior global mode \({\hat{{\Gamma }}}_{0\nu }\) and the marginalized 90% CI limit on Γ0ν.

A list of the systematics and priors is reported in Extended Data Table 1. The efficiencies and the isotopic abundance are multiplicative terms on our expected signal, so the effect of each is reported as a relative effect on Γ0ν. By contrast, the uncertainties on Qββ, the energy bias, and the resolution scaling have a non-trivial relation to the signal rate; therefore, we report the absolute effect of each on Γ0ν. The dominant effect is due to the uncertainty on the energy bias and resolution scaling in physics data. We account for possible correlations between the nuisance parameters by including all of them in the fit simultaneously.

We chose a uniform prior on our physical observable of interest Γ0ν, which means we treat any number of signal events as equally likely. Other possible uninformative choices could be considered appropriate, as well. Because the result of any Bayesian analysis depends to some extent on the choice of the priors, we checked how other reasonable priors affect our result57. We considered: a uniform prior on \(\sqrt{{{\Gamma }}_{0\nu }}\), equivalent to a uniform prior on mββ and also equivalent to using the Jeffreys prior; a scale-invariant uniform prior on logΓ0ν; and a uniform prior on 1/Γ0ν, equivalent to a uniform prior on \({T}_{1/2}^{0\nu }\).

These priors are all undefined at Γ0ν = 0, so we impose a lower cut-off of Γ0ν > 10−27 yr−1, which with the given exposure corresponds to approximately one signal event. The case with a uniform prior on \(\sqrt{{{\Gamma }}_{0\nu }}\) gives the smallest effect, and strengthens the limit by 25%, whereas the flat prior on 1/Γ0ν provides the largest effect, increasing the limit on \({T}_{1/2}^{0\nu }\) by a factor of 4. In fact, all these priors weigh the small values of Γ0ν more. Therefore, our choice of a flat prior on Γ0ν provides the most conservative result.

We compute the 0νββ exclusion sensitivity by generating a set of 104 toy experiments with the background-model, that is, including only the 60Co and linear background components. The toys are split into 15 datasets with exposure and background rates obtained from the background-only fits to our actual data. We fit each toy with the signal-plus-background model, and extract the distribution of 90% CI limits, shown in Extended Data Fig. 4.

We perform the frequentist analysis using the Rolke method61, obtaining a lower limit on the process half-life of \({T}_{1/2}^{0\nu } > 2.6\times {10}^{25}\,{\rm{yr}}\) (90% CI). The profile likelihood function for Γ0ν is retrieved from the full Markov chain produced by the Bayesian analysis tool. The non-uniform priors on the systematic effects in the Bayesian fit are thus incorporated into the frequentist result as well. We extract a 90% confidence interval on Γ0ν by treating −2log as an approximate χ2 distribution with one degree of freedom. The lower limit on \({T}_{1/2}^{0\nu }\) comes from the corresponding upper edge of the confidence interval on Γ0ν. Applying the same method to the toy experiments, we find a median exclusion sensitivity of \({T}_{1/2}^{0\nu } > 2.9\times {10}^{25}\,{\rm{yr}}\).