CUORE Sensitivity to $0\nu\beta\beta$ Decay

TL;DR

This study evaluates CUORE's potential to detect or exclude neutrinoless double beta decay, using Bayesian methods to estimate sensitivity based on background levels and energy resolution over different operational periods.

Contribution

It introduces a Bayesian analysis framework with toy Monte Carlo simulations to assess CUORE's sensitivity to $0 uetaeta$ decay, considering various experimental conditions.

Findings

CUORE can exclude $T_{1/2}^{0 u}$ up to 2×10^{25} years in 3 months.

Sensitivity improves to 9×10^{25} years over 5 years of data.

Discovery potential reaches 4×10^{25} years after 5 years.

Abstract

We report a study of the CUORE sensitivity to neutrinoless double beta ($0\nu\beta\beta$) decay. We used a Bayesian analysis based on a toy Monte Carlo (MC) approach to extract the exclusion sensitivity to the $0\nu\beta\beta$ decay half-life ($T_{1/2}^{0\nu}$) at $90\%$ credibility interval (CI) -- i.e. the interval containing the true value of $T_{1/2}^{0\nu}$ with $90\%$ probability -- and the $3 \sigma$ discovery sensitivity. We consider various background levels and energy resolutions, and describe the influence of the data division in subsets with different background levels. If the background level and the energy resolution meet the expectation, CUORE will reach a $90\%$ CI exclusion sensitivity of $2\cdot10^{25}$ yr with $3$ months, and $9\cdot10^{25}$ yr with $5$ years of live time. Under the same conditions, the discovery sensitivity after $3$ months and $5$ years will be $7\cdot10^{24}$ yr and $4\cdot10^{25}$ yr, respectively.

Tables

Table 1. Table 1: Input parameters for the production of toy-MC spectra.

$BI$ [cts$/($keV$\cdot$kg$\cdot$yr$)$]	$R^{Co}$ [cts$/($kg$\cdot$yr$)$]
$\left(1.02\pm 0.03{\left(\text{stat}\right)}_{-0.10}^{+0.23}\left(\text{syst}\right)\right)\cdot {10}^{\text{-}2}$	$0.428$

Table 2. Table 2: Crystal subsets with different expected α 𝛼 \alpha background in CUORE. The values of B ​ I 𝐵 𝐼 BI and R C ​ o superscript 𝑅 𝐶 𝑜 R^{Co} are taken from  Alduino:2017qet .

Subset	Free	Number of	$BI$ [cts$/$	$R^{Co}$ [cts
Name	Sides	crystals	$($keV$\cdot$kg$\cdot$yr$)$]	$/($kg$\cdot$yr$)$]
Inner	$0$	$528$	$0.82\left(2\right)\cdot {10}^{\text{-}2}$	$0.40$
Middle-1	$1$	$272$	$1.17\left(4\right)\cdot {10}^{\text{-}2}$	$0.47$
Middle-2	$2$	$164$	$1.36\left(4\right)\cdot {10}^{\text{-}2}$	$0.43$
Outer	$3$	$24$	$1.78\left(7\right)\cdot {10}^{\text{-}2}$	$0.59$

Table 3. Table 3: Constants used in H 0 subscript 𝐻 0 H_{0} and H 1 subscript 𝐻 1 H_{1} .

Constant	Symbol	Value
Detector Mass	$m_d$	$741.67$ kg
Avogadro number	$N_A$	$6.022\cdot {10}^{23}$ mol-1
Molar mass	$m_A$	$159.6$ g$/$mol
Live Time	$t_d$	$0.1$–$5$ yr
Efficiency	${\epsilon }_{tot}$	$81.3\%$
130Te abundance	$f_{130}$	$0.34167$
$0\nu \beta \beta$ Q-value	$Q_{\beta \beta }$	$2527.518$ keV
60Co peak position	${\mu }_{Co}$	$\left(2505.692+1.9\right)$ keV
Energy resolution	FWHM	$5,10$ keV

Table 4. Table 4: Median exclusion sensitivity for different energy resolutions and different subset numbers.

FWHM	$N_d$	${\hat{T}}_{1/2}^{\mathrm{\hspace{0.17em}0}\nu }$ at $0.25$ yr	${\hat{T}}_{1/2}^{\mathrm{\hspace{0.17em}0}\nu }$ at $5$ yr
$[$keV$]$		[yr]	[yr]
$5$	$1$	$1.7\cdot {10}^{25}$	$8.9\cdot {10}^{25}$
$10$	$1$	$1.2\cdot {10}^{25}$	$6.1\cdot {10}^{25}$
$5$	$4$	$1.8\cdot {10}^{25}$	$9.1\cdot {10}^{25}$
$10$	$4$	$1.3\cdot {10}^{25}$	$6.2\cdot {10}^{25}$

Equations

$${\Gamma }^{0\nu }=\frac{\ln 2}{T_{1/2}^{0\nu }}.$$

$$P\left(\Gamma^{0\nu},\vec{\theta}\big{|}\vec{E},H_{1}\right)=\\ \frac{P\left(\vec{E}\big{|}\Gamma^{0\nu},\vec{\theta},H_{1}\right)\pi(\Gamma^{0\nu})\pi(\vec{\theta})}{\bigintsss_{\Omega}\bigintsss_{0}^{\infty}P\left(\vec{E}\big{|}\Gamma^{0\nu},\vec{\theta},H_{1}\right)\pi(\Gamma^{0\nu})\pi(\vec{\theta})\,d\vec{\theta}\,d\Gamma^{0\nu}}\ .$$

$$P\left(\Gamma^{0\nu}\big{|}H_{1},\vec{E}\right)=\bigintsss_{\Omega}P\left(\Gamma^{0\nu},\vec{\theta}\big{|}\vec{E},H_{1}\right)d\vec{\theta}\ .$$

$$P\left(\vec{E}\big{|}\Gamma^{0\nu},\vec{\theta},H\right)=\mathcal{L}\left(\vec{E}\big{|}\Gamma^{0\nu},\vec{\theta},H\right)\\ =\prod_{d=1}^{N_{d}}\ \prod_{i=1}^{N}\frac{e^{-\lambda_{di}}\cdot\lambda_{di}^{n_{di}}}{n_{di}!}\ ,$$

$$\ln{\mathcal{L}\left(\vec{E}\big{|}\Gamma^{0\nu},\vec{\theta},H\right)}=\sum_{d=1}^{N_{d}}\ \sum_{i=1}^{N}\left(-\lambda_{di}+\ln{\lambda_{di}^{n_{di}}}\right)\ ,$$

$$f_{bkg}(E)=\frac{1}{E_{max}-E_{min}}$$

$$f_{Co}(E)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{\left(E-{\mu }_{Co}\right)}^2}{2{\sigma }^2}\right).$$

$${\lambda }_{di}^{bkg}={\text{bigintsss}}_{E_{di}^{min}}^{E_{di}^{max}}{M}_d^{bkg}{f}_{bkg}(E)dE$$

$${\lambda }_{di}^{bkg}=B{I}_d\cdot m_d\cdot t_d\cdot \delta E_{di}.$$

$${\lambda }_{di}^{Co}={\text{bigintss}}_{E_{di}^{min}}^{E_{di}^{max}}\frac{M_d^{Co}}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{\left(E-{\mu }_{Co}\right)}^2}{2{\sigma }^2}\right)dE,$$

$$M_d^{Co}=R_d^{Co}\cdot m_d\cdot t_d.$$

$${\lambda }_{di}={\lambda }_{di}^{bkg}+{\lambda }_{di}^{Co}.$$

$${\lambda }_{di}^{0\nu }={\text{bigintss}}_{E_{di}^{min}}^{E_{di}^{max}}\frac{M_d^{0\nu }}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{\left(E-Q_{\beta \beta }\right)}^2}{2{\sigma }^2}\right)dE.$$

$$M_d^{0\nu }={\Gamma }^{0\nu }\frac{N_A}{m_A}\cdot f_{130}\cdot {\epsilon }_{tot}\cdot m_d\cdot t_d,$$

$$T_{1/2}^{\,0\nu}\left(90\%\ CI\right)=T_{1/2}^{\,0\nu}:\\ \bigintsss_{0}^{\ln{2}/T_{1/2}^{\,0\nu}}P\left(\Gamma^{0\nu}\big{|}H_{1},\vec{E}\right)d\Gamma^{0\nu}=0.9\ .$$

$$P\bigl{(}H_{0}\big{|}\vec{E}\bigr{)}+P\bigl{(}H_{1}\big{|}\vec{E}\bigr{)}=1\ .$$

$$\displaystyle\begin{split}P\bigl{(}H_{0}\big{|}\vec{E}\bigr{)}&=\frac{P\bigl{(}\vec{E}\big{|}H_{0}\bigr{)}\pi(H_{0})}{P\bigl{(}\vec{E}\bigr{)}}\\ P\bigl{(}H_{1}\big{|}\vec{E}\bigr{)}&=\frac{P\bigl{(}\vec{E}\big{|}H_{1}\bigr{)}\pi(H_{1})}{P\bigl{(}\vec{E}\bigr{)}}\ .\end{split}$$

$$\begin{split}P\bigl{(}\vec{E}\big{|}H_{0}\bigr{)}=\bigintsss_{\Omega}P\left(\vec{E}\big{|}\vec{\theta},H_{0}\right)\pi(\vec{\theta})\,d\vec{\theta}\\ P\bigl{(}\vec{E}\big{|}H_{1}\bigr{)}=\bigintsss_{\Omega}\bigintsss_{0}^{\infty}P\left(\vec{E}\big{|}\Gamma^{0\nu},\vec{\theta},H_{1}\right)\times\quad\quad\\ \times\pi(\Gamma^{0\nu})\pi(\vec{\theta})\,d\vec{\theta}\,d\Gamma^{0\nu}\ ,\end{split}$$

$$P\bigl{(}\vec{E}\bigr{)}=P\bigl{(}\vec{E}\big{|}H_{0}\bigr{)}\pi(H_{0})+P\bigl{(}\vec{E}\big{|}H_{1}\bigr{)}\pi(H_{1})\ .$$

$$\hat{T}_{1/2}^{0\nu}(3\sigma)=T_{1/2}^{0\nu}:P\bigl{(}H_{0}\big{|}\vec{E}\bigr{)}<0.0027\\ \text{~{}for~{}}N/2~{}\text{experiments.}$$

$$r_{bkg(Co)}=\frac{M^{bkg(Co)}-N^{bkg(Co)}}{{\sigma }_{M^{bkg(Co)}}},$$

Full text

If the background level and the energy resolution meet the expectation, CUORE will reach a $90\%$ CI exclusion sensitivity of $2\cdot 10^{25}$ yr with $3$ months, and $9\cdot 10^{25}$ yr with $5$ years of live time. Under the same conditions, the discovery sensitivity after $3$ months and $5$ years will be $7\cdot 10^{24}$ yr and $4\cdot 10^{25}$ yr, respectively.

††journal: Eur. Phys. J. C

1 Introduction

Neutrinoless double beta decay is a non Standard Model process that violates the total lepton number conservation and implies a Majorana neutrino mass component Schechter:1981bd ; Duerr:2011zd . This decay is currently being investigated with a variety of double beta decaying isotopes. A recent review can be found in Ref. Dell'Oro:2016dbc . The Cryogenic Underground Observatory for Rare Events (CUORE) Artusa:2014lgv ; CUORE-NIM2004 ; Arnaboldi:2003tu is an experiment searching for $0\nu\beta\beta$ decay in 130Te. It is located at the Laboratori Nazionali del Gran Sasso of INFN, Italy. In CUORE, $988$ TeO2 crystals with natural 130Te isotopic abundance and a $750$ g average mass are operated simultaneously as source and bolometric detector for the decay. In this way, the $0\nu\beta\beta$ decay signature is a peak at the $Q$-value of the reaction ($Q_{\beta\beta}$, $2527.518$ keV for 130Te Redshaw:2009cf ; Scielzo:2009nh ; Rahaman:2011zz ). Bolometric crystals are characterized by an excellent energy resolution ($\sim 0.2\%$ Full Width at Half Maximum, FWHM) and a very low background at $Q_{\beta\beta}$, which is expected to be at the $10^{\mbox{-}2}$ cts$/($keV$\cdot$kg$\cdot$yr$)$ level in CUORE Alduino:2017qet .

The current best limit on $0\nu\beta\beta$ decay in 130Te comes from a combined analysis of the CUORE-0 Alduino:2016vjd ; Aguirre:2014lua and Cuoricino data Arnaboldi:2008ds ; Andreotti:2010vj . With a total exposure of $29.6$ kg$\cdot$yr, a limit of $T_{1/2}^{0\nu}>4.0\cdot 10^{24}$ yr ($90\%$ CI) is obtained Alfonso:2015wkao for the $0\nu\beta\beta$ decay half life, $T_{1/2}^{\,0\nu}$.

After the installation of the detector, successfully completed in the summer 2016, CUORE started the commissioning phase at the beginning of 2017. The knowledge of the discovery and exclusion sensitivity to $0\nu\beta\beta$ decay as a function of the measurement live time can be exploited to set the criteria for the unblinding of the data and the release of the $0\nu\beta\beta$ decay analysis results.

In this work, we dedicate our attention to those factors which could strongly affect the sensitivity, such as the Background Index ($BI$) and the energy resolution at $Q_{\beta\beta}$. In CUORE, the crystals in the outer part of the array are expected to show a higher $BI$ than those in the middle Alduino:2017qet . Considering this and following the strategy already implemented by the Gerda Collaboration Agostini:2013mzu ; Agostini:2017iyd , we show how the division of the data into subsets with different $BI$ could improve the sensitivity.

The reported results are obtained by means of a Bayesian analysis performed with the Bayesian Analysis Toolkit (BAT) Caldwell:2008fw . The analysis is based on a toy-MC approach. At a cost of a much longer computation time with respect to the use of the median sensitivity formula Alessandria:2011rc , this provides the full sensitivity probability distribution and not only its median value.

In Section 2, we review the statistical methods for the parameter estimation, as well as for the extraction of the exclusion and discovery sensitivity. Section 3 describes the experimental parameters used for the analysis while its technical implementation is summarized in Section 4. Finally, we present the results in Section 5.

2 Statistical Method

The computation of exclusion and discovery sensitivities presented here follows a Bayesian approach: we exploit the Bayes theorem both for parameter estimation and model comparison. In this work, we use the following notation:

• •

$H$ indicates both a hypothesis and the corresponding model;

• •

$H_{0}$ is the background-only hypothesis, according to which the known physics processes are enough to explain the experimental data. In the present case, we expect the CUORE background to be flat in a $100$ keV region around $Q_{\beta\beta}$, except for the presence of a 60Co summation peak at $2505.7$ keV. Therefore, $H_{0}$ is implemented as a flat background distribution plus a Gaussian describing the 60Co peak. In CUORE-0, this peak was found to be centered at an energy $1.9\pm 0.7$ keV higher than that tabulated in literature Alfonso:2015wkao . This effect, present also in Cuoricino Andreotti:2010vj , is a feature of all gamma summation peaks. Hence, we will consider the 60Co peak to be at $2507.6$ keV.

• •

$H_{1}$ is the background-plus-signal hypothesis, for which some new physics is required to explain the data. In our case, the physics involved in $H_{1}$ contains the background processes as well as $0\nu\beta\beta$ decay. The latter is modeled as a Gaussian peak at $Q_{\beta\beta}$.

• •

$\vec{E}$ represents the data. It is a list of $N$ energy bins centered at the energy $E_{i}$ and containing $n_{i}$ event counts. The energy range is $[2470;2570]$ keV. This is the same range used for the CUORE-0 $0\nu\beta\beta$ decay analysis Alfonso:2015wkao , and is bounded by the possible presence of peaks from 214Bi at $2447.7$ keV and 208Tl X-ray escape at $\sim 2585$ keV Alfonso:2015wkao . While an unbinned fit allows to fully exploit the information contained in the data, it can result in a long computation time for large data samples. Given an energy resolution of $\sim 5$ keV FWHM and using a $1$ keV bin width, the $\pm 3$ sigma range of a Gaussian peak is contained in $12.7$ bins. With the $1$ keV binning choice, the loss of information with respect to the unbinned fit is negligible.

• •

$\Gamma^{0\nu}$ is the parameter describing the $0\nu\beta\beta$ decay rate for $H_{1}$:

[TABLE]

• •

$\vec{\theta}$ is the list of nuisance parameters describing the background processes in both $H_{0}$ and $H_{1}$;

• •

$\Omega$ is the parameter space for the parameters $\vec{\theta}$.

2.1 Parameter Estimation

We perform the parameter estimation for a model $H$ through the Bayes theorem, which yields the probability distribution for the parameters based on the measured data, under the assumption that the model $H$ is correct. In the $0\nu\beta\beta$ decay analysis, we are interested in the measurement of $\Gamma^{0\nu}$ for the hypothesis $H_{1}$. The probability distribution for the parameter set $(\Gamma^{0\nu},\vec{\theta})$ is:

[TABLE]

The numerator contains the conditional probability

P\left(\vec{E}\big{|}\Gamma^{0\nu},\vec{\theta},H_{1}\right) of finding the measured data $\vec{E}$ given the model $H_{1}$ for a set of parameters $(\Gamma^{0\nu},\vec{\theta})$, times the prior probability $\pi$ for each of the considered parameters. The prior probability has to be chosen according to the knowledge available before the analysis of the current data. For instance, the prior for the number of signal counts $\Gamma^{0\nu}$ might be based on the half-life limits reported by previous experiments while the prior for the background level in the region of interest (ROI) could be set based on the extrapolation of the background measured outside the ROI. The denominator represents the overall probability to obtain the data $\vec{E}$ given the hypothesis $H_{1}$ and all possible parameter combinations, $P(\vec{E}|H_{1})$.

The posterior probability distribution for $\Gamma^{0\nu}$ is obtained via marginalization, i.e. integrating

P\left(\Gamma^{0\nu},\vec{\theta}\big{|}\vec{E},H_{1}\right) over all nuisance parameters $\vec{\theta}$:

[TABLE]

For each model $H$, the probability of the data given the model and the parameters has to be defined. For a fixed set of experimental data, this corresponds to the likelihood function james . Dividing the data into $N_{d}$ subsets with index $d$ characterized by different background levels, and considering a binned energy spectrum with $N$ bins and a number $n_{di}$ of events in the bin $i$ of the $d$ subset spectrum, the likelihood function is expressed by the product of a Poisson term for each bin $di$:

[TABLE]

where $\lambda_{di}$ is the expectation value for the bin $di$. The best-fit is defined as the set of parameter values $(\Gamma^{0\nu},\vec{\theta})$ for which the likelihood is at its global maximum. In the practical case, we perform the maximization on the log-likelihood

[TABLE]

where the additive terms $-\ln{(n_{di}!)}$ are dropped from the calculation.

The difference between $H_{0}$ and $H_{1}$ is manifested in the formulation of $\lambda_{di}$. As mentioned above, we parametrize $H_{0}$ with a flat distribution over the considered energy range, i.e. $[2470;2570]$ keV:

[TABLE]

plus a Gaussian distribution for the 60Co peak:

[TABLE]

The expected background counts in the bin $di$ corresponds to the integral of $f_{bkg}(E)$ in the bin $di$ times the total number of background counts $M^{bkg}_{d}$ for the subset $d$:

[TABLE]

where $E^{min}_{di}$ and $E^{max}_{di}$ are the left and right margins of the energy bin $di$, respectively. Considering bins of size $\delta E_{di}$ and expressing $M^{bkg}_{di}$ as function of the background index $BI_{d}$, of the total mass $m_{d}$ and of the measurement live time $t_{d}$, we obtain:

[TABLE]

Similarly, the expectation value for the 60Co distribution on the bin $di$ is:

[TABLE]

where $M^{Co}_{d}$ is the total number of 60Co events for the subset $d$ and can be redefined as function of the 60Co event rate, $R^{Co}_{d}$:

[TABLE]

The total expectation value $\lambda_{di}$ for $H_{0}$ is then:

[TABLE]

In the case of $H_{1}$ an additional expectation value for $0\nu\beta\beta$ decay is required:

[TABLE]

The number of $0\nu\beta\beta$ decay events in the subset $d$ is:

[TABLE]

where $N_{A}$ is the Avogadro number, $m_{a}$ and $f_{130}$ are the molar mass and the isotopic abundance of 130Te and $\varepsilon_{tot}$ is the total efficiency, i.e. the product of the containment efficiency $\varepsilon_{MC}$ (obtained with MC simulations) and the instrumental efficiency $\varepsilon_{instr}$.

2.2 Exclusion Sensitivity

We compute the exclusion sensitivity by means of the $90\%$ CI limit. This is defined as the value of $T_{1/2}^{\,0\nu}$ corresponding to the $90\%$ quantile of the posterior $\Gamma^{0\nu}$ distribution:

[TABLE]

An example of posterior probability for $\Gamma^{0\nu}$ and the relative $90\%$ CI limit is shown in Fig. 1, top. Flat prior distributions are used for all parameters, as described in Sec. 3.

In the Bayesian approach, the limit is a statement regarding the true value of the considered physical quantity. In our case, a $90\%$ CI limit on $T_{1/2}^{\,0\nu}$ is to be interpreted as the value above which, given the current knowledge, the true value of $T_{1/2}^{\,0\nu}$ lies with $90\%$ probability. This differs from a frequentist $90\%$ C.L. limit, which is a statement regarding the possible results of the repetition of identical measurements and should be interpreted as the value above which the best-fit value of $T_{1/2}^{\,0\nu}$ would lie in the $90\%$ of the imaginary identical experiments.

In order to extract the exclusion sensitivity, we generate a set of $N$ toy-MC spectra according to the back-ground-only model, $H_{0}$. We then fit spectra with the background-plus-signal model, $H_{1}$, and obtain the

$T_{1/2}^{\,0\nu}\left(90\%\ CI\right)$ distribution (Fig. 1, bottom). Its median $\hat{T}_{1/2}^{\,0\nu}\left(90\%\ CI\right)$ is referred as the median sensitivity. For a real experiment, the experimental $T_{1/2}^{\,0\nu}$ limit is expected to be above/below $\hat{T}_{1/2}^{\,0\nu}\left(90\%\ CI\right)$ with $50\%$ probability. Alternatively, one can consider the mode of the distribution, which corresponds to the most probable $T_{1/2}^{\,0\nu}$ limit.

The exact procedure for the computation of the exclusion sensitivity is the following:

• •

for each subset, we generate a random number of background events $N_{d}^{bkg}$ according to a Poisson distribution with mean $\lambda^{bkg}_{d}$;

• •

for each subset, we generate $N_{d}^{bkg}$ events with an energy randomly distributed according to $f_{bkg}(E)$;

• •

we repeate the procedure for the 60Co contribution;

• •

we fit the toy-MC spectrum with the $H_{1}$ model (Eq. 2), and marginalize the likelihood with respect to the parameters $BI_{d}$ and $R_{d}^{Co}$ (Eq. 3);

• •

we extract the $90\%$ CI limit on $T_{1/2}^{\,0\nu}$;

• •

we repeat the algorithm for $N$ toy-MC experiments, and build the distribution of $T_{1/2}^{\,0\nu}\left(90\%\ CI\right)$.

2.3 Discovery Sensitivity

The discovery sensitivity provides information on the required strength of the signal amplitude for claiming that the known processes alone are not sufficient to properly describe the experimental data. It is computed on the basis of the comparison between the background-only and the background-plus-signal models. A method for the calculation of the Bayesian discovery sensitivity was introduced in Ref. Caldwell:2006yj . We report it here for completeness.

In our case, we assume that $H_{0}$ and $H_{1}$ are a complete set of models, for which:

[TABLE]

The application of the Bayes theorem to the models $H_{0}$ and $H_{1}$ yields:

[TABLE]

In this case, the numerator contains the probability of measuring the data $\vec{E}$ given the model $H$:

[TABLE]

while the prior probabilities for the models $H_{0}$ and $H_{1}$ can be chosen as $0.5$ so that neither model is favored.

The denominator of Eq. 17 is the sum probability of obtaining the data $\vec{E}$ given either the model $H_{0}$ or $H_{1}$:

[TABLE]

At this point we need to define a criterion for claiming the discovery of new physics. Our choice is to quote the $3~{}\sigma$ (median) discovery sensitivity, i.e. the value of $T_{1/2}^{\,0\nu}$ for which the posterior probability of the back-ground-only model $H_{0}$ given the data is smaller than $0.0027$ in $50\%$ of the possible experiments. In other words:

[TABLE]

The detailed procedure for the determination of the discovery sensitivity is:

• •

we produce a toy-MC spectrum according to the $H_{1}$ model with an arbitrary value of $T_{1/2}^{\,0\nu}$;

• •

we fit the spectrum with both $H_{0}$ and $H_{1}$;

• •

we compute $P(H_{0}|\vec{E})$;

• •

we repeat the procedure for $N$ toy-MC spectra using the same $T_{1/2}^{\,0\nu}$;

• •

we repeat the routine with different values of $T_{1/2}^{\,0\nu}$ until the condition of Eq. 20 is satisfied. The iteration is implemented using the bisection method until a $5\cdot 10^{\mbox{-}5}$ precision is obtained on the median $P(H_{0}|\vec{E})$.

3 Experimental Parameters

The fit parameters of the $H_{1}$ model are $BI$, $R^{Co}$ and $\Gamma^{0\nu}$, while only the first two are present for $H_{0}$. If the data are divided in subsets, different $BI$ and $R^{Co}$ fit parameter are considered for each subset. On the contrary, the inverse $0\nu\beta\beta$ half-life is common to all subsets.

Prior to the assembly of the CUORE crystal towers, we performed a screening survey of the employed materials Alessandria:2011vj ; Barghouty:2010kj ; Wang:2015pxa ; Alessandria:2012zp ; Andreotti:2009dk ; Bellini:2009zw ; Andreotti:2009zza ; giachero . From these measurements, either a non-zero activity was obtained, or a $90\%$ confidence level (C.L.) upper limit was set. Additionally, the radioactive contamination of the crystals and holders was also obtained from the CUORE-0 background model Alduino:2016vtd . We developed a full MC simulation of CUORE Alduino:2017qet , and we used the results of the screening measurements and of the CUORE-0 background model for the normalization of the simulated spectra. We then computed the $BI$ at $Q_{\beta\beta}$ using the output of the simulations. In the present study, we consider only those background contributions for which a non-zero activity is obtained from the available measurements. The largest background consists of $\alpha$ particles emitted by U and Th surface contaminations of the copper structure holding the crystals. Additionally, we consider a 60Co contribution normalized to the $90\%$ C.L. limit from the screening measurement. In this sense, the effect of a 60Co background on the CUORE sensitivity is to be held as an upper limit. Given the 60Co importance especially in case of sub-optimal energy resolution, we preferred to maintain a conservative approach in this regard. In the generation of the toy-MC spectra, we take into account the 60Co half life ($5.27$ yr), and set the start of data taking to January $2017$.

The parameter values used for the production of the toy-MC are reported in Tab. 1. The quoted uncertainty on the $BI$ comes from the CUORE MC simulations Alduino:2017qet . We produce the toy-MC spectra using the best-fit value of the $BI$. In a second time, we repeat the analysis after increasing and decreasing the $BI$ by an amount equivalent to its statistical and systematic uncertainties combined in quadrature.

After running the fit on the entire crystal array as if it were a unique detector, we considered the possibility of dividing the data grouping the crystals with a similar $BI$. Namely, being the background at $Q_{\beta\beta}$ dominated by surface $\alpha$ contamination of the copper structure, the crystals facing a larger copper surface are expected to have a larger $BI$. This effect was already observed in CUORE-0, where the crystals in the uppermost and lowermost levels, which had $3$ sides facing the copper shield, were characterized by a larger background than those in all other levels, which were exposed to coppper only on $2$ sides. Considering the CUORE geometry, the crystals can be divided in $4$ subsets with different numbers of exposed faces. Correspondingly, they are characterized by different $BI$, as reported in Tab. 2.

A major ingredient of a Bayesian analysis is the choice of the priors. In the present case, we use a flat prior for all parameters. In particular, the prior distribution for $\Gamma^{0\nu}$ is flat between zero and a value large enough to contain $>99.999\%$ of its posterior distribution. This corresponds to the most conservative choice. Any other reasonable prior, e.g. a scale invariant prior on $\Gamma^{0\nu}$, would yield a stronger limit. A different prior choice based on the real characteristic of the experimental spectra might be more appropriate for $BI$ and $R^{Co}$ in the analysis of the CUORE data. For the time being the lack of data prevents the use of informative priors. As a cross-check, we performed the analysis using the $BI$ and 60Co rate uncertainties obtained by the background budget as the $\sigma$ of a Gaussian prior. No significant difference was found on the sensitivity band because the Poisson fluctuations of the generated number of background and 60Co events are dominant for the extraction of the $\Gamma^{0\nu}$ posterior probability distribution.

Tab. 3 lists the constant quantities present in the formulation of $H_{0}$ and $H_{1}$. All of them are fixed, with the exception of the live time $t$ and the FWHM of the $0\nu\beta\beta$ decay and 60Co Gaussian peaks. We perform the analysis with a FWHM of $5$ and $10$ keV, corresponding to a $\sigma$ of $2.12$ and $4.25$ keV, respectively. Regarding the efficiency, while in the toy-MC production the $BI$ and $R^{Co}$ are multiplied by the instrumental efficiency111The containment efficiency is already encompassed in $BI$ and $R^{Co}$ Alduino:2017qet ., in the fit the total efficiency is used. This is the product of the containment and instrumental efficiency. Also in this case, we use the same value as for CUORE-0, i.e. $81.3\%$ Alfonso:2015wkao . We evaluate the exclusion and discovery sensitivities for different live times, with $t$ ranging from $0.1$ to $5$ yr and using logarithmically increasing values: $t_{i}=1.05\cdot t_{i-1}$.

4 Fit Procedure

We perform the analysis with the software BAT v1.1.0-DEV Caldwell:2006yj , which internally uses CUBA Hahn:2004fe v4.2 for the integration of multi-dimensional probabilities and the Metropolis-Hastings algorithm Metropolis for the fit. The computation time depends on the number of samples drawn from the considered probability distribution. For the exclusion sensitivity, we draw $10^{5}$ likelihood samples for every toy-MC experiment, while, due to the higher computational cost, we use only $10^{3}$ for the discovery sensitivity. For every combination of live time, $BI$ and energy resolution, we run $10^{5}$ ($10^{3}$) toy-MC experiments for the exclusion (discovery) sensitivity study. In the case of the discovery sensitivity, we chose the number of toy-MC experiments as the minimum for which a $2\%$ relative precision was achievable on the median sensitivity. For the exclusion sensitivity, it was possible to increase both the number of toy-MC experiments and iterations, with a systematic uncertainty on the median sensitivity at the per mil level.

5 Results and Discussion

5.1 Exclusion Sensitivity

The distributions of $90\%$ CI limit as a function of live time with no data subdivision are shown in Fig. 2. For all $BI$ values and all live times, the FWHM of $5$ keV yields a $\sim 45\%$ higher sensitivity with respect to a $10$ keV resolution. The median sensitivity after $3$ months and $5$ years of data collection in the two considered cases are reported in Tab. 4. The dependence of the median sensitivity on live time is typical of a background-dominated experiments: namely, CUORE expects about one event every four days in a $\pm 3\sigma$ region around $Q_{\beta\beta}$. The results in Tab. 4 show also the importance of energy resolution and suggest to put a strong effort in its optimization. As a cross check, we compare the sensitivity just obtained with that provided by the analytical method presented in Alessandria:2011rc and shown in dark green in Fig. 2. The analytical method yields a slightly higher sentitivity for short live times, while the two techniques agree when the data sample is bigger. We justify this with the fact that the uncertainty on the number of background counts obtained with the Bayesian fit is slightly larger than the corresponding Poisson uncertainty assumed in the analytical approach Cowan:2011an , hence the limit on $T_{1/2}^{\,0\nu}$ is systematically weaker222See the discussion of the pulls for a more detailed explanation.. The effect becomes less and less strong with increasing data samples, i.e. with growing live time. With a resolution of $5$ keV, the difference goes from $8\%$ after $3$ months to $<0.1\%$ after $5$ years, while for a $10$ keV FWHM the difference is $\sim 6\%$ after $3$ months and $4\%$ after $5$ years. One remark has to be done concerning the values reported in Alessandria:2011rc : there we gave a $90\%$ C.I. exclusion sensitivity of $9.3\cdot 10^{25}$ yr with $5$ yr of live time. This is $\sim 5\%$ higher than the result presented here and is explained by the use of a different total efficiency: $87.4\%$ in Alessandria:2011rc and $81.3\%$ in this work.

We then extract the exclusion sensitivity after dividing the crystals in $4$ subsets, as described in Sec. 3. The median exclusion sensitivity values after $3$ months and $5$ years of data collection with one and $4$ subsets are reported in Tab. 4. The division in subsets yields only a small improvement (at the percent level) in median sensitivity. Based on this results only, one would conclude that dividing the data into subsets with different $BI$ is not worth the effort. This conclusion is not always true, and strongly relies on the exposure and $BI$ of the considered subsets. As an example, we repeated a toy analysis assuming a $BI$ of $10^{\mbox{-}2}$ cts$/($keV$\cdot$kg$\cdot$yr$)$, and with two subsets of equal exposure and $BI$ $0.5\cdot 10^{\mbox{-}2}$ cts$/($keV$\cdot$kg$\cdot$yr$)$ and $1.5\cdot 10^{\mbox{-}2}$ cts$/($keV$\cdot$kg$\cdot$yr$)$, respectively. In this case, the division of the data in to two subsets yields a $\sim 10\%$ improvement after $5$ yr of data taking. Hence, the data subdivision is a viable option for the final analysis, whose gain strongly depends on the experimental BI of each channel. Similarly, we expect the CUORE bolometers to have different energy resolutions; in CUORE-0, these ranged from $\sim 3$ keV to $\sim 20$ keV FWHM Alduino:2016zrl . In the real CUORE analysis a further splitting of the data can be done by grouping the channels with similar FWHM, or by keeping every channels separate. At the present stage it is not possible to make reliable predictions for the FWHM distribution among the crystals, so we assumed an average value (of $5$ or $10$ keV) throughout the whole work.

Ideally, the final CUORE $0\nu\beta\beta$ decay analysis should be performed keeping the spectra collected by each crystal separate, additionally to the usual division of the data into data sets comprised by two calibration runs Alfonso:2015wkao . Assuming an average frequency of one calibration per month, the total number of energy spectra would be $\sim 6\cdot 10^{4}$. Assuming a different but stationary $BI$ for each crystal, and using the same 60Co rate for all crystals, the fit model would have $\sim 10^{3}$ parameters. This represents a major obstacle for any existing implementation of the Metropolis-Hastings or Gibbs sampling algorithm. A possible way to address the problem might be the use of different algorithms, e.g. nested sampling Feroz:2008xx ; Handley:2015fda , or a partial analytical solution of the likelihood maximization.

We perform two further cross-checks in order to investigate the relative importance of the flat background and the 60Co peak. In the first scenario we set the $BI$ to zero, and do the same for the 60Co rate in the second one. In both cases, the data are not divided into subsets, and resolutions of $5$ and $10$ keV are considered. With no flat background and a $5$ keV resolution, no 60Co event leaks in the $\pm 3\sigma$ region around $Q_{\beta\beta}$ even after $5$ yr of measurement. As a consequence, the $90\%$ CI limits are distributed on a very narrow band, and the median sensitivity reaches $1.2\cdot 10^{27}$ yr after $5$ yr of data collection. On the contrary, if we assume a $10$ keV FWHM, some 60Co events fall in the $0\nu\beta\beta$ decay ROI from the very beginning of the data taking. This results in a strong asymmetry of the sensitivity band. In the second cross-check, we keep the $BI$ at $1.02\cdot 10^{\mbox{-}2}$ cts$/($keV$\cdot$kg$\cdot$yr$)$, but set the 60Co rate to zero. In both cases, the difference with respect to the standard scenario is below $1\%$. We can conclude that the 60Co peak with an initial rate of $0.428$ cts/(kg$\cdot$yr) is not worrisome for a resolution of up to $10$ keV, and that the lower sensitivity obtained with $10$ keV FWHM with respect to the $5$ keV case is ascribable to the relative amplitude of $\lambda^{bkg}_{di}$ and $\lambda^{0\nu}_{di}$ only (Eqs. 9 and 13). This is also confirmed by the computation of the sensitivity for the optimistic scenario without the $1.9$ keV shift of the 60Co peak used in the standard case.

We test the fit correctness and bias computing the pulls, i.e. the normalized residuals, of the number of counts assigned to each of the fit components. Denoting with $N^{bkg}$ and $N^{Co}$ the number of generated background and 60Co events, respectively, and with $M^{bkg}$ and $M^{Co}$ the corresponding number of reconstructed events, the pulls are defined as:

[TABLE]

where $\sigma_{M^{bkg(Co)}}$ is the statistical uncertainty on $M^{bkg(Co)}$ given by the fit.

For an unbiased fit, the distribution of the pulls is expected to be Gaussian with a unitary root mean square (RMS). In the case of exclusion sensitivity, we obtain $r_{bkg}=-0.2\pm 0.4$ and $r_{Co}=0.1\pm 0.5$ for all live times. The fact that the pull distributions are slightly shifted indicates the presence of a bias. Its origin lies in the Bayesian nature of the fit and is that all fit contributions are constrained to be greater than zero. We perform a cross-check, by extending the range of all parameters ($BI$, $R^{Co}$ and $\Gamma^{0\nu}$) to negative values. Under this condition, the bias disappears. In addition to this, an explanation is needed for the small RMS of the pull distributions. This is mainly due to two effects: first, the toy-MC spectra are generated using $H_{0}$, while the fit is performed using $H_{1}$; second, the statistical uncertainties on all parameters are larger than the Poisson uncertainty on the number of generated events. To confirm the first statement, we repeat the fit using $H_{0}$ instead of $H_{1}$ and we obtain pulls with zero mean and an RMS $\sim 0.8$, which is closer to the expected value. Finally, we compare the parameter uncertainty obtained from the fit with the Poisson uncertainty for the equivalent number of counts, and we find that the difference is of $O(20\%)$.

5.2 Discovery Sensitivity

The extraction of the discovery sensitivity involves fits with the background-only and the background-plus-signal models. Moreover, two multi-dimensional integrations have to be performed for each toy-MC spectrum, and a loop over the $0\nu\beta\beta$ decay half-life has to be done until the condition of Eq. 20 is met. Due to the high computation cost, we compute the $3~{}\sigma$ discovery sensitivity for a FWHM of $5$ and $10$ keV with no crystal subdivision. As shown in Fig. 3, with a $5$ keV energy resolution CUORE has a $3~{}\sigma$ discovery sensitivity superior to the limit obtained from the combined analysis of Cuore-0 and Cuoricino data Alfonso:2015wkao after less than one month of operation, and reaches $3.7\cdot 10^{25}$ yr with $5$ yr of live time.

Also in this case, the pulls are characterized by an RMS smaller than expected, but no bias is present due to the use of $H_{1}$ for both the generation and the fit of the toy-MC spectra.

6 Conclusion and Outlook

We implemented a toy-MC method for the computation of the exclusion and discovery sensitivity of CUORE using a Bayesian analysis. We have highlighted the influence of the $BI$ and energy resolution on the exclusion sensitivity, showing how the achievement of the expected $5$ keV FWHM is desirable. Additionally, we have shown how the division of the data into subsets with different $BI$ could yield an improvement in exclusion sensitivity.

Once the CUORE data collection starts and the experimental parameters are available, the sensitivity study can be repeated in a more detailed way. As an example, non-Gaussian spectral shapes for the $0\nu\beta\beta$ decay and 60Co peaks can be used, and the systematics of the energy reconstruction can be included.

Acknowledgments

The CUORE Collaboration thanks the directors and staff of the Laboratori Nazionali del Gran Sasso and the technical staff of our laboratories. CUORE is 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-1404205; The Alfred P. Sloan Foundation; The University of Wisconsin Foundation; Yale University; The US Department of Energy (DOE) Office of Science under Contract Nos. DE-AC02-05CH1-1231, DE-AC52-07NA27344, and DE-SC0012654; The DOE Office of Science, Office of Nuclear Physics under Contract Nos. DE-FG02-08ER41551 and DE-FG03-00ER41138; The National Energy Research Scientific Computing Center (NERSC).