Search for Majorana neutrinos exploiting millikelvin cryogenics with CUORE

The possibility that neutrinos may be their own antiparticles, unique among the known fundamental particles, arises from the symmetric theory of fermions proposed by Ettore Majorana in 19371. Given the profound consequences of such Majorana neutrinos, among which is a potential explanation for the matter–antimatter asymmetry of the universe via leptogenesis2, the Majorana nature of neutrinos commands intense experimental scrutiny globally; one of the primary experimental probes is neutrinoless double beta (0νββ) decay. Here we show results from the search for 0νββ decay of 130Te, using the latest advanced cryogenic calorimeters with the CUORE experiment3. CUORE, operating just 10 millikelvin above absolute zero, has pushed the state of the art on three frontiers: the sheer mass held at such ultralow temperatures, operational longevity, and the low levels of ionizing radiation emanating from the cryogenic infrastructure. We find no evidence for 0νββ decay and set a lower bound of the process half-life as 2.2 × 1025 years at a 90 per cent credibility interval. We discuss potential applications of the advances made with CUORE to other fields such as direct dark matter, neutrino and nuclear physics searches and large-scale quantum computing, which can benefit from sustained operation of large payloads in a low-radioactivity, ultralow-temperature cryogenic environment.


Introduction
The Standard Model (SM) 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 SM: they imply non-vanishing neutrino masses, requiring an extension of the SM, and violate three accidental symmetries connected to the flavor lepton numbers L e , 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 particles 4 .This decay can be mediated by a variety of non-SM mechanisms involving Majorana neutrino masses.A minimal extension of the SM Lagrangian adds heavy Majorana neutrinos that mix with the known neutrinos to produce a set of light Majorana neutrinos, explaining the observed light neutrino masses 5 and at the same time providing a mechanism to explain the baryon asymmetry in the universe 2 .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 76 Ge 6, 7 , 136 Xe 8,9 , and 130 Te.The latter, in the form of TeO 2 cryogenic calorimeters, is employed by the Cryogenic Underground Observatory for Rare Events, CUORE 10,11 .
To fully exploit the potential of TeO 2 crystals as cryogenic calorimeters, the CUORE collaboration designed and built the largest dilution refrigerator ever constructed, capable of cooling ∼1.5 tonne of material to a temperature of ∼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 4 yr measurement campaign and the results of a new high-sensitivity 0νβ β decay search with over 1 tonne•yr of TeO 2 exposure.

The CUORE experiment
CUORE is the culmination of thirty years of 0νβ β decay searches with TeO 2 cryogenic calorimeters 12 . 130Te benefits both from a high natural isotopic abundance of ∼ 34% 13 and a high Q β β of 2527.5 keV 14 , placing the 0νβ β decay region of interest (ROI) above most natural γ-emitting radioactive backgrounds.The detector is an array of 988 nat TeO 2 cubic crystals 15 (Fig. 1) of 5 × 5 × 5 cm 3 size and ∼750 g mass, for a total mass of 742 kg, which corresponds to 206 kg of 130 Te.The array is arranged as 19 towers, each comprised of 13 floors containing 4 crystals.The crystals are operated as cryogenic calorimeters 16 at a temperature of ∼10 mK.To achieve this low-temperature environment, a novel cryogenic infrastructure -the CUORE cryostat -has been realized.
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. 1) is instrumented with a neutrontransmutation-doped germanium thermistor (NTD) 18  The detector after installation.The plastic ring was used during assembly for radon protection.Bottom right: 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 polytetrafluorethylene (PTFE) supports and the gold wires instrumenting the NTD and the heater provide the thermal link between the crystal and the heat bath, i.e. the Cu frames 17 .converts thermal pulses into electric signals and a heater 19 to inject reference heat pulses for thermal gain stabilization 20 .Thermal signals can be induced by electrons emitted in 0νβ β decays but also other background radiation, e.g.γ 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 ∼3600 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 210 Pb-depleted ancient lead recovered from a Roman shipwreck 21 (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 of 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 highefficiency cleaning procedures and manipulation protocols 22 .

Cryogenic innovation and performance
Dilution refrigerator technology was originally proposed in the '50s 23  '80s driven also by the application of cryogenic calorimeters for single particle detection 24 .Gradually, experimental volumes of the order of tens of liters capable of hosting cold masses of up to 60 kg at 10 mK temperature 17 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 breakthrough in cryogenic technology, reaching an experimental volume of ∼1 m 3 and a cold mass of 1.5 tonne (detectors, holders, shields) at 10 mK, which correspond 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 mK 24 .
The CUORE detector is hosted in a multistage cryogenfree cryostat 25 (Fig. 1), equipped with 5 pulse tube (PT) 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. 2 (top and right), 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 ∼ 90%.Fig. 2 (bottom) shows that the temperature stability achieved is at the level of 0.2% (±3σ range) over a period in excess of 1 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 ∼ 0.5 keV, so any process inducing heat dissipation ≥ 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 PTs at the 1.4 Hz operational frequency and its harmonics are particularly relevant.In CUORE, we actively tune the PT relative phases for vibration cancellation 26 (Fig. 3).This solution is transferable to any cryogenic application involving signals in the same bandwidth of the PT-induced noise.
CUORE now collects sensitive exposure with 984 out of 988 calorimeters, at a rate which is, to our knowledge, unprecedented for cryogenic particle detectors.Below, we describe the data treatment and 0νβ β decay search with >1 tonne•yr of TeO 2 exposure.Figure 3. PT phase optimization.Top: frequency spectrum of the noise for a random combination of the PT phases (orange) and after the active phase tuning (teal).The frequency spectrum of physical signals is also reported for reference.Bottom: integral of the power spectrum at the PT frequency (1.4 Hz) and its harmonics before and after active noise cancellation.

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 232 Th and/or 60 Co sources.
The voltage across each NTD is amplified, passed through an anti-aliasing filter, and continuously digitized with a 1 kHz sampling frequency 27,28 .We identify thermal pulses in the data stream and compute the pulse amplitudes, applying optimum filters that maximize the frequencydependent signal-to-noise ratio 29 .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 2615 keV 208 Tl calibration line for the remainder.Drift-stabilized pulse amplitudes are converted to energy using the regularly acquired source calibration data 30 .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 impacted by high noise or failed data processing, which amounts to 5% of the exposure.Furthermore calorimeters with > 19 keV full width at half maximum (FWHM) energy resolution at the 2615 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 (MC) simulations show that ∼ 88% of 0νβ β decay events release their full energy in a single crystal 31 .Hence, we apply an anti-coincidence (AC) cut that excludes events depositing energy in more than one crystal.Finally, we use pulse shape discrimination (PSD) to eliminate pulses 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 Tab. 1, with details provided in Methods.
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 2615 keV calibration line 3 .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, i.e. 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 result 32 , while 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 < 0.7 keV.We summarise all the relevant analysis quantities in Tab. 1. Fig. 4 shows the full energy spectrum along with the [2490, 2575] keV fit region, which contains only one background peak at 2505.7 keV from the simultaneous absorption of two coincident γ rays from 60 Co in the same crystal.We estimate ∼90% of the continuum background consists of degraded α particles from radioactive contaminants of the support structure surface, while the other ∼10% are multi-Compton scattered 2615 keV γ events 31,33 .
We run an unbinned Bayesian fit with uniform nonnegative priors on the background and 0νβ β decay rates.The fit with a background-only model, i.e. excluding the 0νβ β component, yields a mean background rate of (1.49 ± 0.04) • 10 −2 counts/(keV•kg•yr) at Q β β for a corre- 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.Top right inset: the ROI 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.
sponding median exclusion sensitivity of T 0ν 1/2 > 2.8 • 10 25 yr (90% credibility interval (CI)).The fit with the signal-plusbackground 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 0ν 1/2 > 2.2 • 10 25 yr (90% CI).Systematic uncertainties are included as nuisance parameters and affect both the best fit and the limit by 0.8% (Extended Data Tab.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 130 Te, having improved in accordance with our increased exposure since our previous result 32 , 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 130 Te 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 calculations [34][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 130 Te exposure of 1000 kg•yr.

Impact
We have demonstrated for the first time 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 utilize the same cryogenic in-frastructure as CUORE, replacing the TeO 2 crystals with scintillating Li 100 2 MoO 4 crystals and exploiting the scintillation light for > 100-fold active suppression of the α background 42,43 .In parallel, the AMoRE collaboration aims to build a large mass calorimetric 0νβ β decay experiment in Korea 44 .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, such as dark matter searches like SuperCDMS 45 and CRESST 46 , and low-energy observatories exploiting CEνNS for solar and supernova neutrino studies 47 and neutrino flux monitoring of nuclear reactors 48 .
A quite serendipitous impact 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 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 qubits 49 suggests the one-time niche, low-radioactivity ultra-low temperature cryogenics pioneered for 0νβ β decay may become mainstream in large-scale quantum computing 50 .

Online Content
Methods, additional Extended Data and Source Data are available in the online version of the paper.hai Jiao Tong University; Shanghai Laboratory for Particle Physics and Cosmology, Shanghai 200240, China 22 Wright Laboratory, Department of Physics, Yale University, New Haven, CT 06520, USA 23 Dipartimento di Fisica e Astronomia, Alma Mater Studiorum -Università di Bologna, Bologna I-40127, Italy 24 IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France 25 Lawrence Livermore National Laboratory, Livermore, CA 94550, USA 26 Department of Nuclear Engineering, University of California, Berkeley, CA 94720, USA 27

Optimum trigger and analysis threshold
The continuous data stream of CUORE is triggered with the optimum trigger (OT), an algorithm based on the optimum filter 51 characterized by a lower threshold than a more standard derivative trigger 32 .A lower threshold allows 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 components 52,53 .
The OT 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-tonoise 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 OT threshold.This approach allows monitoring of the threshold during data collection, and is based on the assumption that the signal shape is not energy dependent, i.e. that the average pulse obtained from high energy γ events is a good template also 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 ROI through the AC cut.

Efficiencies
The total efficiency is the product of the reconstruction, AC, 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 events 54 , which are a good approximation of signal-like events.
The AC efficiency is the probability that a true singlecrystal event correctly passes our AC cut, instead of being wrongly vetoed due to an accidental coincidence with an unrelated event.It is extracted as the acceptance of fully absorbed γ events at 1460 keV from the electron capture decays of 40 K, which provide a reference sample of single-crystal events.
The PSD efficiency is obtained as the average acceptance of events in the 60 Co, 40 K, and 208 Tl γ 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 result 32 , 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 MC simulations 55 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-Mo 56 , and has been adapted for use in CUORE 57 .This technique replaces the algorithm employed in previous CUORE results, which was based on 6 pulse shape variables 30 .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 9/16 vector and attempt to reconstruct the 10-second waveform using only this leading component.For any waveform x and leading PCA component w with length n, we define the reconstruction error as: 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 β β ROI.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 package 58 .The model parameters are the 0νβ β decay rate (Γ 0ν ), a linearly sloped background, and the 60 Co sum peak amplitude.Γ 0ν and the 60 Co rate are common to all datasets, with the 60 Co 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•kg•yr), is dataset-dependent, while the linear slope to the background is shared among all datasets since any structure to the shape of the background should not vary between datasets.Γ 0ν , the 60 Co rate, and the background rate parameters have uniform priors constrained to non-negative values, while the linear slope to the background has a uniform prior allowing 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 resolution 59,60 , the value of Q β β , the natural isotopic abundance of 130 Te, and the reconstruction, AC, 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 Γ0ν and the marginalized 90% CI limit on Γ 0ν .
A list of the systematics and priors is reported in Extended Data Tab. 1.The efficiencies and the isotopic abundance are multiplicative terms on our expected signal, so the impact of each is reported as a relative effect on Γ 0ν .In 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.Since the result of any Bayesian analysis depends to some extent on the choice of the priors, we checked how other reasonable priors affect our result 57 .We considered: • a uniform prior on √ Γ 0ν , equivalent to a uniform prior on m β β and also equivalent to using the Jeffreys prior; • a scale-invariant uniform prior on log Γ 0ν ; • a uniform prior on 1/Γ 0ν , equivalent to a uniform prior on T 0ν 1/2 .
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 √ Γ 0ν gives the smallest effect, and strengthens the limit by 25%, while the flat prior on 1/Γ 0ν provides the largest effect, increasing the limit on T 0ν 1/2 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 10 4 toy experiments with the background-model, i.e. including only the 60 Co 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-plusbackground model, and extract the distribution of 90% CI limits, shown in Extended Data Fig. 4.
We perform the frequentist analysis using the Rolke method 61 , obtaining a lower limit on the process half-life of T 0ν 1/2 > 2.6 • 10 25 yr (90% CI).The profile likelihood function L 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 −2 log L as an approximate χ 2 distribution with one degree of freedom.The lower limit on T 0ν 1/2 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 0ν 1/2 > 2.9 • 10 25 yr.

Data Availability
The data generated during this analysis and shown in paper figures are available in ASCII format (CSV) as Source Data in the repository Extended Data Table 1.Systematics affecting the 0νβ β decay analysis.The total analysis efficiency is the product of all the efficiencies listed in Tab. 1 except containment.The PSD efficiency refers to its additional systematic uncertainty described in the text.The first four systematics are multiplicative effects and the impact of each is presented as a percentage.The last two systematics have a non-trivial effect on Γ 0ν , hence we report the absolute effect.We report the variation induced on the marginalized 90% CI limit (third column) and the posterior global mode Γ0ν (last column).In calibration data, the AC is not applied in order to maximize the statistics on the γ peaks, and the PSD mostly removes pileup events (events with more than one energy deposit in the time window).In physics data, the PSD mostly eliminates random noise events, which can correspond to either physical events with excessive noise or to noise-induced events with non-physical pulse shapes.Such events appear randomly across the energy spectrum, so the cut mostly acts on the continuum.The distribution contains only events that passed the other base cuts.The second order polynomial fit is shown in orange.Right: two example pulses from this calorimeter.The actual pulse is drawn in teal, and the corresponding reconstruction obtained by the PCA is drawn in orange.The top pulse deviates from the expected shape of a good pulse and is rejected, while the bottom one conforms to the expected response and is accepted.

Figure 1 .
Figure 1.The CUORE detector.Left: Rendering of the 6-stage cryostat, with the pulse tubes and dilution unit, the internal low-radioactivity modern and Roman lead shields, and the array of 988 TeO 2 crystals (light blue).Top right:The detector after installation.The plastic ring was used during assembly for radon protection.Bottom right: 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 polytetrafluorethylene (PTFE) supports and the gold wires instrumenting the NTD and the heater provide the thermal link between the crystal and the heat bath, i.e. the Cu frames17 .

Figure 2 .
Figure 2. Cryogenic performance.Top: The exposure accumulated by CUORE (left, teal), along with the exposure used for this analysis (left, orange).Part of 2017 and 2018 was dedicated to maintenance and optimization of the cryogenic setup; since then, CUORE has been operating stably with a 90% duty cycle (March 2019-October 2020) (right).Middle: Examples of temperature instabilities induced by external causes, e.g.blackouts and earthquakes, or human intervention, such as regular maintenance or the insertion of calibration sources.Bottom: 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).

Figure 4 .
Figure 4. Physics spectrum for 1038.4 kg•yr of TeO 2 exposure.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.Top right inset: the ROI 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.

Extended Data Figure 1 .
Working principle of the cryogenic calorimeter.Left: simplified calorimeter thermal model.The detector is modeled as a single object with heat capacity C coupled to the heat bath (with constant temperature T 0 ) through the thermal conductance G.The NTD thermistor for signal readout is glued to the absorber.Right: Example of a CUORE pulse from the 2615 keV calibration line: T 0 corresponds to the baseline height, the pulse amplitude is proportional to the deposited energy, and the decay time depends on the C/G ratio.

Extended Data Figure 3 .
©Aquarium -Archivio Soprintendenza Archeologica della Sardegna Extended Data Figure 2. Roman lead.Top left: the lead bricks recovery from the Sardinian sea.Bottom left: the ingot inscriptions were cut and preserved, while the ingot bodies were used for the CUORE internal lead shield.Right: Lateral view of the internal lead shield 21 .Pulse Shape Discrimination.Effect of the PSD cut on calibration data around the 2615 keV line (left) and on physics data near Q β β (right).

Extended Data Figure 4 . 5 .
Optimum trigger and statistical analysis.Top left: Distribution of energy thresholds at 90% trigger efficiency for the OT algorithm in a single dataset.The 40 keV analysis threshold is indicated by the vertical line.Top right: 90% CI exclusion limits on T 0ν 1/2 from an ensemble of 10 4 toy experiments generated with the background-only model, with background rates obtained from the background-only fit to the data.The median exclusion sensitivity is indicated by the orange line.Bottom left: Posterior probability distribution for Γ 0ν obtained from the Bayesian fit, with the 90% CI highlighted.Bottom right: ∆χ 2 values obtained from the profile likelihood of Γ 0ν , with ∆χ 2 = 0 being the most-favored value.The frequentist limit at 90% confidence level (CL) is indicated.PCA performance.Left: example of a normalization fit of the PCA reconstruction error vs energy for a single calorimeter and dataset.

Table 1 .
Main parameters for the 0νβ β analysis.The resolution and efficiencies are exposure-weighted average values.