of the 2νββ Decay Half-Life of ^{130}Te with CUORE.

We measured two-neutrino double beta decay of 130 Te using an exposure of 300.7 kg · yr accumu-lated with the CUORE detector. Using a Bayesian analysis to ﬁt simulated spectra to experimental data, it was possible to disentangle all the major background sources and precisely measure the two-neutrino contribution. The half-life is in agreement with past measurements with a strongly reduced uncertainty: T 2 ν 1 / 2 = 7 . 71 +0 . 08 − 0 . 06 (stat . ) +0 . 12 − 0 . 15 (syst . ) × 10 20 yr. This measurement is the most precise determination of the 130 Te 2 νββ decay half-life to date.


INTRODUCTION
Two-neutrino double beta (2νββ) decay is a nuclear transition with the longest lifetime experimentally measured. This process occurs when two neutrons in a nucleus simultaneously decay emitting two anti-neutrinos and two electrons. This decay can act as background for a hypothetical process called neutrinoless double beta (0νββ) decay [1], which may occur if the neutrino were a Majorana fermion [2], and which violates lepton number conservation [3]. 0νββ decay would be new physics and could explain the origin and nature of the neutrino mass states [4][5][6][7].
The CUORE experiment primarily searches for neutrinoless double beta decay (0νββ) of 130 Te → 130 Xe + 2e − [22,23], however other searches are possible [24][25][26][27]. This letter will outline first the CUORE detector, the data collection, and the analysis of the 2νββ decay of 130 Te. Thereafter, we provide a description of the technique used to fit the experimental data (comprised of events from both the 2νββ decay and the background sources), and finally we present a discussion of the fit results.

CUORE DETECTOR
CUORE is located underground at the Gran Sasso National Laboratory of INFN, Italy, with ∼3600 meter water equivalent overburden to shield from cosmogenic backgrounds [28,29]. The CUORE detector consists of 988 TeO 2 crystals [30] arranged in 19 towers made of 13 floors of 4 crystals, read out as individual channels. Each crystal is 5×5×5 cm 3 in size and ∼750 g in mass, with a natural abundance of ∼34% for 130 Te. The crystals are cooled to ∼10 mK at which point they have a low heat capacity and can be operated as cryogenic bolometers. The energy deposited in the crystal by particle interaction causes a temperature increase which is measured via a neutron transmutation doped (NTD) Ge thermistor [31]. The signal rise time is ∼100 ms, and the high energy resolution is second only to that of Ge detectors. Si heaters are used to inject regular reference pulses for a detector-based correction of thermal gain drifts from long term temperature variations.
The detector is housed in a large cryogen-free cryostat, cooled to ∼10 mK by a dilution refrigerator [32]. The cryostat and the detector were constructed with strict radiopurity controls [30] in order to reduce the α and γ backgrounds seen in Cuoricino [33], and further reduce the γ background seen in CUORE-0 [34]. The cryostat is equipped with two lead shields: a 6-cm thick shield of ancient Roman lead [35] at ∼4 K around and below the detector, and a 30-cm thick shield at ∼50 mK located above the detector. An additional external shield is comprised of a 25-cm thick layer of lead surrounded by a 20-cm thick layer of polyethylene. A layer (2 cm thick) of boric acid that absorbs thermalized neutrons is located between the two shields.
When 2νββ decay occurs inside a bolometer the neutrinos escape without interacting, thus we detect the two electrons sum kinetic energy forming a continuous, βlike spectrum from 0 keV up to the Q-value of the decay (Q ββ = 2527 keV [36][37][38]). Background contributions to this spectrum originate from radioactivity in the detector and cryostat components. These backgrounds can be disentangled and quantified via careful analysis of the observed spectral shape and topological information in the segmented CUORE detector in comparison to a detailed background model [39,40].

DATA COLLECTION
CUORE began taking data in early 2017. The data collected through mid 2019 are analyzed in this work and are grouped into 7 datasets (physics data bounded by calibration data). The 2νββ decay analysis requires high quality data over the whole energy range, specifically channels need to be both well-performing and wellcalibrated. This led us to exclude 2 datasets from the analysis given that the large majority of channels did not satisfy these criteria. With this choice of dataset-channel we have 300.7 kg·yr of TeO 2 exposure (102.7 kg·yr of 130 Te exposure).
The data itself are a collection of events, corresponding to a triggered waveform on a single bolometer. The modularity of the detector allows us to reconstruct the event topology via a time based coincidence analysis. Events are grouped into multiplets, M i (i = number of triggered bolometers) if they occur within a ±30 ms window on bolometers that are ≤ 15 cm apart from each other, with a minimum energy of 70 keV. Given the extremely The observed M1 spectrum (black ) compared with its reconstruction as obtained by the background model (blue). The reconstructed 2νββ decay component is shown in yellow for comparison. We observe that from 900 to 2000 keV more than 50% of the M1 spectrum counts originate from the 2νββ decay process. The experimental data and the spectra reconstructed by the fit have been converted back to 1 keV binning for illustrative purposes. Selected gamma lines from background contaminants are labeled.
low trigger rate of CUORE bolometers (∼1 mHz) and the distance requirement, multiplets with i > 1 contain practically no accidental coincidence events and are mainly induced by particles depositing energy in multiple crystals.
We split the data into three types of spectra: a multiplicity 1 (M 1 ) spectrum comprised of events where energy was deposited into a single bolometer, a multiplicity 2 (M 2 ) spectrum comprised of the single energies detected by each of the two bolometers simultaneously triggered, and a Σ 2 spectrum comprised of the sum energy of the M 2 events. The energy of a 2νββ decay event is deposited into a single bolometer with a probability obtained from Monte Carlo simulations of ∼ 90%. The majority of backgrounds deposit energy across two or more bolometers (such as γ's that scatter from one crystal into another or α decays that occur on a surface between two crystals), making the M 2 and Σ 2 spectra useful for understanding backgrounds. Events with multiplicity higher than 2 are not considered in this analysis since they do not add new information.

SPECTRAL FIT
We analyze the events with energies from a threshold of 350 keV to 2.8 MeV, where the M 1 spectrum is dominated by 2νββ decay (between 900 to 2000 keV the contribution exceeds 50% of total M 1 events) along with γ/β emissions from radioactive contaminants. To disentangle the 2νββ decay signal we construct a background model (BM) that describes the data via a comprehensive list of possible sources. Guidelines for this work are taken from the CUORE-0 BM [39] and the CUORE background budget [41]. The background sources are radioactive contaminations located both in the bulk of the detector and cryostat components, on the surfaces of crystals, and materials with a line of sight to them. We also include cosmogenic muons.
We developed a Geant4 [42] Monte Carlo (MC) simulation [39,41] which outputs the spectra produced by each source in the detector, reproducing all relevant features of the experimental data (e.g., multiplicity, time resolution, energy dependent trigger efficiencies, etc.). The elements of the CUORE experiment, including the cryostat, are grouped into 9 geometric entities used in the model as background source positions: the crystals, the copper structure holding the towers, the copper ves-sel enclosing the detector, the Roman lead, the internal and external shields made of modern lead, the cryostat thermal shields (grouped into two elements: the copper shields inside the Roman lead and those outside), and finally the internal lead suspension system. The BM uses 62 simulated sources to fit the experimental data. One is the 2νββ decay in the crystals, and 60 others refer to different contaminants in the 9 elements listed above. These include bulk and surface 238 U and 232 Th contaminations (allowing for secular equilibrium breaks), bulk 60 Co, 40 K, and a few other long lived isotopes, as indicated in Fig. 1. All these isotopes are identified from the presence of one or more characteristic γ lines in the observed spectra. The only exception is 90 Sr, a long-lived pure β emitter, that could be present due to a hypothesized contamination by radioactive fallout. The remaining simulation (number 62) is the cosmogenic muon flux. As described in [39] a variable binning is applied to all the spectra: the minimum bin size is 15 keV, and bins with less than 30 counts are merged. All counts belonging to a single γ line are combined into a single bin to avoid systematics from the modeling of the γ peak shapes in MC simulations. Finally, the trigger efficiency vs. energy and the efficiency of quality cuts are included in the analysis as global parameters.
The observed spectrum is reconstructed by simultaneously fitting a linear combination of the 62 MC simulated spectra to the M 1 , M 2 , and Σ 2 data. The fit is done with a Bayesian approach using a Markov-Chain Monte Carlo (MCMC), implemented in the JAGS software package, to sample the joint posterior probability density function (PDF) of the fit parameters [43][44][45][46]. The likelihood is a product, over bins and spectra, of Poisson distributions that give the probability of drawing the experimental counts as a function of the MC spectra normalizations. To prevent bias while tuning data quality cuts and setting the fitting procedure, the MC normalization coefficient was blinded to keep the extracted 2νββ decay half-life in terms of a nonphysical ratio that could not be compared to previous results.
For each source, except cosmogenic muons, a uniform prior is used. For muons, additional information is gained from the high multiplicity spectra (M > 5) where muons become dominant, and is used to extract a Gaussian prior for the BM fit. The fit result is a joint posterior PDF for the 62 parameters from which we extract the marginalized posterior PDF for the 2νββ decay rate. The fitting procedure closely follows the description in [39] and a paper detailing the CUORE BM is in preparation.

MODEL SYSTEMATICS
The background model is able to reproduce the major features of the observed spectra (see Fig. 1) with a global χ 2 /d.o.f of 681/365 (χ 2 red = 1.87). The sub-optimal agreement between the data and the MC likely arises from an imperfect modelling of source position and distribution. Increased statistics from more data will allow for refinement of the background model by better identifying source locations or additional sub-dominant contaminants.
In order to check the stability of the 2νββ decay halflife result we run multiple fits over the whole dataset varying aspects of the background model. In particular we test two different models for the 2νββ decay spectral shape, we alter the list of background sources used and we remove the 90 Sr source. As additional probes of our sensitivity to various aspects of the BM sources we fit subsets of data in which we split the detector in half in different ways (see Fig. 2) and perform the fit on single datasets. -2νββ Decay Model There are two competing models for 2νββ decay, yielding slightly different spectral shapes [10,11,47]. The Single State Dominance (SSD) is the default used in this analysis. It results in a better fit quality, and might be the first experimental hint for SSD dominance in 130 Te 2νββ decay. The alternate mechanism, Higher State Dominance (HSD), yields a slightly worse fit and a few percent increase of the 2νββ half-life. As this is the only shape test we perform on the 2νββ decay spectrum, we conservatively assume this systematic to be double-sided and estimate the uncertainty as 68% of the difference between the fits with HSD and SSD: ±1.3%.
-Energy Threshold The energy threshold used in this analysis is 350 keV. If we vary this threshold in the range of 300-800 keV we observe an increase of the 2νββ decay rate. We assume this systematic to be uniform between the best fit and the value that deviates the most. We symmetrize this uncertainty around the best fit to account for possible deviations given by untested threshold values, with a result of ±0.4%.
-Geometrical effect All contaminants in the model are uniformly distributed in the 9 simulated elements. To investigate possible biases we compare fits done by splitting the detector according to crystal, floor, or tower number (see Fig. 2). Each pair of results are statistically compatible with each other to within two sigma. We take the pair with the largest splitting, subtract the statistical error and interpret the result as the 1σ uncertainty of a flat distribution: ±0.8% -90 Sr As mentioned, a background source due to possible crystal contamination with 90 Sr was introduced. This is the only long-lived pure β emitter produced by fission that produces a background (via its daughter 90 Y) extending up to 2.2 MeV, without any associated gamma emission that would allow us to constrain its activity [17]. The 2νββ decay result is weakly sensitive to this contaminant, as upon its removal the counts ascribed to 2νββ decay increase resulting in a slightly shorter halflife. Since 90 Sr has no clear signature we use it as a proxy for the removal or addition of components to the BM. We take the systematic to be symmetric and 68% of the difference between the best fit and the fit without, giving ±0.3%.
-Datasets We investigated the 2νββ decay result stability in time by fitting separately each of the 5 datasets used in this analysis and observed only statistical variations in the fit result. We also fit the two excluded datasets and use the result to quantify the bias introduced by their removal. This yields an asymmetric uncertainty of +0.3% and -1.1%.
-List of background sources In the reference fit we use 62 sources, selected to be comprehensive, as extracted from the CUORE-0 background model [39]. Given the limited statistics used in this work and the choice of fitting only the γ region some of the sources could be degenerate with each other while others cannot be identified easily. To check how the fit performs with a different background source list, all components with contributions compatible with 0, including 90 Sr, are removed. This reduces the number of distinct background source components down to 25. This has an impact on the 2νββ decay fit result compatible to that resulting from the removal of the 90 Sr alone.
As a result of this study we can conclude that all the systematics we explored are at most in the range of 1%. The dominant contribution comes from the uncertainty in the decay model (SSD vs HSD) which may be improved with increased statistics or theoretical input. Finally, other sources of uncertainties such as the efficiency of our coincidence selection, the chance of mixing up M 1 with M 2 events, or the efficiency of pulse shape cuts, have an overall impact on the final error that is lower than 0.1%. The Monte Carlo statistics, though optimized to yield negligible error, is properly accounted for in the fitting procedure.

2νββ DECAY RESULTS AND DISCUSSION
To extract a robust estimate of the 2νββ decay halflife, we combine our systematics in quadrature. Through the unblinding of the correct normalization coefficient for the MC spectrum, we obtain the measurement of the 2νββ decay half-life of 130 Te. Though the posterior for 90 Sr is compatible with null activity, the insertion of 90 Sr in the BM does weakly distort the 2νββ decay posterior. Removal of the 90 Sr source results in a symmetric posterior, however we choose to include this source due to the high anti-correlation with 2νββ decay.

+0.12
−0.15 (syst.) × 10 20 yr, a value consistent with previous measurements (see Table I). This result is the most precise measurement of the 2νββ decay half-life of 130 Te to date and one of the most precise measurements of a 2νββ decay half-life. It represents a substantial improvement over previous measurements from NEMO-3 [49] and CUORE-0 [39] owing to the CUORE strict radiopurity controls, the improved signal-to-noise ratio, the increased statistics, and the robust background model.

CONCLUSION
In this paper we described the analysis of the 130 Te 2νββ decay measured with CUORE. We exploit the ge-ometry of the CUORE detector to tag single scatter and multiple scatter events to obtain separate spectra dominated by 2νββ decay and background events, respectively. The 130 Te 2νββ decay half-life is measured to be T 2ν 1/2 = 7.71 +0.08 −0.06 (stat.)

+0.12
−0.15 (syst.) × 10 20 yr. Compared to previous results (Table I) this is the most precise determination of the 2νββ decay half-life in 130 Te. The present result is dominated by a ∼2% systematic uncertainty. Further improvement will require a better understanding of the background sources localization, as indicated by the observed BM variations with different geometrical detector splittings. This refinement, as well as an improved study of the SSD vs HSD models, is feasible in the near future given the increased statistics being collected by CUORE.

ACKNOWLEDGMENTS
The CUORE Collaboration thanks the directors and staff of the Laboratori Nazionali del Gran Sasso and the technical staff of our laboratories. This work makes use of both the DIANA data analysis and APOLLO data acquisition software packages, which were developed by the CUORICINO, CUORE, LUCIFER and CUPID-0 Collaborations.

Supplemental Materials
Here we present additional plots to illustrate key results from the CUORE background model fits with respect to the 2νββ decay half-life result. In particular we show the fit results for the 3 spectra used in the fits (M 1 , M 2 , Σ 2 ) which show good agreement between data and model. We also include a comparison of the 2νββ decay posterior with and without the 90 Sr source to illustrate that, while there is some slight distortion, the overall impact is quite negligible. This is further supported by the posterior of the 90 Sr contribution.
The M 1 spectrum (Fig. 3) is comprised of singlecrystal events which contain a significant contribution from 2νββ decay. The fit residuals show that in the region of 1-2 MeV the reconstructed spectrum matches the observed data quite well.
In Figs. 4 and 5 we show two views of the M 2 data: a spectrum from the individual components of the M 2 multiplets (i.e. the M 2 spectrum), and a spectrum from the sum of the two components (i.e., the Σ 2 spectrum). The M 2 spectrum shown in Fig. 4 displays the energy spectra from individual components of M 2 multiplets. Events in the M 2 spectrum provide useful information on the localization of contaminations, as the strict coincidence criterion described earlier limits these to the nearest neighbor crystals. The Σ 2 spectrum in Fig. 5 shows the summed energy of each member of an M 2 multiplet.
The Σ 2 events show more clearly the γ lines that produce the physical interactions in two crystals. We see here, that in the γ region the background model accurately reconstructs the observed spectrum. Minor disagreements in some of the peaks is attributed to a potential for further improvement in source localization throughout the detector.
There is a possible contribution to the background from a fission product, 90 Sr. This is a pure β emitter with a decay energy of 0.564 MeV to 90 Y , which is another nearly pure β emitter with a 2.28 MeV endpoint allowing the decay chain to contribute to the background up through this energy. The net effect is an anti-correlation in the fit between the 90 Sr rate and the 2νββ decay rate. Systematics checks (described earlier) show that toggling this potential source off alters the 2νββ decay rate by ∼0.6%. The inclusion of the 90 Sr source causes a slight asymmetry in the 2νββ decay posterior (Fig. 6(a)). By fitting this distribution with a 2-sided Gaussian we get the left and right uncertainty range. Without 90 Sr the posterior is far more symmetric requiring only a single Gaussian to fit. The half-life from this analysis with the 90 Sr source is 7.71 +0.08 −0.06 ×10 20 yr, compared to the result without the source: 7.67±0.05×10 20 yr. As there is an anti-correlation between the 2νββ and 90 Sr decay rates we keep this source in the final fit to get a conservative result. An examination of the 90 Sr posterior itself ( Fig. 6(b)) shows that this is a conservative approach as the posterior indicates a most probable value of zero.

FIG. 3. Top:
The measured M1 spectrum (blue) and its reconstruction (red ). The spectra are binned with an adaptive binning to contain peaks into a single bin (to avoid dependence on the peak shape), while also achieving good resolution of the continuum shape. Bottom: The ratio of the data to the reconstructed model with 1σ, 2σ and 3σ error bars. It is clear from the data that we are able to faithfully reconstruct the continuum and peaks from sources. FIG. 4. The M2 spectrum provides an indication of the localization of sources as they are populated by coincident events from neighboring crystals. Top: The measured M2 spectrum (blue) and its reconstruction (red ). The spectra are binned with an adaptive binning to contain peaks into a single bin (to avoid dependence on the peak shape), while also achieving good resolution of the continuum shape. Bottom: The ratio of the data to the reconstructed model with 1σ, 2σ and 3σ error bars. FIG. 5. The Σ2 spectrum better indicates the contributions to the background that originate from various γ lines as a result of interactions with two crystals. Top: The measured Σ2 spectrum (blue) and its reconstruction (red ). The spectra are binned with an adaptive binning to contain peaks into a single bin (to avoid dependence on the peak shape), while also achieving good resolution of the continuum shape. Bottom: The ratio of the data to the reconstructed model with 1σ, 2σ and 3σ error bars.   6. (a) Posterior of the 2νββ decay reference fit (red line) compared with the posterior of the 2νββ decay fit with the 90 Sr source turned off (blue line). The latter posterior is far more symmetric due to the anti-correlation between the two sources causing a distortion in the reference fit. (b) Posterior for the contribution from 90 Sr in the reference fit, with the 90% C.I. in yellow. The posterior peaks at a value consistent with 0 activity, indicating that the contribution from this source on the 2νββ decay half-life measurement is negligible. Since there is a slight anti-correlation and distortion of the 2νββ posterior we make a conservative choice to include it in the model for the 2νββ decay fit result.