∎

11institutetext:  Department of Physics and Astronomy, University of South
Carolina, Columbia, SC 29208, USA 22institutetext: Department of Physics and
Astronomy, University of California, Los Angeles, CA 90095, USA
33institutetext: INFN – Laboratori Nazionali di Legnaro, Legnaro (Padova)
I-35020, Italy 44institutetext: INFN – Sezione di Bologna, Bologna I-40127,
Italy 55institutetext: Dipartimento di Fisica, Sapienza Università di Roma,
Roma I-00185, Italy 66institutetext: INFN – Sezione di Roma, Roma I-00185,
Italy 77institutetext: INFN – Laboratori Nazionali del Gran Sasso, Assergi
(L’Aquila) I-67100, Italy 88institutetext: INFN – Sezione di Milano Bicocca,
Milano I-20126, Italy 99institutetext: Dipartimento di Fisica, Università di
Milano-Bicocca, Milano I-20126, Italy 1010institutetext: Center for Neutrino
Physics, Virginia Polytechnic Institute and State University, Blacksburg,
Virginia 24061, USA 1111institutetext: INFN – Sezione di Genova, Genova
I-16146, Italy 1212institutetext: Dipartimento di Fisica, Università di
Genova, Genova I-16146, Italy 1313institutetext: Massachusetts Institute of
Technology, Cambridge, MA 02139, USA 1414institutetext: Key Laboratory of
Nuclear Physics and Ion-beam Application (MOE), Institute of Modern Physics,
Fudan University, Shanghai 200433, China 1515institutetext: Department of
Physics, University of California, Berkeley, CA 94720, USA 1616institutetext:
Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA
94720, USA 1717institutetext: INFN – Gran Sasso Science Institute, L’Aquila
I-67100, Italy 1818institutetext: Wright Laboratory, Department of Physics,
Yale University, New Haven, CT 06520, USA 1919institutetext: INFN – Laboratori
Nazionali di Frascati, Frascati (Roma) I-00044, Italy 2020institutetext:
Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France
2121institutetext: Physics Department, California Polytechnic State
University, San Luis Obispo, CA 93407, USA 2222institutetext: INPAC and School
of Physics and Astronomy, Shanghai Jiao Tong University; Shanghai Laboratory
for Particle Physics and Cosmology, Shanghai 200240, China 2323institutetext:
Dipartimento di Fisica e Astronomia, Alma Mater Studiorum – Università di
Bologna, Bologna I-40127, Italy 2424institutetext: Service de Physique des
Particules, CEA / Saclay, 91191 Gif-sur-Yvette, France 2525institutetext:
Lawrence Livermore National Laboratory, Livermore, CA 94550, USA
2626institutetext: Department of Nuclear Engineering, University of
California, Berkeley, CA 94720, USA 2727institutetext: Dipartimento di
Ingegneria Civile e Meccanica, Università degli Studi di Cassino e del Lazio
Meridionale, Cassino I-03043, Italy 2828institutetext: Department of Physics
and Astronomy, The Johns Hopkins University, 3400 North Charles Street
Baltimore, MD, 21211 2929institutetext: INFN – Sezione di Padova, Padova
I-35131, Italy 3030institutetext: Engineering Division, Lawrence Berkeley
National Laboratory, Berkeley, CA 94720, USA

fn1Deceased

# Search for Double-Beta Decay of $\mathrm{{}^{130}Te}$ to the $0^{+}$ States
of $\mathrm{{}^{130}Xe}$ with CUORE

D. Q. AdamsUSC  C. AlduinoUSC  K. AlfonsoUCLA  F. T. Avignone IIIUSC  O.
AzzoliniINFNLegnaro  G. BariINFNBologna  F. BelliniRoma,INFNRoma  G.
BenatoLNGS  M. BiassoniINFNMiB  A. BrancaMilano,INFNMiB  C.
BrofferioMilano,INFNMiB  C. BucciLNGS  J. CamilleriVirginiaTech  A.
CaminataINFNGenova  A. CampaniGenova,INFNGenova  L. CanonicaMIT,LNGS  X. G.
CaoFudan  S. CapelliMilano,INFNMiB  L. CappelliLNGS,BerkeleyPhys,LBNLNucSci
L. CardaniINFNRoma  P. CarnitiMilano,INFNMiB  N. CasaliINFNRoma  E.
CeliGSSI,LNGS  D. ChiesaMilano,INFNMiB  M. ClemenzaMilano,INFNMiB  S.
CopelloGenova,INFNGenova  C. CosmelliRoma,INFNRoma  O. CremonesiINFNMiB  R. J.
CreswickUSC  A. D’AddabboGSSI,LNGS  I. DafineiINFNRoma  C. J. DavisYale  S.
Dell’OroMilano,INFNMiB  S. Di DomizioGenova,INFNGenova  V. DompèGSSI,LNGS  D.
Q. FangFudan  G. FantiniRoma,INFNRoma  M. FaverzaniMilano,INFNMiB  E.
FerriMilano,INFNMiB  F. FerroniGSSI,INFNRoma  E. FioriniINFNMiB,Milano  M. A.
FranceschiINFNFrascati  S. J. FreedmanLBNLNucSci,BerkeleyPhys,fn1  S.H.
FuFudan  B. K. FujikawaLBNLNucSci  A. GiacheroMilano,INFNMiB  L.
GironiMilano,INFNMiB  A. GiulianiParis-Saclay  P. GorlaLNGS  C. GottiINFNMiB
T. D. GutierrezCalPoly  K. HanSJTU  K. M. HeegerYale  R. G. HuangBerkeleyPhys
H. Z. HuangUCLA  J. JohnstonMIT  G. KeppelINFNLegnaro  Yu. G.
KolomenskyBerkeleyPhys,LBNLNucSci  C. LigiINFNFrascati  L. MaUCLA  Y. G.
MaFudan  L. MariniBerkeleyPhys,LBNLNucSci  R. H. MaruyamaYale  D. MayerMIT  Y.
MeiLBNLNucSci  N. MoggiBolognaAstro,INFNBologna  S. MorgantiINFNRoma  T.
NapolitanoINFNFrascati  M. NastasiMilano,INFNMiB  J. NikkelYale  C.
NonesSaclay  E. B. NormanLLNL,BerkeleyNucEng  A. NucciottiMilano,INFNMiB  I.
NutiniMilano,INFNMiB  T. O’DonnellVirginiaTech  J. L. OuelletMIT  S. PaganYale
C. E. PagliaroneLNGS,Cassino  L. PagnaniniGSSI,LNGS  M.
PallaviciniGenova,INFNGenova  L. PattavinaLNGS  M. PavanMilano,INFNMiB  G.
PessinaINFNMiB  V. PettinacciINFNRoma  C. PiraINFNLegnaro  S. PirroLNGS  S.
PozziMilano,INFNMiB  E. PrevitaliMilano,INFNMiB  A. PuiuGSSI,LNGS  C.
RosenfeldUSC  C. RusconiUSC,LNGS  M. SakaiBerkeleyPhys  S. SangiorgioLLNL  B.
SchmidtLBNLNucSci  N. D. ScielzoLLNL  V. SharmaVirginiaTech  V.
SinghBerkeleyPhys  M. SistiINFNMiB  D. SpellerJHU  P.T. SurukuchiYale  L.
TaffarelloINFNPadova  F. TerranovaMilano,INFNMiB  C. TomeiINFNRoma  K. J.
VetterBerkeleyPhys,LBNLNucSci  M. VignatiINFNRoma  S. L.
WagaarachchiBerkeleyPhys,LBNLNucSci  B. S. WangLLNL,BerkeleyNucEng  B.
WelliverLBNLNucSci  J. WilsonUSC  K. WilsonUSC  L. A. WinslowMIT  S.
ZimmermannLBNLEngineering  S. ZucchelliBolognaAstro,INFNBologna

(Received: date / Accepted: date)

###### Abstract

The CUORE experiment is a large bolometric array searching for the lepton
number violating neutrino-less double beta decay ($0\nu\beta\beta$) in the
isotope $\mathrm{{}^{130}Te}$. In this work we present the latest results on
two searches for the double beta decay (DBD) of $\mathrm{{}^{130}Te}$ to the
first $0^{+}_{2}$ excited state of $\mathrm{{}^{130}Xe}$: the $0\nu\beta\beta$
decay and the Standard Model-allowed two-neutrinos double beta decay
($2\nu\beta\beta$). Both searches are based on a 372.5 kg$\times$yr TeO2
exposure. The de-excitation gamma rays emitted by the excited Xe nucleus in
the final state yield a unique signature, which can be searched for with low
background by studying coincident events in two or more bolometers. The
closely packed arrangement of the CUORE crystals constitutes a significant
advantage in this regard. The median limit setting sensitivities at 90%
Credible Interval (C.I.) of the given searches were estimated as
$\mathrm{S^{0\nu}_{1/2}=5.6\times 10^{24}\>\mathrm{yr}}$ for the
${0\nu\beta\beta}$ decay and $\mathrm{S^{2\nu}_{1/2}=2.1\times
10^{24}\>\mathrm{yr}}$ for the ${2\nu\beta\beta}$ decay. No significant
evidence for either of the decay modes was observed and a Bayesian lower bound
at $90\%$ C.I. on the decay half lives is obtained as:
$\mathrm{(T_{1/2})^{0\nu}_{0^{+}_{2}}>5.9\times 10^{24}\>\mathrm{yr}}$ for the
$0\nu\beta\beta$ mode and $\mathrm{(T_{1/2})^{2\nu}_{0^{+}_{2}}>1.3\times
10^{24}\>\mathrm{yr}}$ for the $2\nu\beta\beta$ mode. These represent the most
stringent limits on the DBD of 130Te to excited states and improve by a factor
$\sim 5$ the previous results on this process.

††journal: Eur. Phys. J. C

##  1 Introduction

Double beta decay (DBD) is an extremely rare nuclear process where a
simultaneous transmutation of a pair of neutrons into protons converts a
nucleus (A, Z) into an isobar (A, Z+2), with the emission of two electrons and
two anti-neutrinos. This two-neutrino decay mode ($2\nu\beta\beta$) is
predicted in the Standard Model and was detected in several nuclei. The
neutrinoless mode of the decay ($0\nu\beta\beta$) is a posited Beyond Standard
Model process that could shed light on many open aspects of modern particle
physics and cosmology such as the existence of lepton number violation and
elementary Majorana fermions, the neutrino mass scale, and the baryon
asymmetry in the Universe Racah:1937qq ; Furry:1939qr ; Pontecorvo:1968 ;
Schechter:1981bd ; Fukugita:1986hr . Both DBD modes can proceed through
transitions to the ground state as well as to various excited states of the
daughter nucleus. While the former can be easier to detect through their
shorter half-lives, the latter leaves a unique signature which may be detected
with significantly reduced backgrounds. The excited state decays also provide
powerful tests of the nuclear physics of DBD and can shed light on nuclear
matrix element calculations as well as the ongoing discussion on the quenching
of the effective axial coupling constant $g_{A}$; eventually, they could even
be used to disentangle the mechanism of $0\nu\beta\beta$ decay Avignone:2007fu
.

So far, $2\nu\beta\beta$ decay to the first $0^{+}$ excited state has been
observed in only 2 isotopes: ${}^{100}\mathrm{Mo}$ Barabash:1995fn  and
${}^{150}\mathrm{Nd}$ Barabash:2004dv , with half lives of
$(\mathrm{T}_{1/2})^{2\nu}_{0^{+}}=6.1^{+1.8}_{-1.1}\times
10^{20}\>\mathrm{yr}$ and
$(\mathrm{T}_{1/2})^{2\nu}_{0^{+}}=1.4^{+0.4}_{-0.2}\mathrm{(stat.)}\pm
0.3\mathrm{(syst.)}\times 10^{20}\>\mathrm{yr}$, respectively. Searches for
the same process in other isotopes has yielded lower limits from $3.1\times
10^{20}\>\mathrm{yr}$ to $8.3\times 10^{23}\>\mathrm{yr}$ at $90\>\%$
Confidence Level (C.L.) (see Ref. Barabash:2017bgb  for a review).

In this work, we focus on the search for $0\nu\beta\beta$ and $2\nu\beta\beta$
decays of $\mathrm{{}^{130}Te}$ to the first $0^{+}$ excited state of
$\mathrm{{}^{130}Xe}$ with the CUORE experiment. Presently, the strongest
limits on the decay to excited states half-life of ${}^{130}\mathrm{Te}$ come
from a combination of Cuoricino Andreotti:2011in  and CUORE-0 Adams:2018yrj
data: the latter (not included in Ref. Barabash:2017bgb ) was published
recently and includes the combination of the predecessor’s results. The
obtained limits are:

$\displaystyle(\mathrm{T}_{1/2})_{0^{+}_{2}}^{0\nu}$ $\displaystyle>1.4\times
10^{24}\>\mathrm{yr}\;(90\%\>\mathrm{C.L.})$ (1a)
$\displaystyle(\mathrm{T}_{1/2})_{0^{+}_{2}}^{2\nu}$ $\displaystyle>2.5\times
10^{23}\>\mathrm{yr}\;(90\%\>\mathrm{C.L.})$ (1b)

for the $0\nu$ and $2\nu$ process, respectively. Theoretical predictions
PhysRevC.91.054309  on the $\mathrm{(T_{1/2})}^{2\nu}_{0^{+}_{2}}$ observable
in the $2\nu\beta\beta$ decay channel are based on the QRPA approach and favor
the following range:

${}^{th}(\mathrm{T}_{1/2})^{2\nu}_{0^{+}_{2}}=(7.2-16)\times
10^{24}\>\mathrm{yr}$ (2)

where the range depends on the precise treatment of $g_{A}$. The lower bound
assumes a constant function of the mass number A, and the upper bound assumes
a value of $g_{A}=0.6$ PhysRevC.91.054309 ; SUHONEN2013153 ; SUHONEN20141 .

A data driven estimate of the $2\nu\beta\beta$ ground state to excited state
decay rate in the IBM-II framework based on Ref. Barea:2013bz ; Kotila:2012zza
; Barabash:2010ie  is reported in Ref. Lehnert2015  as

${}^{th}(\mathrm{T}_{1/2})^{2\nu}_{0^{+}_{2}}=2.2\times 10^{25}\>\mathrm{yr}.$
(3)

In this regard, as stated before, both a measurement or a more stringent limit
with respect to Ref. Adams:2018yrj  are informative from the point of view of
refining and validating the theoretical computations.

The decay to excited states has a unique signature. The double-beta decay
emits two electrons, which share kinetic energy up to 734 keV. The subsequent
decay of the excited daughter nucleus typically emits two or three high energy
gamma rays in cascade. Due to the emission of such coincident de-excitation
$\gamma$ rays, both $0\nu\beta\beta$ and $2\nu\beta\beta$ decay channels allow
a significant background reduction with respect to the corresponding
transitions to the ground state. This holds especially in an experimental
setup that exploits a high detector granularity, such as the CUORE experiment.

##  2 Detector and Data Production

The Cryogenic Underground Observatory for Rare Events (CUORE) Arnaboldi:2002du
; Brofferio:2019yoc  is a ton-scale cryogenic detector located at the
underground Laboratori Nazionali del Gran Sasso (LNGS) in Italy. CUORE is
designed to search for the $0\nu\beta\beta$ decay of $\mathrm{{}^{130}Te}$ to
the ground state of $\mathrm{{}^{130}Xe}$ Alduino:2017ehq ; Adams:2019jhp ,
and has a low background rate near the $0\nu\beta\beta$ decay region of
interest (ROI), an excellent energy resolution, and a high detection
efficiency. The CUORE detector consists of a close-packed array of 988 TeO2
crystals operating as cryogenic bolometers Fiorini:1983yj ; Enss:2008x ;
Arnaboldi:2010fj  at a temperature of $\sim$10 mK. The CUORE crystals are
$5\times 5\times 5$ cm3 cubes weighing 750 g each, arranged in 19 towers: each
consisting of a copper structure with 13 floors and 4 crystals per floor. A
custom-made 3He/4He dilution refrigerator, which represents the state of the
art for this cryogenic technique, is used to cool down the CUORE cryostat,
where the entire array is contained and shielded. Alduino:2019xia ;
Alessandria:2013ufa ; Buccheri:2014bma ; Arnaboldi:2017aek ; DiDomizio:2018ldc
; Benato:2017kdf ; DAddabbo:2017efe .

Each CUORE crystal records thermal pulses via a neutron-transmutation doped
(NTD) germanium thermistor Haller:1984ntd  glued to its surface. Any energy
deposition in the crystal causes a sudden rise in temperature and can indicate
the emission of a particle inside, the crossing of a particle through, or some
environmental thermal instability (e.g. earthquakes).

The data acquisition and production of CUORE event data used in this work
closely follows the procedure used in Alduino:2019xia  and is described in
detail in Alduino:2016zrl . We briefly review the basic process here and
highlight the differences.

The NTD converts the thermal signal to a voltage output, which is amplified,
filtered through a 6-pole Bessel anti-aliasing filter, and sampled
continuously at 1kHz. The data are stored to disk and triggered offline with
an algorithm based on the optimum filter (OF) Gatti:1986cw ; DiDomizio:2010ph
; campani2020lowering .

For each triggered pulse, a 10 second window around each trigger (3 seconds
before and 7 seconds after) is processed through a series of analysis steps,
with the aim of extracting the physical quantities associated to the pulse.
The waveform is filtered using an OF built from a pulse template, and the
measured noise power spectrum.

The signal amplitudes are then evaluated from the OF filtered waveforms and
those amplitudes are corrected for small changes in the detector gain due to
temperature drifts. We calibrate each bolometer individually using dedicated
calibration runs with $\mathrm{{}^{232}Th}$ and $\mathrm{{}^{60}Co}$ gamma
sources deployed around the detector array. These calibration runs typically
last a few days every two months.

We impose a pulse shape selection (PSA) based on 6 pulse shape parameters.
This cut removes noisy events, pileup events, and non-physical events.

Unlike the decay to ground state search described in Ref. Alduino:2019xia ,
the physics search described in the present work focuses on coincident energy
depositions in multiple crystals. In particular, we are focusing on events
where energy is deposited in either two or three bolometers. As the
reconstructed time difference between events on nearby bolometers is affected
by differences in pulse rise times, a bolometer-by-bolometer correction is
applied. Sets of coincident energy releases in $M$ bolometers within a $\pm
5\>$ ms time window are grouped together as multiplets of multiplicity $M$.

CUORE started its data taking in May 2017 and, after two significant
interruptions for important maintenance of the cryogenic system, is now seeing
its exposure grow at an average rate of $\sim$ 50 kg$\times$yr/month. The
CUORE data collection is organized in datasets: a dataset begins with a gamma
calibration campaign that typically lasts 2-3 days, followed by 6-8 weeks of
uninterrupted background data taking, and ends with another gamma calibration.

Recently, the CUORE collaboration released the results of the search for
$0\nu\beta\beta$ decay to the g.s. on the the accumulated exposure of 372.5
kg$\cdot$y, setting an improved limit on the half-life of
$\mathrm{{}^{130}Te}$ of $(\rm{T_{1/2}})^{0\nu}_{0^{+}_{1}}>3.2\times 10^{25}$
yr Adams:2019jhp .

##  3 Analysis

In this section we describe the analysis steps that are specific to the search
for ${}^{130}\mathrm{Te}$ decay to the excited states of
${}^{130}\mathrm{Xe}$. The de-excitation of the ${}^{130}\mathrm{Xe}$ nucleus
follows one of three possible patterns, i.e. paths through states of
decreasing energy from the $0^{+}_{2}$ to the $0^{+}_{1}$ ground state (Figure
1). Details about the probability of each de-excitation pattern, referred in
the following as A, B and C (in decreasing order of probability), and the
energy of the emitted $\gamma$ rays are reported in Table 1.

The simultaneous emission of DBD betas and de-excitation gammas produces
coincidence multiplets, i.e. sets of simultaneous pulses in $M$ bolometers,
grouped by the coincidence algorithm. We search for events with full
containment of the final state gammas in the crystals: more specifically we
try to avoid multiplets where one or more of the final state $\gamma$s escape
the source crystal and are absorbed by some non-active part of the
experimental apparatus, or Compton scattering events, where the energy of a
single de-excitation gamma is split among two or more detectors. We place
energy selection cuts to find these events, which are listed in Table 3 and
described in more details in Sec. 3.4.

Partitions are defined as unique groupings of energy depositions that pass a
particular set of energy selection cuts. For a fixed multiplicity $M$ and a
source pattern, they are identified by all possible ways of partitioning the
final state particles in $M$ different crystals. Finally we define signatures
as partitions from different patterns that are indistinguishable. Single-site
($M=1$) signatures are not taken into account, as the $0\nu\beta\beta$ decay
channel would be indistinguishable from the same decay on the ground state of
${}^{130}\mathrm{Xe}$, while the $2\nu\beta\beta$ decay channel, instead,
would suffer from high background from the decay to the ground state.
Therefore, there remain $8$ partitions for patterns A and B, and $14$ for
pattern C. Each of the partitions is labelled with strings of $3$ characters
with the following convention

[Multiplicity] [Pattern] [Index]

where Multiplicity = 2,3,4 indicates the number of involved crystals, Pattern
= A,B,C stores the originating de-excitation pattern, and Index is a unique
integer counter to distinguish the various combinations of energy groupings
for that pattern and multiplicity. Partitions sharing the same expected energy
release are indistinguishable and are merged as signatures. For this reason,
instead of handling a total number of $22$ partitions we are left with $15$
signatures excitedTAUP19 . An example of indistinguishable partitions is given
in Table 3 by the 2A0-2B1 signature. In one of the crystals a gamma energy
release of 1257 keV is expected. This can be either due (see Table 1) to
$\gamma_{1}$ from pattern A or the simultaneous absorption of
$\gamma_{1}+\gamma_{2}$ from pattern B.

Figure 1: The decay scheme of ${}^{130}\mathrm{Te}$ is shown with details about the involved excited states of ${}^{130}\mathrm{Xe}$ up to its first $0^{+}$ excited state. The nomenclature $0^{+}_{1},...,0^{+}_{n}$ indicates states with the same angular momentum in increasing order of excitation energy. An energy scale is shown (right) where the ${}^{130}\mathrm{Xe}$ ground state is taken as reference SINGH200133 . Pattern | BR [%] | Energy $\gamma_{1}$ | Energy $\gamma_{2}$ | Energy $\gamma_{3}$
---|---|---|---|---
A | $86\%$ | $1257\>\mathrm{keV}$ | $536\>\mathrm{keV}$ | -
B | $12\%$ | $\>\>671\>\mathrm{keV}$ | $586\>\mathrm{keV}$ | $536\>\mathrm{keV}$
C | $2\%$ | $1122\>\mathrm{keV}$ | $671\>\mathrm{keV}$ | -
Table 1: The de-excitation $\gamma$ rays emitted by ${}^{130}\mathrm{Xe}^{*}$
in the transition from the $0^{+}_{2}$ to the ground state. Each row
corresponds to a different path through intermediate states. The energies of
the emitted $\gamma$s are listed, in order of energy, along with the branching
ratio (BR) of each pattern SINGH200133 .

###  3.1 Monte Carlo Simulations

We use Monte Carlo (MC) simulations to compute the detection efficiency (Sec.
3.2) and the expected background (Sec. 3.3) associated with each signature, to
rank experimental signatures and eventually fine tune selection cuts on the
most sensitive ones (Sec. 3.4).

CUORE uses a Geant4-based MC to simulate energy depositions in the detector.
The Geant4 software G4_AGOSTINELLI2003250  simulates particle interactions in
the various volumes and materials of a modeled detector geometry. A separate
post-processing step converts the resulting energy depositions into an output
as close to the output of the data production as possible. We refer to Ref.
Alduino_2017  for further details about the CUORE MC simulations.

Signal simulations, that are simulations of the double beta decay to excited
states and the subsequent de-excitation gammas, are produced separately for
each process ($0\nu$,$2\nu$) and pattern (A,B,C). Gamma energies are generated
as monochromatic. Angular correlations induce a negligible effect on the
containment efficiency of the experimental signatures listed in Table 3 as
opposed to isotropic gamma emission and compared with the dominant systematic
uncertainty described in Sec. 5.1. Beta energies are randomly extracted from
the beta spectrum of the corresponding decay Haxton:1985am ; Doi:1981mi ;
Doi:1981mj ; Kotila:2012zza  in the HSD hypothesis 111The $2\nu\beta\beta$ is
called single state dominated (SSD) if it is governed by the lowest $1^{+}$
energy level. In the higher-state dominated (HSD) case the calculation is
simplified by summing over all the virtual intermediate states and assuming an
average closure energy. In the SSD hypothesis the cumulative sum energy
distribution of emitted electrons Kotila:2012zza  differs by $<0.3\%$ with
respect to HSD. We note though that this analysis cannot infer the shape of
the beta spectrum because in $2\nu\beta\beta$ signatures the fit is performed
just on the gamma peaks.

Background simulations take as input the CUORE background model CUORE2nuBB ,
and include contaminations in the crystal and several other parts of the CUORE
setup222The CUORE background model we refer to is still preliminary, however
the estimates of the background activities are good enough to understand what
will be the expected contribution to the present search. In the final fit the
exact values are floated., such as: the copper tower structure, the closest
copper vessel enclosing the detector, the Roman lead, the internal and
external modern lead shields and the internal lead suspension system. The
contaminants include bulk and surface 238U and 232Th chains with different
hypotheses on secular equilibrium breaks, bulk 60Co, 40K, and a few other long
lived isotopes. Additional sources of background included are the cosmic muon
flux and the $2\nu\beta\beta$ decay of 130Te to the ground state. Both signal
and background simulated energy spectra are convolved with a Gaussian
resolution that has a width of $5~{}\mathrm{keV}$ full width at half maximum,
a standard choice for our simulation studies Alduino_2017 .

###  3.2 Efficiency Evaluation

The detection efficiency of a given signature consists of two components: the
containment efficiency and the analysis efficiency. Given a signature $s$ and
a set of energy selection cuts on the involved bolometers, the corresponding
containment efficiency $\varepsilon_{s}$ represents the probability that the
energy released by a nuclear decay of ${}^{130}\mathrm{Te}$ to the $0^{+}_{2}$
state of ${}^{130}\mathrm{Xe}$ matches the topology of the signature. We
evaluate this efficiency component from the signal MC simulations described in
Sec. 3.1, by summing over the contributions of all patterns $p$ populating the
signature $s$

$\varepsilon_{s}=\Bigg{[}\sum_{p}\mathrm{BR}_{p}\cdot\frac{\big{[}N^{(sel)}_{MC}\big{]}^{(s)}_{p}}{\big{[}N^{(tot)}_{MC}\big{]}_{p}\>\>\>}\Bigg{]}$
(4)

where $\mathrm{BR}_{p}$ is the branching ratio of pattern $p$,
$\big{[}N^{(sel)}_{MC}\big{]}^{(s)}_{p}$ and
$\big{[}N^{(tot)}_{MC}\big{]}_{p}$ are respectively the selected and total
number of simulated decays in the de-excitation pattern of interest. For
$0\nu\beta\beta$ decay signatures the signal is monochromatic in all the
involved crystals, so the signal region is expected to lie around a specific
point in the M-dimensional space of coincident energy releases. A selection is
enforced, in simulations, with a box cut, i.e. a selection interval for energy
releases in each crystal, defined as

$|E_{i}-Q_{i}|<5\>\mathrm{keV}\quad\mathrm{where}\quad i=1...M$ (5)

where $E_{i}$ is the reconstructed energy release in the ordered energy
space333The energy releases of each M-bolometers multiplet are ordered in
descending order so that $E_{i}>E_{i+1}$. and $Q_{i}$ is the corresponding
expected energy release. For $2\nu\beta\beta$ decay signatures the same
selections apply except the one crystal where the energy release from the
$\beta\beta$ is expected. Since the emitted neutrinos carry away an unknown
(on an event basis) amount of undetected energy, the expected energy release
is not monochromatic. It is instead expected to vary from $Q_{j}^{min}$ to
$Q_{j}^{max}$ where $j$ indicates the bolometer where the $\beta\beta$ release
their energy. For that bolometer, in each multiplet the following selection is
applied

$Q_{j}^{min}-5\>\mathrm{keV}<E_{j}<Q_{j}^{max}+5\>\mathrm{keV}$ (6)

Selection cuts need to be further tuned at a later stage to optimize the
sensitivity to signal peaks. We do this including the widest possible
sidebands around each signal peak in order to best constrain the underlying
continuous background. We try to avoid including background peaks in the fit
range, in order to minimize systematics due to their modeling. This process
yields the selections listed in Table 3 (see Sec. 3.4). We then update our
computation of signal efficiencies using Eq. 4, where $N^{sel}_{MC}$ is
replaced by the result of a Gaussian fit to the distribution of selected MC
signal events.

The second efficiency contribution, namely the analysis efficiency, is the
combination of the probability of correctly detecting and reconstructing the
energy deposited in each bolometer (cut efficiency, $\varepsilon^{cut}$), and
the probability of assigning the correct multiplicity and avoiding an
accidental coincidence (accidentals efficiency, $\varepsilon^{acc}$).

The cut efficiency term is named after the data processing cuts needed to
select triggered events that pass the base and PSA cuts (see Sec. 2). The
method used to calculate this efficiency follows closely what was used in
Adams:2019jhp . We measure the efficiency of correctly triggering,
reconstructing the pulse energy and the pile-up contribution (base cuts) from
heater pulses. The base cut efficiency is computed on each bolometer-dataset
pair given the large number of available heater events, and then averaged to
obtain a per-dataset value. The PSA cut efficiency is extracted from two
independent samples of events: either coincident double-site events where the
total energy released is compatible with known prominent $\gamma$ lines, or
single-site events due to fully absorbed $\gamma$ lines. The first sample
includes events whose energy spans a wide range, and allows the determination
of the PSA cut efficiency dependence on energy. The second sample has a higher
statistics but provides a measurement at the energies of the selected $\gamma$
peaks only, rather than on a continuum. We evaluate for each dataset the PSA
cut efficiency term as the average of the two efficiencies obtained from such
samples. We treat the difference as a systematic effect. The
$\varepsilon^{cut}$ term must be raised to the $M^{th}$ power because it
models bolometer-related efficiencies and a multiplet is selected if and only
if all of the involved bolometers pass the selection cuts.

The accidentals efficiency term $\varepsilon^{acc}$ is obtained, separately
for each dataset, as the survival probability of the 40K $\gamma$ line at 1460
keV. A fully absorbed 40K line in CUORE is uncorrelated to any other physical
event because it follows an electron capture of a $\sim 3$ keV shell, which is
below threshold.

Summarizing, the total signal detection efficiency for signature $s$ is

$\varepsilon^{tot}=\varepsilon_{s}\times(\varepsilon^{cut})^{M}\times\varepsilon^{acc}.$
(7)

Since the cut efficiency and accidentals term are evaluated separately for
each dataset, the total efficiency term in Eq. 7 must be thought of as the
signal efficiency for signature $s$ for a specific dataset. A summary of the
relevant efficiency values is provided in Table 2, where the per-dataset
values are exposure weighted over all datasets. The containment efficiency is
the dominant term.

Table 2:  We report the efficiency terms that appear in Eq. 7 separately for
the $0\nu\beta\beta$ (top) and $2\nu\beta\beta$ (bottom) analyses. The
containment term dominates the efficiency. We report the cut efficiency raised
to power $M$ according to the signature it refers to. We quote effective
values computed as exposure weighted mean for the cut and accidentals
efficiency terms. All values are percentages, the uncertainty on the last
digit is included in round brackets.

$0\nu\beta\beta$
---
| 2A0-2B1 | 2A1-2B2 | 3A0
Containment | $4.6(2)$ | $2.9(1)$ | $2.5(1)$
Cut | $78.7(2)$ | $69.8(3)$
Accidentals | $98.7(1)$
Total | $3.5(2)$ | $2.3(1)$ | $1.7(1)$

$2\nu\beta\beta$
---
| 2A0-2B1 | 2A1-2B2 | 3A0
Containment | $4.2(2)$ | $2.4(1)$ | $0.19(1)$
Cut | $78.7(2)$ | $69.8(3)$
Accidentals | $98.7(1)$
Total | $3.2(1)$ | $1.9(1)$ | $0.13(1)$

###  3.3 Background Contributions

Radioactive decays and particle interactions other than $\mathrm{{}^{130}Te}$
decay to $\mathrm{{}^{130}Xe}$ excited state, may mimic the process we search
for. We estimate this background contribution by means of background MC
simulations described in Sec 3.1. We combine background simulations of
different sources, according to the CUORE background model, and from the
simulated background spectra we compute the expected number of background
counts for each signature $B_{s}$, by counting the expected events from each
source included in the background model, and summing the contributions from
all sources. We apply the same tight selection cuts around the signal region
defined in Eqs. 5 and 6.

We use $B_{s}$ to evaluate an approximate sensitivity for each signature and
ultimately select the ones that will enter the analysis (see Sec. 3.4). Once
the signatures that enter the analysis are selected, we optimize the selection
cuts around the signal region in order to reject background structures while
leaving the widest possible sidebands around the expected signal position. In
this way we can parameterize the background with an appropriate analytical
function, whose shape is dictated by background simulations, and use that to
perform the final analysis (see Sec. 4.1). With this method we infer the
number of reconstructed background events in each signature from data, rather
than relying just on simulations.

###  3.4 Experimental Signature Ranking

The 15 unique signatures under analysis have different signal efficiencies and
backgrounds, and thus different detection sensitivities of the signal. In this
section we evaluate an approximate sensitivity of each signature and reduce
the 15 signatures down to the most sensitive subset.

We analytically evaluate the discovery sensitivity of signature $s$ starting
from a background-only model for the total number of counts observed in a
single bin centered at the expected signal position. In background-free
signatures $B_{s}\ll 1$ we assume an exponentially decaying prior
P$(\mu)=e^{-\mu}$ where $\mu$ is the true value of the number of background
counts. In background-limited ones $B_{s}\gg 1$ we assume a Gaussian prior
whose mean and variance are $B_{s}$. We define the discovery sensitivity as
the minimum number of observed counts $N_{s}$ such that the probability of
observing $N>N_{s}$ counts in the background-only model is smaller than a
given threshold $p_{th}$. Then, from $N_{s}$, we extract the corresponding
half life sensitivity

$\tilde{S}_{1/2}(\varepsilon_{s},B_{s})=\bigg{[}\frac{\ln(2)\>M\Delta
t\>N_{A}\>\eta(^{130}\mathrm{Te})}{m(\mathrm{TeO}_{2})}\bigg{]}S(\varepsilon_{s},B_{s})$
(8)

where $M$ is the detector mass, $\Delta t$ its live time, $N_{A}$ the Avogadro
constant, $\eta(^{130}\rm{Te})=(34.167\pm 0.002)\%$ Fehr:2004  the isotopic
abundance of ${}^{130}\mathrm{Te}$ in natural tellurium,
$m(\mathrm{TeO}_{2})=159.6$ g/mol the molecular mass of a tellurium dioxide
molecule Fehr:2004  and $S(\varepsilon_{s},B_{s})$ is a score function

$S(\varepsilon_{s},B_{s})=\begin{cases}\frac{\varepsilon_{s}}{-\ln(p_{th})}&B_{s}<B_{th}\\\
\frac{\varepsilon_{s}}{n_{\sigma}(p_{th})\sqrt{B}_{s}}&B_{s}\geq
B_{th}\end{cases}$ (9)

where $n_{\sigma}(p_{th})$ is the number of Gaussian sigma which correspond to
$p_{th}$, and $B_{th}$ sets the transition from the background-free
approximation to the background-limited approximation making
$S(\varepsilon_{s},B_{s})$ continuous. For $n_{\sigma}=5$, $p_{th}\sim 3\times
10^{-7}$ and $B_{th}\sim 9$. We note though that all signatures have a number
of expected background counts either $<1$ or $>10$ and their ranking would not
be affected by a different choice of the $B_{th}$ threshold.

We compute the relative score of each signature $s$ as

$R_{s}\doteq\frac{S(\varepsilon_{s},B_{s})}{\sum_{s^{\prime}}S(\varepsilon_{s^{\prime}},B_{s^{\prime}})}$
(10)

where $s^{\prime}$ is an index running on all experimental signatures. The
efficiency term $\epsilon_{s}$ only includes the MC-based containment
efficiency term. This is acceptable for the computation of an approximate
analytical score function, since the containment term by far dominates the
overall efficiency (Tab. 2).

We set a threshold of $R_{s}>5\%$ and we identify three signatures both for
the $0\nu\beta\beta$ and $2\nu\beta\beta$ decay search, namely the 2A0-2B1,
2A2-2B3 and 3A0 listed in Table 3. The selected experimental signatures
account for a majority of the sensitivity contributions among studied
signatures. The sum of their scores accounts for $84\%\>(87\%)$ total score of
all signatures in the $0\nu\beta\beta\>(2\nu\beta\beta)$ search respectively.

Table 3: Selected experimental signatures for DBD search on the $0^{+}_{2}$
excited state of ${}^{130}\mathrm{Xe}$ in the $0\nu\beta\beta$ (top) and
$2\nu\beta\beta$ (bottom) channel are listed. For each signature the
corresponding Regions Of Interest (ROI, i.e. the applied selection cuts) are
listed in terms of the ordered energy releases $E_{1}\geq E_{2}\geq E_{3}$.
The component that will be used for the fit is highlighted with a ∗
superscript. For each signature the partition of the secondaries expected to
contribute are listed. For each secondary we report in round brackets the
expected energy release in keV. The last row reports the corresponding
relative score (Eq. 10).

$0\nu\beta\beta$
---
2A0-2B1 |  2A2-2B3 |  3A0
$E_{2}:\gamma(1257)$ |  $E_{2}:\gamma(536)$ |  $E_{3}:\gamma(536)$
$E_{1}:\beta\beta(734)+\gamma(536)$ |  $E_{1}:\beta\beta(734)+\gamma(1257)$ |  $E_{2}:\beta\beta(734)$
|  |  $E_{1}:\gamma(1257)$
$1247<E_{2}<1280$ |  $523<E^{*}_{2}<573$ |  $526<E_{3}<546$
$1247<E_{1}^{*}<1280$ |  $1981<E_{1}<2001$ |  $700<E^{*}_{2}<760$
|  |  $1247<E_{1}<1267$
$39\%$ |  $25\%$ |  $20\%$

$2\nu\beta\beta$
---
2A0-2B1 |  2A2-2B3 |  3A0
$E_{2}:\beta\beta(0-734)+\gamma(536)$ |  $E_{2}:\gamma(536)$ |  $E_{3}:\beta\beta(0-734)$
$E_{1}:\gamma(1257)$ |  $E_{1}:\beta\beta(0-734)+\gamma(1257)$ |  $E_{2}:\gamma(536)$
|  |  $E_{1}:\gamma(1257)$
$620<E_{2}<1150$ |  $523<E^{*}_{2}<573$ |  $400<E_{3}<523$
$1220<E^{*}_{1}<1300$ |  $1360<E_{1}<1990$ |  $523<E^{*}_{2}<573$
|  |  $1779<E_{1}+E_{2}<1807$
$40\%$ |  $22\%$ |  $25\%$

##  4 Physics Extraction

We use a phenomenological parameterization of the background in the fitting
regions (as opposed to using the predicted spectra from the MC), hence real
data are required to tune the fit. To avoid biasing our results, we build the
fit (Sec. 4.1) on blinded data using the blinding procedure described in Sec.
4.2.

###  4.1 Fitting Technique

We extract the $0\nu\beta\beta$ decay rate and $2\nu\beta\beta$ decay rate of
${}^{130}\mathrm{Te}$ to the $0^{+}_{2}$ excited state of
${}^{130}\mathrm{Xe}$ using two separate Bayesian fits. For the single process
($0\nu$ or $2\nu$) the fit is run simultaneously on all the involved
signatures. Every multiplet of multiplicity $M$ can be represented as a point
$\vec{E}_{ev}$ in a $M$-dimensional space of reconstructed energies. The
energy releases are ordered so that $E_{i}>E_{i+1}\>\>\forall i=1,..,M-1$. For
each signature, one of the components of the $\vec{E}_{ev}$ vector is selected
to perform the fit. This is referred to as projected energy, and indicated
with a ∗ superscript in Table 3. In the following we will denote this energy
as $E_{ev}$.

An unbinned Bayesian fit is implemented with the BAT software package
Caldwell:2008fw . It allows simultaneous sampling and maximization of the
posterior probability density function (pdf) via Markov Chain Monte Carlo. The
likelihood can be decomposed, for each signature and dataset, as follows:

$\begin{split}\log\mathcal{L}_{s,ds}&=-(\lambda_{S\>s,ds}+\lambda_{B\>s,ds})+\\\
&+\sum_{ev\in(s,ds)}\log\bigg{[}\lambda_{S\>s,ds}\xi_{bo,ds}f_{S}(E_{ev})+\\\
&+\lambda_{B\>s,ds}\xi_{bo,ds}f_{B}(E_{ev})\bigg{]}\end{split}$ (11)

where the subscripts $s$, $bo$, $ds$ will be used to refer to a specific
signature, bolometer, or dataset respectively. The form of Eq. 11 is the same
for $0\nu\beta\beta$ and $2\nu\beta\beta$, the $\lambda_{S}$ and $\lambda_{B}$
terms are the expected number of signal and background events respectively,
$\xi_{bo,ds}$ is the ratio between the exposure of bolometer $bo$ to the
exposure of dataset $ds$, $f_{S}$ and $f_{B}$ are the normalized signal and
background pdfs. They depend just on the projected energy variable $E_{ev}$.

The response function of CUORE bolometers to monochromatic energy releases has
a functional form defined phenomenologically for each bolometer-dataset pair
Alexey2018  Laura2018  as the superposition of 3 Gaussian components to
account for non-Gaussian tails. A correction for the bias in the energy scale
reconstruction is implemented together with the resolution dependence on
energy (see Ref. Alduino:2017ehq  for more details). The signal term
$f_{S}(E_{ev})$ models such shape in the bolometer-dataset pair the projected
energy $E_{ev}$ was released in. The expected number of signal counts can be
written as

$\begin{split}\lambda_{S\>s,ds}&=\Gamma^{(p)}_{\beta\beta}[\mathrm{yr}^{-1}]\bigg{[}\frac{N_{A}\>10^{3}\>\eta(^{130}\mathrm{Te})}{m(\mathrm{TeO}_{2})\>[\rm{g/mol}]}\bigg{]}\epsilon_{s}\cdot\\\
&\cdot(M\Delta t)_{ds}\>\mathrm{[kg\cdot
yr]}\>(\epsilon^{cut})_{ds}^{M}\>(\epsilon^{acc})_{ds}\end{split}$ (12)

where $\Gamma^{(p)}_{\beta\beta}$ is the decay rate of process $p$ and the
other parameters were introduced following Eq. 8. The
$\Gamma_{\beta\beta}^{(p)}$ parameter describes the rate of the process
$p=0\nu,2\nu$ under investigation and is given in both cases a uniform
physical prior, $\Gamma_{\beta\beta}>0$.

The background term $f_{B}(E_{ev})$ is parameterized as

$f_{B}(E_{ev})=\frac{1}{\Delta
E}\bigg{[}1+m_{s}\big{(}E_{ev}-E^{(s)}_{0}\big{)}\bigg{]}$ (13)

where $\Delta E=E^{max}_{s}-E^{min}_{s}$ is the width of the region of
interest, $E^{(s)}_{0}$ is the center of it, and $m_{s}$ describes the slope
of the background for signature $s$. The normalization of the background term
represents the number of expected background counts

$\lambda_{B\>s,ds}=\mathrm{BI}_{s}\>(M\Delta
t)_{ds}\>(E^{max}_{s}-E^{min}_{s})$ (14)

where BIs is the background index for signature $s$. Background simulations
suggest that a uniform event distribution is enough to describe the continuous
background in all signatures except the $2\nu\beta\beta$ 2A0-2B1. For this
reason the $m_{s}$ parameter is included only when necessary. The background
is fully described by the BIs and $m_{s}$ which, together, make 4 (3) nuisance
parameters in the $2\nu\beta\beta$ $(0\nu\beta\beta)$ case respectively, that
will be marginalized over. The prior for background indices BIs and slopes
$m_{s}$ is uniform.

The combined log-likelihood reads

$\log\mathcal{L}(\mathcal{D}|\mathrm{H}_{S+B})=\sum_{s,ds}\log\mathcal{L}_{s,ds}$
(15)

where $\mathrm{H}_{S+B}$ indicates that the likelihood is written in the
signal-plus-background model hypothesis HS+B, i.e. that the existence of the
process of interest is assumed. The background-only hypothesis HB is a
particular case that can be obtained by setting $\Gamma_{\beta\beta}=0$.

###  4.2 Blinding and Sensitivity

We blind the data by injecting simulated signal events into the experimental
spectrum. We inject a random and unknown number of fake signal events that
would correspond to an event rate larger than the current 90 % upper limit
Adams:2018yrj ; Andreotti:2011in . Then we compute the expected number of
counts in each signature, according to known efficiencies and exposures, for
each dataset. Each generated signal event is randomly assigned a bolometer,
according to its exposure within the considered dataset. Finally the projected
energy of the signal event is generated according to the detector response
function $f_{S}(E_{ev}|Q_{s})$ centered at the expected position $Q_{s}$ of
the monochromatic energy release in the projected energy space.

The injection rate of the simulated signal events is comprised between:

$\Gamma^{min}_{p}=1\cdot
10^{-23}\>[\mathrm{1/yr}]\quad\mathrm{and}\quad\Gamma^{max}_{p}=5\cdot
10^{-23}\>[\mathrm{1/yr}]$ (16)

We then fit the blinded datasets to get data driven estimates of the
background levels in each fitting window. These background estimates are used
as inputs to our sensitivity studies in the next section.

The results of the fits to the blinded data are reported in Table 4. We see a
non-null background for both the $\mathrm{2A0-2B1}$ and $\mathrm{2A2-2B3}$
signatures. No background is expected for the $\mathrm{3A0}$ signature.

To extract the median half-life sensitivity for each decay, we generate
$10^{4}$ background-only Toy Monte Carlo simulations (ToyMC), using the
numbers in Tab. 4. A background-only ToyMC simulation is an ensemble of
simulated datasets, according to the following procedure which is iterated
$N_{toy}$ times, to produce the same number of ToyMC ensembles. We define a
set of signatures, together with the multiplicity and cuts in the ordered
energy variables that identify candidate events. For each signature, we set a
functional form for the background pdf, either constant or linear, and sample
a value from the posterior pdf of the corresponding blinded fit. We compute
the number of expected background events for each signature and dataset
according to Eq. 14 and sample the actual number of background events from a
Poisson distribution with expectation value equal to the number of expected
counts.

We store each simulated ToyMC event as a vector of ordered energy releases
$\vec{E}_{ev}$ and related bolometers $\vec{\mathrm{ch}}_{ev}$, where the
bolometers are randomly extracted from the active bolometers of each dataset
according to their exposure in the data, while the energies are generated
according to the selected shape of the background pdf computed with the
parameters (e.g. background index) generated according to the posterior pdfs
obtained with the blinded fit to the data.

We then fit each ToyMC with the signal-plus-background model
$\mathrm{H}_{S+B}$ and compute the lower limit for the decay half life from
the $90\%$ quantile of the marginalized posterior pdf for the decay rate
parameter. We show the distribution of such limits in Figure 2 for both the
$0\nu\beta\beta$ and $2\nu\beta\beta$ decay process. We quote the half-life
sensitivity as the median limit of the ToyMCs (Table 4). They are
respectively: $\mathrm{S^{0\nu}_{1/2}=(5.6\pm 1.4)\times
10^{24}\>\mathrm{yr}}$ for the ${0\nu}$ decay and
$\mathrm{S^{2\nu}_{1/2}=(2.1\pm 0.5)\times 10^{24}\>\mathrm{yr}}$ for the
${2\nu}$ decay where the uncertainty is the MAD of the corresponding
distribution.

Figure 2:  Distribution of $90\%$ C.I. marginalized upper limits on
$\mathrm{T_{1/2}}=\log 2/\Gamma$ for $0\nu\beta\beta$ decay (top) and
$2\nu\beta\beta$ decay (bottom) obtained from Toy MC simulations. We obtain a
median sensitivity of $S^{0\nu}_{1/2}=2.1\times 10^{24}$ yr and
$S^{2\nu}_{1/2}=5.6\times 10^{24}$ yr (black dashed line), compared to the
$90\%$ C.I. limit from this analysis (red solid line).  Table 4: Results of
the blinded fit to $0\nu\beta\beta$ (top) and $2\nu\beta\beta$ (bottom)
candidate events in the signatures of Table 3. For each parameter the mean and
standard deviation of the corresponding marginalized posterior distribution
are reported. These values are used only as input to the sensitivity studies
and fit validation. Final results reported in Table 5.

$0\nu\beta\beta$
---
Observable | Blinded Fit Value | Units
Blinded $\mathrm{\Gamma_{\beta\beta}^{0\nu}}$ | $2.8\pm{0.1}$ | $\mathrm{10^{-23}\>[yr^{-1}]}$
$\mathrm{BI_{2A0-2B1}}$ | $6.1\pm{3.6}$ | $\mathrm{10^{-4}\>[counts/(keV\>kg\>yr)]}$
$\mathrm{BI_{2A2-2B3}}$ | $2.7\pm{1.5}$ | $\mathrm{10^{-4}\>[counts/(keV\>kg\>yr)]}$
$\mathrm{BI_{3A0}}$ | $8.7\pm{8.3}$ | $\mathrm{10^{-5}\>[counts/(keV\>kg\>yr)]}$

$2\nu\beta\beta$
---
Observable | Blinded Fit Value | Units
Blinded $\mathrm{\Gamma_{\beta\beta}^{2\nu}}$ | $\>5.1\pm{0.1}$ | $\mathrm{10^{-23}\>[yr^{-1}]}$
$\mathrm{BI_{2A0-2B1}}$ | $\>3.3\pm{0.4}$ | $\mathrm{10^{-3}\>[counts/(keV\>kg\>yr)]}$
$\mathrm{m_{2A0-2B1}}$ | $-5.5\pm{4.3}$ | $\mathrm{10^{-3}\>[1/keV]}$
$\mathrm{BI_{2A2-2B3}}$ | $\>4.0\pm{0.6}$ | $\mathrm{10^{-3}\>[counts/(keV\>kg\>yr)]}$
$\mathrm{BI_{3A0}}$ | $\>6.9\pm{6.8}$ | $\mathrm{10^{-5}\>[counts/(keV\>kg\>yr)]}$

##  5 Results

Table 5:  We report here mean and standard deviation of the marginalized
posterior distributions for the decay rate and background parameters for each
signature, derived from unblinded data. The S1/2 parameter indicates the
median expected sensitivity for limit setting at $90\%$ C.I. on the T1/2
parameter together with the MAD of its distribution. The last row reports the
marginalized $90\%$ C.I. Bayesian lower limit on the decay half life. All
results come from the combined fit with systematics.

$0\nu\beta\beta$
---
Parameter | Final Fit Value | Units
$\mathrm{\Gamma_{\beta\beta}^{0\nu}}$ | $5.8\pm 4.5$ |  $10^{-26}$ [yr-1]
BI2A0-2B1 | $2.1\pm 1.4$ |  $10^{-4}$ [counts/(keV kg yr)]
BI2A2-2B3 | $2.7\pm 1.2$ |  $10^{-4}$ [counts/(keV kg yr)]
BI3A0 | $5.8\pm 5.5$ |  $10^{-5}$ [counts/(keV kg yr)]
$\mathrm{S^{0\nu}_{1/2}}$ | $5.6\pm 1.4$ |  $10^{24}$ [yr]
$\mathrm{T}_{1/2}^{90\%}$ | $>5.9$ |  $10^{24}$ [yr]

$2\nu\beta\beta$
---
Parameter | Final Fit Value | Units
$\mathrm{\Gamma_{\beta\beta}^{2\nu}}$ | $2.8\pm 1.8$ |  $10^{-25}$ [yr-1]
BI2A0-2B1 | $3.0\pm 0.3$ |  $10^{-3}$ [counts/(keV kg yr)]
$m_{\mathrm{2A0-2B1}}$ | $-5.2\pm 4.2$ |  $10^{-3}$ [keV-1]
BI2A2-2B3 | $4.3\pm 0.5$ |  $10^{-3}$ [counts/(keV kg yr)]
BI3A0 | $5.4\pm 5.4$ |  $10^{-5}$ [counts/(keV kg yr)]
$\mathrm{S^{2\nu}_{1/2}}$ | $2.1\pm 0.5$ |  $10^{24}$ [yr]
$\mathrm{T}_{1/2}^{90\%}$ | $>1.3$ |  $10^{24}$ [yr]

Figure 3: Result of the unbinned fit plotted on binned data. Error bars are
just a visual aid and correspond to the square root of the bin contents. We
show the best fit curve (blue solid), its signal component (blue dashed) at
the global mode of the posterior for $0\nu\beta\beta$ (left) and
$2\nu\beta\beta$ (right), and the 90 % C.I. marginalized limit on the decay
rate (red solid).

We show in Figure 3 and Table 5 the results of the fit to unblinded data for
both $0\nu\beta\beta$ and $2\nu\beta\beta$. Though the data are binned for
graphical reasons, the analysis is unbinned. Including contributions from all
sources of systematic uncertainty listed in Table 6, no significant signal is
observed in either decay mode. The global mode of the joint posterior pdf for
the rate parameter is

$\displaystyle\hat{\rm{\Gamma}}^{0\nu}_{\beta\beta}=4.0_{-4.0}^{+3.0}\times
10^{-26}\>\rm{yr}^{-1}$ (17a)
$\displaystyle\hat{\rm{\Gamma}}^{2\nu}_{\beta\beta}=2.2_{-1.9}^{+1.7}\times
10^{-25}\>\rm{yr}^{-1}$ (17b)

whereas we quote the uncertainty as the marginalized $68~{}\%$ smallest
interval. We report the marginalized posterior pdf and the $90\%$ C.I. in
Figure 4. We include the marginalized posterior pdfs for individual background
parameters in Figure 5 $(0\nu\beta\beta)$ and Figure 6 $(2\nu\beta\beta)$.
Their means agree with the corresponding results from blinded data within one
standard deviation. We observe a slight negative background fluctuation (i.e.
limit stronger than expected) in the $0\nu\beta\beta$ decay analysis and a
positive one (i.e. limit looser than expected) in the $2\nu\beta\beta$ decay
analysis with respect to the median $90\%$ C.I. limit. The probability of
setting an even stronger (looser) limit in the $0\nu\beta\beta$
($2\nu\beta\beta$) decay analysis is respectively $45.1\%$ and $10.4\%$. Null
signal rates are included in the $1\sigma$ ($2\sigma$) smallest C.I.444We
refer here to the Gaussian case to define the probability content of any
$n\sigma$ interval. of the marginalized pdf for the
$\rm{\Gamma}^{0\nu}_{\beta\beta}$ ($\rm{\Gamma}^{2\nu}_{\beta\beta}$) rate
parameter respectively. The following Bayesian lower bounds on the
corresponding half life parameters are set:

$\displaystyle\big{(}\mathrm{T_{1/2}}\big{)}^{0\nu}_{0^{+}_{2}}$
$\displaystyle>5.9\times 10^{24}\>\mathrm{yr}\quad(\mathrm{90\%\>\>C.I.})$
(18a) $\displaystyle\big{(}\mathrm{T_{1/2}}\big{)}^{2\nu}_{0^{+}_{2}}$
$\displaystyle>1.3\times 10^{24}\>\mathrm{yr}\quad(\mathrm{90\%\>\>C.I.})$
(18b)

The results reported in this Article represent the most stringent limits on
the DBD of 130Te to excited states and improve by a factor $\sim 5$ the
previous results on this process.

###  5.1 Systematic Uncertainties

The major sources of systematic uncertainty are: signal efficiencies, the
detector response function, energy calibration, and the uncertainty of the
isotopic abundance of 130Te (see Table 6). Each systematic uncertainty can be
introduced as a set of nuisance parameters in the fit with a specified prior
distribution. Each set of nuisance parameters can be activated independently
with the priors listed in Table 6. We individually monitor the effect of
activating each source of systematic uncertainty repeating the fit and
comparing the 90% C.I. Bayesian limit on the half life with respect to the
minimal model, where we describe all sources of systematics with constants
rather than fit parameters. Finally, we repeat the fit activating all
additional nuisance parameters at once.

We include, for each dataset, two separate parameters to model different
sources of uncertainty in the cut efficiency evaluation, and replace the
$\varepsilon^{cut}$ constant with the sum of
$\varepsilon^{cut\>(I)}+\varepsilon^{cut\>(II)}$. We refer to cut efficiency I
to parameterize the uncertainty due to the finiteness of the samples of pulser
events and $\gamma$ decays used to extract the cut efficiency. Its prior is
Gaussian with mean equal to $\varepsilon^{cut}$ and width equal to the
corresponding uncertainty (see Table. 2). The additional cut efficiency II is
uniformly distributed, with zero mean. It models the systematic uncertainty
due to the PSA efficiency, shifting the cut efficiency by at most $0.7\%$. The
accidentals efficiency contributes to the systematic uncertainty just with the
uncertainty due to limited statistics in the 40K peak. We add one nuisance
parameter per dataset with a Gaussian distributed prior to model this effect.
The containment efficiency is instead affected by uncertainty due to the
simulation of Compton scattering events at low energy. The uncertainty due to
the number of simulated signal events is negligible. We take the ratio of the
Compton scattering attenuation coefficient to reference data Allison:2016lfl
as a measure of the relative uncertainty on this efficiency term. We account
for this, for each signature, introducing a nuisance containment efficiency
parameter with a Gaussian prior. The 130Te natural isotopic abundance on
natural Te is modeled with a single global nuisance parameter with a Gaussian
prior $\eta=(34.167\pm 0.002)\%$. Both the detector response function shape
and the energy scale bias are evaluated from data, as anticipated in Sec. 4.1.
Each effect is separately parameterized with a 2nd order polynomial as a
function of energy, whose coefficients are evaluated with a fit to the 5-7
most visible peaks in each dataset. The uncertainty and correlations among
such parameters are themselves a source of systematic uncertainty, and are
included as 2 independent sets of correlated parameters per dataset, with a
multivariate normal prior distribution. In this way, the detector response
function width and position are allowed to float within their uncertainty.

Uncertainty in modeling the detector response function leads to the dominant
systematic effect on the limit, which is below a 1% shift. Sub-dominant
effects come from the energy scale bias and the containment efficiency (Table
6).

Source | Prior | Effect on $\mathrm{T}_{1/2}^{90\%}$
---|---|---
$0\nu\beta\beta$ | $2\nu\beta\beta$
Cut efficiency I | Gaussian | $0.2\%$ | $<0.1\%$
Cut efficiency II | Uniform | $0.1\%$ | $0.1\%$
Accidentals efficiency | Gaussian | $0.2\%$ | $<0.1\%$
Containment efficiency | Gaussian | $0.3\%$ | $0.1\%$
${}^{130}\mathrm{Te}$ isotopic abundance | Gaussian | $<0.1\%$ | $0.1\%$
Energy scale bias | Multiv. Gaussian | $0.3\%$ | $0.1\%$
Detector resolution | Multiv. Gaussian | $0.5\%$ | $0.4\%$
Combined | Multivariate | $0.4\%$ | $0.1\%$
Table 6:  Systematic uncertainties. We report each effect separately and their
combination on the marginalized $90\%$ C.I. T${}_{1/2}^{90\%}$ limit on the
$0\nu\beta\beta$ and $2\nu\beta\beta$ decay half life.  Figure 4: Marginalized
decay rate posterior pdf for $0\nu\beta\beta$ (top) and $2\nu\beta\beta$
(bottom) from the combined fit with all systematics. We show the $90\%$ C.I.
in gray.

##  6 Conclusions

We have presented the latest search for DBD of 130Te to the first $0^{+}$
excited state of 130Xe with CUORE based on a 372.5 kg$\cdot$yr TeO2 exposure.
We found no evidence for either $0\nu\beta\beta$ nor $2\nu\beta\beta$ decay
and we placed a Bayesian lower bound at $90\%$ C.I. on the decay half lives of
$\mathrm{(T_{1/2})^{0\nu}_{0^{+}_{2}}>5.9\times 10^{24}\>\mathrm{yr}}$ for the
$0\nu\beta\beta$ mode and $\mathrm{(T_{1/2})^{2\nu}_{0^{+}_{2}}>1.3\times
10^{24}\>\mathrm{yr}}$ for the $2\nu\beta\beta$ mode.

The median limit setting sensitivity for the $2\nu\beta\beta$ decay of
$2.1\times 10^{24}$yr is starting to approach the $7.2\times 10^{24}$yr lower
bound of the QRPA theoretical predictions half life range for this decay mode.
The CUORE experiment is steadily taking data at an average rate of $\sim 50$
kg$\cdot$yr/month and by the end of its data taking the collected exposure is
expected to increase by an order of magnitude. Work is ongoing to improve the
sensitivity by extending the analysis to not-fully-contained events,
leveraging the information from the topology of higher dimension coincident
signal multiplets to further reduce the background, and improving the signal
efficiency by lowering the threshold of the pulse shape discrimination
algorithm.

###### Acknowledgements.

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 Fisica Nucleare (INFN); the
National Science Foundation under Grant Nos. NSF-PHY-0605119, NSF-PHY-0500337,
NSF-PHY-0855314, NSF-PHY-0902171, NSF-PHY-0969852, NSF-PHY-1307204, NSF-
PHY-1314881, NSF-PHY-1401832, and NSF-PHY-1913374; and Yale University. This
material is also based upon work supported by the US Department of Energy
(DOE) Office of Science under Contract Nos. DE-AC02-05CH11231 and DE-
AC52-07NA27344; by the DOE Office of Science, Office of Nuclear Physics under
Contract Nos. DE-FG02-08ER41551, DE-FG03-00ER41138, DE-SC0012654, DE-
SC0020423, DE-SC0019316; and by the EU Horizon2020 research and innovation
program under the Marie Sklodowska-Curie Grant Agreement No. 754496. This
research used resources of the National Energy Research Scientific Computing
Center (NERSC). 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. We acknowledge contributions from
J. Kotila on dedicated computation of beta spectra.

## References

  * [1] G. Racah.  Sulla simmetria tra particelle e antiparticelle.  Nuovo Cim., 14:322–328, 1937.
  * [2] W.H. Furry.  On transition probabilities in double beta-disintegration.  Phys. Rev., 56:1184–1193, 1939.
  * [3] B. Pontecorvo.  Neutrino experiments and the problem of conservation of leptonic charge.  Sov. Phys. JETP, 26:984, 1968.
  * [4] J. Schechter and J.W.F. Valle.  Neutrinoless Double Beta Decay in SU(2) x U(1) Theories.  Phys. Rev. D, 25:2951, 1982.
  * [5] M. Fukugita and T. Yanagida.  Baryogenesis Without Grand Unification.  Phys. Lett. B, 174:45–47, 1986.
  * [6] F. T. Avignone, S. R. Elliott, and J. Engel.  Double Beta Decay, Majorana Neutrinos, and Neutrino Mass.  Rev. Mod. Phys., 80:481–516, 2008.
  * [7] A. S. Barabash et al.  Two neutrino double beta decay of 100Mo to the first excited $0^{+}$ state in 100Ru.  Phys. Lett., B345:408–413, 1995.
  * [8] A. S. Barabash, F. Hubert, P. Hubert, and V. I. Umatov.  Double beta decay of 150Nd to the first $0^{+}$ excited state of 150Sm.  JETP Lett., 79:10–12, 2004.  [Pisma Zh. Eksp. Teor. Fiz.79,12(2004)].
  * [9] A. S. Barabash.  Double beta decay to the excited states: Review.  AIP Conf. Proc., 1894(1):020002, 2017.
  * [10] E. Andreotti et al.  Search for double-$\beta$ decay of 130Te to the first 0+ excited state of 130Xe with CUORICINO.  Phys. Rev., C85:045503, 2012.
  * [11] C. Alduino et al.  Double-beta decay of ${}^{130}\hbox{Te}$ to the first $0^{+}$ excited state of ${}^{130}\hbox{Xe}$ with CUORE-0.  Eur. Phys. J., C79(9):795, 2019.
  * [12] Pekka Pirinen and Jouni Suhonen.  Systematic approach to $\beta$ and $2\nu\beta\beta$ decays of mass $a=100-136$ nuclei.  Phys. Rev. C, 91:054309, 2015.
  * [13] Jouni Suhonen and Osvaldo Civitarese.  Probing the quenching of $g_{A}$ by single and double beta decays.  Phys. Lett. B, 725(1):153 – 157, 2013.
  * [14] J. Suhonen and O. Civitarese.  Single and double beta decays in the a=100, a=116 and a=128 triplets of isobars.  Nucl. Phys. A, 924:1 – 23, 2014.
  * [15] J. Barea, J. Kotila, and F. Iachello.  Nuclear matrix elements for double-$\beta$ decay.  Phys. Rev. C, 87(1):014315, 2013.
  * [16] J. Kotila and F. Iachello.  Phase space factors for double-$\beta$ decay.  Phys. Rev. C, 85:034316, 2012.
  * [17] A.S. Barabash.  Precise half-life values for two neutrino double beta decay.  Phys. Rev. C, 81:035501, 2010.
  * [18] Bjoern Lehnert.  Excited state transitions in double beta decay: A brief review.  EPJ Web of Conferences, 93:01025, 01 2015.
  * [19] C. Arnaboldi et al.  CUORE: A Cryogenic underground observatory for rare events.  Nucl. Instrum. Meth. A, 518:775–798, 2004.
  * [20] C. Brofferio, O. Cremonesi, and S. Dell’Oro.  Neutrinoless Double Beta Decay Experiments With TeO2 Low-Temperature Detectors.  Front. in Phys., 7:86, 2019.
  * [21] C. Alduino et al.  First Results from CUORE: A Search for Lepton Number Violation via $0\nu\beta\beta$ Decay of 130Te.  Phys. Rev. Lett., 120(13):132501, 2018.
  * [22] D.Q. Adams et al.  Improved Limit on Neutrinoless Double-Beta Decay in 130Te with CUORE.  Phys. Rev. Lett., 124(12):122501, 2020.
  * [23] E. Fiorini and T.O. Niinikoski.  Low Temperature Calorimetry for Rare Decays.  Nucl. Instrum. Meth. A, 224:83, 1984.
  * [24] C. Enss and D. Mc Cammon.  Physical principles of low temperature detectors: Ultimate performance limits and cur- rent detector capabilities.  J. Low Temp. Phys., 151:5, 2008.
  * [25] C. Arnaboldi et al.  Production of high purity TeO2 single crystals for the study of neutrinoless double beta decay.  J. Cryst. Growth, 312(20):2999–3008, 2010.
  * [26] C. Alduino et al.  The CUORE cryostat: An infrastructure for rare event searches at millikelvin temperatures.  Cryogenics, 102:9–21, 2019.
  * [27] F. Alessandria et al.  The 4K outer cryostat for the CUORE experiment: Construction and quality control.  Nucl. Instrum. Meth. A, 727:65–72, 2013.
  * [28] E. Buccheri, M. Capodiferro, S. Morganti, F. Orio, A. Pelosi, and V. Pettinacci.  An assembly line for the construction of ultra-radio-pure detectors.  Nucl. Instrum. Meth. A, 768:130–140, 2014.
  * [29] C. Arnaboldi, P. Carniti, L. Cassina, C. Gotti, X. Liu, M. Maino, G. Pessina, C. Rosenfeld, and B.X. Zhu.  A front-end electronic system for large arrays of bolometers.  JINST, 13(02):P02026, 2018.
  * [30] S. Di Domizio, A. Branca, A. Caminata, L. Canonica, S. Copello, A. Giachero, E. Guardincerri, L. Marini, M. Pallavicini, and M. Vignati.  A data acquisition and control system for large mass bolometer arrays.  JINST, 13(12):P12003, 2018.
  * [31] G. Benato et al.  Radon mitigation during the installation of the CUORE 0$\nu\beta\beta$ decay detector.  JINST, 13(01):P01010, 2018.
  * [32] A. D’Addabbo, C. Bucci, L. Canonica, S. Di Domizio, P. Gorla, L. Marini, A. Nucciotti, I. Nutini, C. Rusconi, and B. Welliver.  An active noise cancellation technique for the CUORE Pulse Tube Cryocoolers.  Cryogenics, 93:56–65, 2018.
  * [33] E. E. Haller, N. P. Palaio, M. Rodder, W. L. Hansen, and E. Kreysa.  NTD germanium, A novel material for low temperature bolometers.  Neutron Transmutation Doping of Semiconductor Materials. edited by R. D. Larrabee (Springer, Boston), 1984.
  * [34] C. Alduino et al.  Analysis techniques for the evaluation of the neutrinoless double-$\beta$ decay lifetime in 130Te with the CUORE-0 detector.  Phys. Rev. C, 93(4):045503, 2016.
  * [35] E. Gatti and P.F. Manfredi.  Processing the Signals From Solid State Detectors in Elementary Particle Physics.  Riv. Nuovo Cim., 9N1:1–146, 1986.
  * [36] S. Di Domizio, F. Orio, and M. Vignati.  Lowering the energy threshold of large-mass bolometric detectors.  JINST, 6:P02007, 2011.
  * [37] A. Campani et al.  Lowering the Energy Threshold of the CUORE Experiment: Benefits in the Surface Alpha Events Reconstruction.  J. Low Temp. Phys., pages 1–10, 2020.
  * [38] G. Fantini et al.  Sensitivity to double beta decay of $\mathrm{{}^{130}Te}$ to the first $0^{+}$ excited state of $\mathrm{{}^{130}Xe}$ in CUORE.  J. of Phys.: Conference Series, 1468, 02 2020.
  * [39] Balraj Singh.  Nuclear data sheets for a= 130.  Nuclear Data Sheets, 93(1):33–242, 2001.
  * [40] S. Agostinelli et al.  Geant4 a simulation toolkit.  Nucl. Instrum. Meth. A, 506(3):250 – 303, 2003.
  * [41] C. Alduino et al.  The projected background for the CUORE experiment.  Eur. Phys. J. C, 77(8), Aug 2017.
  * [42] W.C. Haxton and G.J. Stephenson.  Double beta Decay.  Prog. Part. Nucl. Phys., 12:409–479, 1984.
  * [43] M. Doi, T. Kotani, H. Nishiura, K. Okuda, and E. Takasugi.  Neutrino Mass, the Right-handed Interaction and the Double Beta Decay. 1. Formalism.  Prog. Theor. Phys., 66:1739, 1981.  [Erratum: Prog.Theor.Phys. 68, 347 (1982)].
  * [44] M. Doi, T. Kotani, H. Nishiura, K. Okuda, and E. Takasugi.  Neutrino Mass, the Right-handed Interaction and the Double Beta Decay. 2. General Properties and Data Analysis.  Prog. Theor. Phys., 66:1765, 1981.  [Erratum: Prog.Theor.Phys. 68, 348 (1982)].
  * [45] D.Q. Adams et al.  Measurement of the $2\nu\beta\beta$ Decay Half-life of 130Te with CUORE.  in preparation, 2020.
  * [46] Manuela A Fehr, Mark Rehkämper, and Alex N Halliday.  Application of MC-ICPMS to the precise determination of tellurium isotope compositions in chondrites, iron meteorites and sulfides.  Int. J. of Mass Spectrometry, 232:83–94, 2004.
  * [47] Allen Caldwell, Daniel Kollar, and Kevin Kroninger.  BAT: The Bayesian Analysis Toolkit.  Comput. Phys. Commun., 180:2197–2209, 2009.
  * [48] A Drobizhev.  Searching for $0\nu\beta\beta$ decay of 130Te with the ton-scale CUORE bolometer array.  PhD thesis, University of California, Berkeley, 2018.
  * [49] L. Marini.  The CUORE experiment: from the commissioning to the first $0\nu\beta\beta$ limit.  PhD thesis, Università degli studi di Genova, 2018.
  * [50] J. Allison et al.  Recent developments in Geant4.  Nucl. Instrum. Meth. A, 835:186–225, 2016.

Figure 5: Marginalized posterior pdf for the background parameters of the
$0\nu\beta\beta$ model from the combined fit with all systematics. We show the
$68.3\%,95.5\%,99.7\%$ smallest C.I. in green, yellow and red respectively.
Figure 6: Marginalized posterior pdf for the background parameters of the
$2\nu\beta\beta$ model from the combined fit with all systematics. We show the
$68.3\%,95.5\%,99.7\%$ smallest C.I. in green, yellow and red respectively.