Measurement of the Two-Neutrino Double Beta Decay Half-life of $^{130}$Te with the CUORE-0 Experiment

We report on the measurement of the two-neutrino double beta decay half-life of $^{130}$Te with the CUORE-0 detector. From an exposure of 33.4 kg$\cdot$y of TeO$_2$, the half-life is determined to be $T_{1/2}^{2\nu}$ = [8.2 $\pm$ 0.2 (stat.) $\pm$ 0.6 (syst.)] $\times$ 10$^{20}$y. This result is obtained after a detailed reconstruction of the sources responsible for the CUORE-0 counting rate, with a specific study of those contributing to the $^{130}$Te neutrinoless double beta decay region of interest.


Introduction
Double-beta decay is a rare nuclear process in which two nucleons simultaneously decay and emit two electrons. The allowed Standard Model version of this process emits two (anti-)neutrinos and is called two-neutrino double-beta decay (2νββ). This decay is interesting in a E-mail: cuore-spokesperson@lngs.infn.it b Deceased its own right as the slowest process ever directly observed [1,2]. Moreover it may represent an important source of background for the neutrinoless double-beta decay (0νββ), i.e. a related process with no neutrino emission [3]. 0νββ manifestly violates lepton number and therefore its discovery would point to new physics beyond the Standard Model. Experiments searching for 0νββ have made great leaps forward in sensitivity by using a variety of techniques and isotopes [1,2,4,5].
The Cryogenic Underground Observatory for Rare Events (CUORE) [6,7] is the latest and most massive in a family of bolometric detectors designed to search for the 0νββ decay of 130 Te (Q-value= 2528 keV [8,9,10]). The detector combines the excellent energy resolution achievable with the bolometric technique (∼ 5 keV at 2615 keV) with the exceptionally high natural abundance of 130 Te (34%). The first phase of CUORE, named CUORE-0, was 1/19 th the size of CUORE and operated at the Laboratori Nazionali del Gran Sasso (LNGS), in Italy, between 2013 and 2015. In addition to being a competitive 0νββ experiment [11,5,12], CUORE-0 is a test of the assembly protocols for CUORE: the reconstruction of background sources responsible for the CUORE-0 counting rate enables us to verify that the necessary background requirements for CUORE are fulfilled.
In bolometers, the neutrinos emitted in 2νββ are not detected, and the summed kinetic energy of the two electrons forms a continuous spectrum from 0 keV up to the Q-value of the decay. Conversely, 0νββ produces no neutrinos and the experimental signature is a sharp peak at the Q-value of the decay, broadened by the energy resolution of the detector. This broadening smears 2νββ events into the 0νββ region of interest (ROI) around the Q-value and forms an irreducible background to the 0νββ signal; thus a good energy resolution is key to mitigating this background. The other background contributions in the ROI come from naturally occurring radioactivity in the detector components. These background sources can be disentangled and described quantitatively by carefully analyzing the shape of the measured spectrum and constructing a detailed background model, including both physics processes and instrumental effects. This paper reports the CUORE-0 background model as well as a new precision measurement of the 130 Te 2νββ half-life. We use a detailed Geant4-based simulation and a Bayesian fitting algorithm with a priori constraints from materials assay to reconstruct the experimental data and form a posteriori estimates of the background source activities. A frequentist analysis, using the same model, is presented in [13]. We present the experimental details, including the data acquisition and analysis chain, in Sec. 2. Background constraints from external measurements are summarized in Sec. 3. The construction of the Monte Carlo simulation code is presented in Sec. 4. The validation of Monte Carlo simulations, by comparing external radioactive source calibration spectra with the data, is presented in Sec. 5. The set of identified background sources and their effects on the experimental data is discussed in Sec. 6. The Bayesian fitting tool is introduced in Sec. 7. We present the fit results and discuss their uncertainties in Sec. 8. Finally, the 2νββ decay half-life evaluation is presented in Sec. 9 and the 0νββ ROI reconstruction is discussed in Sec. 10.

Experiment
The CUORE-0 detector is one tower of 52 nat TeO 2 crystals arranged in 13 floors, each with 4 crystals [14]. Each crystal is 750 g and is operated as an independent bolometer at ∼10 mK. At these temperatures, the interaction of a particle with the crystal generates a measurable temperature rise proportional to the energy deposited. The total detector mass is 39 kg of TeO 2 , or 10.8 kg of 130 Te. The tower is situated in a dilution refrigerator that provides the cooling power necessary to keep the TeO 2 crystals at their working temperature. The copper and PTFE support structure holds the crystals and provides the thermal link to the refrigerator. The tower is surrounded by several layers of shielding, including low-background Roman lead and an anti-radon box. A schematic of the experiment is shown in Fig. 1 (Left). The crystals temperatures are measured continuously using neutron transmutation doped (NTD) germanium thermistors, coupled to the crystals using epoxy. The thermistors convert temperature variations to a voltage output, which is digitized at a rate of 125 Samples/s. A software trigger is used to identify events and collect them in 5 s windows. Each window is divided into two periods: 1 s before the trigger and 4 s after it. The period before the trigger is used to establish the baseline temperature of the crystal and the remaining 4 s is used to determine the pulse amplitude. Together these are used to extract the energy deposited. A silicon resistor is also coupled to each crystal with epoxy and is used to generate reference thermal pulses every 300 s. These are used to stabilize the gain of the bolometer against temperature fluctuations. Forced triggers are used for noise and threshold studies. The threshold (in the default high-energy triggered data used for this analysis) is bolometer dependent and is approximately between 30 and 120 keV. The energy response is calibrated by inserting thoriated tungsten wires inside the outer vessel of the cryostat, and using the γ lines from the 232 Th decay chain to calibrate each bolometer independently. The details of the CUORE-0 detector design, operation and performance are described in [14]. New protocols for material selection, cleaning, and handling were developed for the detector crystals and tower support structure. The dilution refrigerator, shielding, and other cryostat components are those from the Cuoricino experiment [15,16].

Data Production
This analysis uses data collected with CUORE-0 from March 2013 to March 2015. The data are grouped in datasets, which last approximately one month. Each dataset has approximately three days of calibration data at both the beginning and end, collected while the detector was exposed to a radioactive calibration source. In total, 20 datasets were used in the 0νββ analysis: however, in the present analysis we exclude the first dataset because of a substantial (∼15 keV) miscalibration in the highest energy region (above 3 MeV). The physics data (i.e. excluding the calibration data of each dataset) use a total of 33.4 kg·y of TeO 2 exposure, or  9.3 kg·y of 130 Te exposure, after all data quality selections. The data production converts the data from a series of triggered waveforms into a calibrated energy spectrum. The details of the data production can be found in [12], but we outline the general procedure here. The entire waveform is used to extract the amplitude of the pulse. Then, each bolometer in each dataset is calibrated independently. A time-coincidence analysis is performed to search for events that deposit energy across multiple bolometers. Finally, the 0νββ ROI (2470-2580 keV) is blinded for analysis. Once the data are blinded, we implement a series of event selection cuts to maximize our sensitivity to physics events. We exclude periods of cryostat instability and malfunction. We enforce a pile-up cut of 7.1 s around each event: the 3.1 s before and 4 s after the event. We apply a series of pulse shape cuts to reject deformed or non-physical events mostly contaminating the low energy region below ∼400 keV. Double-beta decay events are usually confined within the crystal they originated from. However, many background sources deposit energy in multiple crystals within the response time of the detector. We use this information in the analysis by forming multiplets of events that occur within a coincidence window of ±5 ms in different crystals. Since the event rate is approximately 1 mHz, the probability of accidental (i.e. causally unrelated) coincidences is extremely small (∼ 10 −5 ). We then build energy spectra from these multiplets: -M spectrum is the energy spectrum of all events, with no coincidence criteria applied; -M 1 spectrum is the energy spectrum of the events with the requirement that only one bolometer triggered (multiplicity 1 or M 1 events); -M 2 spectrum is the energy spectrum of the events with the requirement that two bolometers triggered (multiplicity 2 or M 2 events); -Σ 2 spectrum is the energy spectrum associated to M 2 multiplets, each multiplet produces an entry with an energy E(Σ 2 ) that is the sum of the energies of the two events composing the multiplet.
Higher-order multiplets are used only to evaluate the contribution of muons to the background. The signal cut efficiency as a function of energy is defined as the fraction of true signal events that pass all the event cuts. In [12], we calculated the efficiency of these cuts by measuring their effect on γ-peaks in the energy spectrum. This takes advantage of the fact that the events in the γ-peaks are a nearly pure sample of true signal events. In the present analysis we use a new technique, which takes advantage of the coincidence analysis and allows a better reproduction of the energy dependence of the efficiency. Since accidental coincidences are negligible, M 2 events provide a pure sample of good events on the whole energy spectrum. The cut efficiency ε C is modeled as an exponential function of the energy, ε C (1−e −E/Ec ), which is fitted to M 2 events on a bin-by-bin basis; it rapidly reaches a stable value of ε C = 0.943 ± 0.002 at energies above E c ∼ 100 keV. This is consistent with the signal cut efficiency quoted in [12] of 0.937 ± 0.007. Above E ∼ 7 MeV the M 2 spectra have insufficient statistics to give meaningful fit results. Therefore, in this analysis we do not include events with energy above 7 MeV.

Background Sources
In rare event searches, some backgrounds are ubiquitous due to the natural decay of 40 K and the daughters of the 232 Th and 238 U decay chains, including the surface implantation of 210 Pb from environmental 222 Rn. Based on the results of previous bolometric experiments [15,16,17], we expect these radioactive contaminants to be located in the whole experimental setup, including the detector itself. The cosmogenic activation, especially of copper and tellurium, resulting in 60 Co is also of concern. A small contribution from cosmic muons [18], environmental gamma rays (γ) [19], and neutrons interacting directly in the detector is expected [19]. With this in mind, extreme care was taken into the selection of the materials used to build the CUORE-0 detector [14], and in the cleaning of all surfaces facing each TeO 2 crystal [20,21]. A lot of effort was also devoted to: the detector design, to minimize both the total mass of the inactive parts of the detector tower and the surface area facing the array; the detector assembly procedure under controlled atmosphere [14]; and the optimization of the production protocol of every detector component to limit exposure to cosmic rays. Tables 1 and 2 report the bulk and surface activities from the screening of the CUORE-0 components. A discussion about the different assay techniques that we used to derive the quoted contamination limits can be found in [17,20]. The Cuoricino results (this experiment used the same cryostat and shield as CUORE-0) and screening results guide the construction of the background model and the definition of the priors on material contaminants.

Monte Carlo
The background sources are simulated using a Geant4based Monte Carlo code called MCuoreZ . The code generates and propagates primary and any secondary particles through the CUORE-0 geometry until they are detected in the TeO 2 crystals. The code outputs the energy and time of the energy depositions (time is used to properly take into account correlations in nuclear decay chains). A second program takes the output of MCuoreZ and applies a detector response function and incorporates other read-out features.  <2·10 −9 <9·10 −9 <1·10 −6 Si heaters [23,24] <3·10

Monte Carlo Simulation
MCuoreZ is implemented in Geant4 version 4.9.6.p03. α, β and γ particles, nuclear recoils, neutrons and muons are propagated down to keV energies, with an optimization done on the different volumes to balance simulation accuracy and speed. We have chosen the Livermore physics list, and particles can be generated and propagated in the bulk and on the surface of all components. The surface contamination is modeled, according to diffusion processes, with an exponential density profile and a variable depth parameter.
Single radioactive decays as well as the 238 U and 232 Th decay chains have been implemented using the G4RadioactiveDecay database. 2νββ is parameterized according to [25]. The generation of external muons, neutrons and γ is described in [18,19]. Due to the low counting statistics in CUORE-0 data some components are indistiguishable (i.e. exhibit degenerate spectral shapes) 1 and can be grouped, provided that the prior on their material contamination is properly evaluated. This simplification reduces the number of free parameters in the final fit of the simulation to the data. A similar case holds for components made of the same material (i.e., characterized by an identical contamination): their simulations can be grouped once scaled by mass or surface to properly equalize the contamination densities. The detector as built and the geometry as implemented in MCuoreZ are shown in Fig. 1. The outermost volume included in CUORE-0 geometry is the external neutron shield. Although its contamination is negligible, its effect on the propagation of external neutrons and muons is important. The next layer is the external lead shield made of modern lead (ExtPb). Only the inner 10 cm are simulated, since self-shielding is so high that the contribution from the outer volume is negligible.
The cryostat is then modeled as composed by two volumes: the Cryostat External Shields (CryoExt) and the Cryostat Internal Shields (CryoInt).
The CryoExt groups three components made of copper that has been underground for more than 25 years. These components exhibit degenerate spectra and the only prior considered for their contamination is that on 60 Co derived from their identical history. Included in this volume are a small amount of superinsulation material (Mylar and Al) and the gap between the OVC (cryostat Outer Vacuum Chamber) and the ExtPb. The ExtPb is enclosed in acrylic glass and flushed with nitrogen (anti-radon box) to prevent radon from filling this gap and entering the experimental volume.
The CryoInt groups three shields made of the same newer copper, underground since 2002: the IVC (cryostat Inner Vacuum Chamber), the 600 mK and 50 mK shields. The roman lead [26] shield (IntPb) is inserted between the IVC and the 600 mK shield. The Holder is the structure that holds the crystals. It is made from NOSV copper, a special copper alloy from Aurubis company suitable for cryogenic use, cleaned according to the CUORE protocols [14]. It has two main parts: the frame that supports the crystals and the surrounding cylindrical box. Finally, the Crystals are designed as TeO 2 cubes with identical contaminations.
In Table 3, the components listed as Small Parts are not included in the following analysis because of their small mass and negligible contamination (see Tables 1  and 2). Only the PTFE spacers could provide a sizable contribution to the background, however their spectra are degenerate with the Holder ones, therefore their contribution is included in the latter element.

Monte Carlo Data Production
In order to make the Monte Carlo reproduce experimental data -as already anticipated -a second code is used to recreate the detector time and energy response. We account for the timing resolution of each crystal by combining energy depositions that occur in the same crystal within a window of ±5 ms. The absolute time of events is assigned based on a random distribution with an event rate of 1 mHz (the experimental rate during physics runs) and any physical coincidences from MCuoreZ are preserved.
Once the simulated events are correlated correctly in time, the resulting energy depositions are smeared with a Gaussian energy response function. The width of the function varies linearly with energy and is based on the FWHM resolution measured on γ peaks in the M spectrum, between 511 and 2615 keV.
We reproduce also the energy dependence of the trigger efficiency. This is experimentally evaluated for each bolometer using heater-generated pulses with variable amplitudes. For each amplitude (then converted into a particle equivalent energy) the efficiency is defined as the number of triggered signals over the generated ones. The efficiency vs. energy curve is interpolated with an Erf function and used in simulation data production.
As is done in the experiment, events with coincidences in multiple bolometers within a ±5 ms window are combined into multiplets. Pile-up events (i.e. events occurring in the same bolometer within the pile-up rejection window set by the analysis -see Section 2) are removed to account for dead time.
The simulation properly reproduce the γ particle energies, since we calibrate CUORE-0 spectra using 232 Th γ lines. However, we observe in our data an energy quenching effect for α particles compared to γ, which makes the measured energies higher than the known energies [27,28]. We account for this effect in the simulations by shifting α energy depositions by 0.8%, which is the average shift observed for 232 Th and 238 U α peaks visible in CUORE-0 background spectrum.

Reconstruction of the Calibration Source
We use calibration data to test the MCuoreZ simulation. The thoriated wire calibration source is deployed between the CryoExt and the ExtPb, therefore its simulation involves all volumes internal to the ExtPb. The calibration data rate is about 100 times higher than in physics runs, so pile-up effects become sizeable. Fig. 2 shows that with proper treatment of pile-up effects, we achieve a good data-simulation agreement. We observe small deviations in the low energy region between 100 and 300 keV and on a few peaks. The deviation in the low energy region (the origin is still unknown) is a potential source of systematic errors that will be considered in Sec. 8. For what concerns peaks, the largest deviation observed is a 7% difference on the 228 Ac line at 1153 keV (branching ratio 0.1%). We traced this effect to the Geant4 version 2 ; anyway it has no consequences on the background reconstruction since it involves only few very low intensity lines.

Construction of Background Model
The really critical part of the background model construction is the determination a priori of the most relevant sources to be included in the fit. Omitting a relevant source could lead to a poor fit or, even worse, a good fit to the wrong model. The degeneracy of many of the sources energy spectra introduces a further complication. The analysis reported in this section uses lines in the M 1 and Σ 2 spectra to identify the radioactive isotopes contributing to the background counting rate, their possible location and any other characteristic feature that needs to be considered in the construction of the background model. Sources such as the external muons, neutrons, and gammas are straight-forward to model since their fluxes are well known from independent laboratory measurements. The radioactive decays are more difficult to disentangle. The location and distribution of the contaminants must be specified, and in the case of decay chains a choice must be made between the secular equilibrium hypothesis and a break in the chain due to material handling. We use the external screening measurements presented in Sec. 3, estimates of cosmogenic exposure, and distinctive features in the data itself to select the most probable sources. These data also provide the priors for the fit. Figure 3 shows the M 1 spectrum. The spectrum below the 2.615 MeV γ line of 208 Tl includes many γ lines. Peaks are grouped in a single bin, while in the continuum the minimum bin size is 15 keV. Bottom panel: Bin-by-bin ratio between counts in the experimental spectrum and counts in the reconstructed one. We observe a less than 10% discrepancy in the rate at low energies (< 300 keV) that is more likely due to small errors in the geometry reconstruction.  Tables 4 and 5. The lines from 40 K and other contaminants are also present; see Table 6. The 232 Th activity is basically unchanged relative to Cuoricino, especially at high energies. This suggests that most of the observed 232 Th contamination is located in the cryostat and lead shields (i.e. the two structures present in both experiments). There is also a contribution from close sources, most likely the copper Holder , since we do not observe a large attenuation of the low-energy γ lines relative to the high-energy ones. Moreover, if the tower data are divided into floors of crystals, we do observe a floor by floor dependence that is not consistent with the one predicted by simulations.
In Monte Carlo simulations this behavior can be reproduced by a point source close to the crystal that recorded the highest γ counting rate. However, the floor by floor dependence is not seen in the α region, so the source cannot be directly facing the detectors, and is most likely located in the 50 mK copper shield. A similar anomaly was present in Cuoricino. The two strongest γ emitters of the 238 U chain are the 222 Rn daughters: 214 Pb and 214 Bi. All of the observed 238 U peaks are attributed to these two isotopes. Compared to Cuoricino, the average 238 U activity is a factor two to three lower in CUORE-0, revealing a better cleaning of the detector components and/or a better 222 Rn control. Indeed, the intensity of the observed lines, particular those from 214 Bi, vary from dataset to dataset by as much as a factor of 3. Since 222 Rn is a gas, it emanates from contaminated materials in the laboratory environment and its air concentration changes depending on ventilation conditions. As described in Section 4, the most likely place for 222 Rn to enter the experimental setup is in the gap between the OVC and the ExtPb. Occasionally, the nitrogen flushing system inside the acrylic glass box malfunctions and the radon level increases in this volume. This additional varying component of the 238 U contamination is well modeled by an increase in the contamination of the CryoExt.
In addition, we observe a line at 803 keV due to 210 Po not in equilibrium with the 238 U chain. This is fully explained by the 210 Pb (τ 1/2 =22.2 y) bulk contamination of the ExtPb. Before it was installed 25 y ago, this lead had a 210 Pb activity of 16 Bq/kg, which is now reduced by a factor ∼2. This same line is also present in Cuoricino data. The activity reported in Table 6 for the 803 keV peak includes a contribution from a 214 Bi line at 806 keV from the 238 U chain, that cannot be disentangled. 40 K is a unique long-lived isotope that decays via both β − and electron capture, with a negligible β + branching. The signature of this contaminant is the single γ line at 1.46 MeV, which is produced in ∼10% of the decays. The single γ line is not sufficient to determine the activity in each volume, however an asymmetry in the rate along the CUORE-0 tower suggests that in addition to bulk contaminants there is a 40 K (extended) source at the bottom of the cryostat. A number of isotopes produced by cosmogenic activation are also identified. The most critical of these contaminants is 60 Co. The coincidence of its two γ lines in a single crystal produces a peak in the 0νββ ROI (at 2505.7 keV). 60 Co is primarily the result of fast neutron interactions on copper [29], but it can also be produced in tellurium [30,31]. Exploiting a coincident γ analysis we can set an upper limit of 3 × 10 −7 Bq/kg for the 60 Co concentration in the crystals, i.e. well below the level needed to explain the observed rate. Therefore most of the 60 Co contamination is located in the copper components.
The main copper components are the Holder , Cryo-Int, and CryoExt. Following [29], we derive upper limits on their 60 Co activities based on the time they spent at sea-level and underground. The time spent underground before the start of data taking was a few months for the Holder , 14 y for the CryoInt, and 25 y for the CryoExt. Assuming the 60 Co to be in saturation (i.e. the condition in which the production rate equals the decay rate) in both CryoInt and CryoExt just before the underground storage, its activity today would be about 180 µBq/kg for CryoInt and 42 µBq/kg for Cry-oExt. These are used as upper limits for the real contamination. We have only a very rough estimate for the cosmic ray exposure of the copper Holder (since this is made of many parts, each with its own history of exposure and underground storage), which results in an expected activity of ∼50 µBq/kg. As expected, we also observe 54 Mn and 57 Co from the activation of copper [32]. Due to their short half-lives, they can only be located in the Holder . The nuclear fallout product 137 Cs is also present as it was in Cuoricino, therefore it is located in the cryostat copper structures (since it can't contaminate ancient Roman lead). Due to the low intensity of the peak CryoInt and CryoExt spectra are degenerate with each other, and we use the first to account for this contamination. The only expected contaminant due to the cosmogenic activation of the crystals themselves is 125 Sb and its daughter 125m Te; in fact, 125 Sb is the only long lived isotope produced in tellurium with high cross-section [30]. There is evidence for cosmogenic activation of the Roman lead in the IntPb. We observe three γs emitted by 108m Ag (see Table 6), a silver isotope produced by the neutron activation of silver impurities in the Roman lead with τ 1/2 = 438 y. We observe a peak at 439 keV that is ascribed to 202 Tl, a daughter of the long-lived cosmogenic activation product 202 Pb (τ 1/2 =52.5 × 10 3 y). The presence of 202 Tl contaminant is confirmed by a gamma spectroscopy measurement of a IntPb sample. Finally, we observe the two high-energy γ peaks at 1063.7 keV and 1770.2 keV from the decay of 207 Bi without the more intense line at 569 keV, suggesting that this fallout product contaminates the ExtPb, as confirmed by MCuoreZ simulations.

CUORE-0 α Region Analysis
The peaks observed in the α region of CUORE-0 ( Fig. 5) come from the 238 U and 232 Th decay chains and 190 Pt. Due to the short range of α particles and recoiling nuclei, the sources are contaminants close to or inside the bolometers: the Holder and the Crystals. The energy, multiplicity and intensity of the peaks can be used to efficiently constrain 238 U and 232 Th activities in these close components, simplifying the reconstruction of the more complicated γ region. Furthermore they allow us to differentiate bulk and surface as well as Holder and Crystals contaminations. Contaminants in the bulk of the Crystals produce M 1 events with Gaussian peaks centered at the Q-value, since both the α and nuclear recoil are detected in the same crystal. M 2 events are a clear indication of the presence of contaminants on the surfaces of Crystals; they are produced when the α or the recoiling nucleus escapes the source crystal to enter one of the neighboring ones. These events reconstruct in the Σ 2 spectrum at the decay Q-value while in the M 2 spectra they produce two peaks each: one at the recoil energy (∼70-100 keV) and one at the α energy ( ∼70-100 keV below Q-value); see Fig. 6. As the contaminants become deeper, the α loses more energy in the material where it originates from and a low-energy tail on the α peaks (or high-energy tail on the recoil peaks) becomes more pronounced. Contaminants on surfaces of Crystals also contribute to the M 1 spectrum when both the α and recoiling nucleus are stopped within the source crystal or the escaping particles are not detected by another bolometer. Contaminants in the Holder produce M 1 events with at most the energy of the α. Again, depending on the depth, the line shape can vary from a peak with a strong low energy tail (shallow depth), to a flat continuum with no noticeable structures and extending far below the energy of the α (deep surface or bulk contamination). In general, the non-ideal behavior of the detector complicates the reconstruction of α spectra. Some effects are well modeled in simulation: the energy threshold depletes M 2 in favor of M 1 , and αs are quenched with respect to γs. Some effects are not fully understood, e.g. the broadening of all the peaks produced by surface contaminants, possibly explained by a different response of the bolometer to surface and bulk energy depositions. We choose a variable range binning to minimize these uncertainties as described in Section 7. Table 7 lists the isotopes producing the most prominent α peaks identified in CUORE-0 spectra. Surface contaminants in the Crystals are identified on the basis of the Σ 2 spectrum in Fig. 5. All the visible peaks belong to the 232 Th or 238 U chains. We observe breaks in secular equilibrium, particularly evident in the case of the 5.   [22] due to the chemical similarities of tellurium and polonium. 232 Th and 230 Th (a 238 U chain product) produce sharp lines that can be explained by a bulk contamination of the Crystals. Most likely, these two long-lived isotopes are the only survivors of a 232 Th and 238 U contamination in the materials used for crystal production due to chemical affinity of thorium and oxygen. Finally, we find evidence for a bulk 190 Pt contamination of Crystals. The line in the CUORE-0 spectrum is shifted 15 keV above the isotope Q-value of 3249 keV, a much higher shift than the one observed for all the other lines in the α region. The crystals are grown in a platinum crucible. This contamination can be introduced during the crystal growth as a small grain of platinum locally modifying the detector response.

Source list
A list of the 57 sources used for the present analysis is shown in Table 8. The components in the first column are the 6 defined in the fifth column of Table 3. The priors in the fourth column derive mainly from Table 1. The priors for the 60 Co activities are discussed in Sec. 6.1. There are two components included in the analysis but not included in Table 3. One is the 50 mK Spot, a point-like source located on the internal surface of the 50 mK thermal shield facing tower floor 10 (see Fig 1). The other is the Bottom Plate ExtPb, a disc-like source placed on the internal bottom plate of the external lead shield ExtPb. These two sources model the 232 Th excess observed on the floor number 10 and the 40 K excess observed on bottom floors, which are discussed in Sec 6.1. In order to properly reproduce the shape of the α peaks in M 1 and M 2 spectra, Crystals and Holder surface contamination with a few representative depths are included in the analysis (Fig. 7). The exact choice of contaminant depth is treated as a systematic uncertainty. This is discussed in Sec. 8. Violations of secular equilibrium in the 238 U and 232 Th decay chains are only considered when they produce a distinguishable signature in the spectrum. A background contribution due to muons is included in the analysis; contributions from neutrons and external γs are negligible. Muons contribute mainly through γs produced by interactions in the detector components. A prior is obtained from high multiplicity data, which is compatible with the muon flux measured by other experiments.

Bayesian Fit Construction
The activities of the sources used for the background model are determined by fitting the observed M 1 , M 2 , and Σ 2 spectra with a linear combination of simulated source spectra. The expectation value of the counts in the i-th bin of the experimental spectrum is given by: where C M C ij,α is the expectation value for the i-th bin of the simulated spectrum for the j-th source and N j is the unknown activity of the j-th source. To ensure sufficient statistics, both the experimental and simulated spectra are summed over all the active crystals. The fit is performed with a Bayesian approach using JAGS (Just Another Gibbs Sampler) [33,34,35]. JAGS exploits Markov Chain Monte Carlo simulations to sample the joint posterior probability distribution function (pdf) of the model parameters. Following Bayes' theorem the posterior pdf is evaluated combining the likelihood and the prior distributions. The available data to define the likelihood s are the observed counts in the bins of the experimental (C exp i,α ) and simulated (C M C ij,α ) spectra, both of which obey Poisson statistics. Therefore, the joint posterior pdf is: The prior s for N j , which describe our prior knowledge about source activities, are specified in Table 8. In the case of a measured activity, we adopt a Gaussian prior centered at the measured value with the measurement uncertainty as the width of the gaussian. For upper limits, we adopt a half gaussian with a width such that our 90% upper limit is the 90% value of the prior. In all the other cases, we use a uniform non-informative prior with an activity that ranges from 0 to an upper limit higher than the maximum activity compatible with the CUORE-0 data. Similarly, we use uniform prior s over wide ranges for C M C ij,α . We chose a variable binning of the spectra to maximize the information content while minimizing the effects of statistical fluctuations and detector non-ideal behavior. Therefore, to avoid systematic uncertainties due to the lineshape, all the counts belonging to the same γ or α peak are included in a single bin. The minimum bin size in the continuum is 15 keV, and bins with less than 30 counts are merged with their immediate neighbor. The fit extends from 118 keV to 7 MeV. The threshold at 118 keV is set to exclude the low-energy noise events (contaminating few datasets) and the nuclear recoil peak (that is mis-calibrated). In building the Σ 2 spectrum, we require that the energy of each event is above threshold. An exception is set for events with E > 2.7 MeV in coincidence with events below the fit threshold, to correctly build-up the Q-value peaks in the α region of Σ 2 spectrum.

Reference Fit and Systematics
The reference fit is the result of the fit to data from the total 33.4 kg·yr TeO 2 exposure. The reconstructions of the experimental spectra are shown in Figs. 8, 9, and 10 for the M 1 , M 2 , and Σ 2 spectra, respectively. The fit results for the 57 free parameters are summarized in Table 8. The marginalized posterior distributions are used to evaluate the central values and the statistical uncertainties of the activities, or to calculate 90% upper limits for undetermined contaminations. The correlation matrix of the 57 sources is illustrated in Fig. 11. In general, most of the components used to fit the α region are not correlated to those used to reconstruct the γ region. As expected, the same contaminants in neighboring components of the experimental setup are highly correlated due to the similar spectra. For the M 1 , M 2 , and Σ 2 spectra, the normalized fit residuals have a Gaussian-like distribution with mean 0 and standard deviation 1. The reduced chi-square with 57 parameters and 478 degrees of freedom is 1.36. We do not expect perfect statistical agreement between the data and reconstruction since the uncertainties associated with the simulated spectra account for only the statistical fluctuation in the bin counts. They do not include the systematic uncertainties in the MC simulations.
To check the stability of the background model, the dependence on priors, and the systematic uncertainties, especially those affecting the 2νββ half-life, we run Table 8 List of the sources used to fit the CUORE-0 background data. The columns show (1) the name of the contaminated element, (2) the source index (j in eq. 1), and (3) the contaminant. If not otherwise specified, 232 Th, 238 U, and 210 Pb refer to the whole decay chains in secular equilibrium, while the label "only" indicates that only the decay of the specified isotope is generated. For surface contaminants, the simulated depth is indicated in µm. Column (4) reports the prior used in the fit, when not specified a non informative prior is used (see text for details). Column (5) reports the posterior with the statistical error (limits are 90% C.L.). Column (6) reports the range of systematic uncertainties (limits are 90% C.L.). In the case of crystal sources, systematic uncertainties can arise from non-uniform contaminants in the different crystals.

Component
Index    Table 8. Note the anti-correlation between the 2νββ and the 40 K activity in Crystal bulk. 2νββ activity changes by less than 1%. These tests also cover the systematics due to miscalibration [12].
-Fit Energy Threshold. We run the fit with different energy thresholds ranging from 118 to 500 keV, covering the region where the reconstruction of the calibration source is slightly worse (Sec. 5). The quality of fit reconstruction is unchanged and the 2νββ activity variations are below 2%.
-Contamination Depth Uncertainty. We fit the α region with different depths to model surface contaminations. Several models perform similarly to the reference fit, however the results of the background reconstruction, particularly the 2νββ rate, are unaffected.
-Dependence on Prior distributions. To evaluate the systematic uncertainty related to the Prior choice, we perform two JAGS fits. In the first fit, the half-Gaussian priors used in the case of upper limits on source activities are changed to uniform priors with the minimum at 0 and the maximum at 3σ above the upper limit. In the second fit, uniform non-informative priors are used for all components. In both cases, the global fit reconstruction is good and the 2νββ result changes by ∼1%.
-Selection of background sources. In the reference fit there are 14 undetermined sources whose activity is quoted as upper limit. To check the fit stability against the removal of these sources, we run a minimum model fit with only 43 sources. Once more, the global fit reconstruction and the 2νββ result are not affected. In the tests detailed above, the overall goodness of the fit remains stable, while we observe variations in the activities of the individual sources. These variations are used as an evaluation of the systematic uncertainty on the 57 source activities (Table 8, sixth column). There are caveats using the reference fit results as an exact estimation of the material contamination. Indeed, degenerate source spectra allow us to use a single source to represent a group of possible sources. Examples are: the Holder that also accounts for the contribution of the Small Parts, surface contaminants in close components that are modeled with few representative depths, or bulk contamination in far components that also include surface ones.
9 130 Te 2νββ decay The background reconstruction allows us to measure the 2νββ of 130 Te with high accuracy. Fig. 12 shows the fit result compared with the CUORE-0 M 1 . 2νββ produces (3.27 ± 0.08) × 10 4 counts, corresponding to ∼10% of the events in the M 1 γ region from 118 keV to 2.7 MeV. As shown in Fig. 13, removing the 2νββ component results in a dramatically poorer fit in this region. The 2νββ activity is (3.43 ± 0.09) × 10 −5 Bq/kg, with a statistical uncertainty that is amplified by the strong anti-correlation to the 40 K contamination in crystal bulk (but not to other 40 K sources). Indeed, this is the only case where the β spectrum of 40 K (having a shape that resembles that of 2νββ) contributes to the detector counting rate. For all the other 40 K sources, only the EC decay (branching ratio 89%) contributes to the detector counting rate through the 1460 keV line and its Compton tail. The Posterior for the 2νββ activity as obtained from the reference fit is shown in Fig 14. Also shown is the Posterior associated to the fit bias. This is derived from systematic studies discussed in Sec. 8 and is represented as a flat distribution. Fig 14 also shows the 68% Confidence Intervals associated to statistical and systematic errors.
The half-life value obtained for 2νββ is  The signature of 130 Te 0νββ is a Gaussian line centered at 2528 keV, the transition energy of the isotope, in M 1 . Modeling the shape of the background in this region, especially possible subdominant peaks, and identifying the main sources of background is relevant not only for CUORE-0, but also for the future evolution of 0νββ searches with bolometers. It is useful to group the sources used for the fit into three major classes: the two elements that will be identical (though replicated 19 times) in CUORE, Holder  and Crystals, and the element that will change, i.e. the cryogenic and radioactive shield systems (the sum of the CryoInt, CryoExt, IntPb, and ExtPb). The contribution from these elements to the 0νββ region of interest (2470-2570 keV) in M 1 are shown in Fig. 15 and listed in Table 9. The largest contribution comes from the shields. This is mainly 232 Th contamination. The Holder is the second largest contributor due to degraded αs from 238 U and 232 Th deep surface contaminants. Bulk and shallow-depth contaminants account for less than 0.3% of the background. A very small fraction of the background comes from 238 U, 232 Th, and 210 Pb Crystals surface contaminants, and from muon interactions. The systematic uncertainties are negligible.
Finally, Fig. 16 shows a wider region centered around the ROI. This plot is produced by tagging the energy depositions where at least 90% of the energy was deposited by α particles. We found ∼ 24% of the ROI background was produced by α events. After reducing γ backgrounds from the shields, these α events are ex-  Table 9 Sources contributing to the 0νββ ROI. The flat counting rate in this region (i.e. excluding the 60 Co sum peak) is 0.058±0.006 counts/(keV kg y) [5]. Column (2) reports the contribution of the different sources. "Shields" here stands for the sum of CryoInt, CryoExt, IntPb, and ExtPb. pected to dominate the ROI rate in CUORE. This motivates the development of α particle discrimination for future bolometer-based experiments; see [36] and references therein.

Conclusion
In this paper, we successfully reconstruct the CUORE-0 background using 57 sources modeled using a detailed Monte Carlo simulation. We find that 10% of the M 1 counting rate in the range [118−2700] keV is unequivocally due to 130 Te 2νββ decay. We measure its half-life to be: T 2ν 1/2 = [8.2 ± 0.2 (stat.) ± 0.6 (syst.)] × 10 20 y. Compared to previous results obtained from MiDBD [6.1 ± 1.4 (stat.) +2.9 −3.5 (syst.)] × 10 20 y [37], from Cuoricino [38], and from NEMO [7.0 ± 0.9 (stat.) ± 1.1 (syst.)] × 10 20 y [39], this is the most precise measurement to date. We find that the background rate in the 130 Te 0νββ region of interest is dominated by the shields. This result gives us confidence that we are on track to achieve the requirements for 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 was supported by the Istituto Nazionale di