CUORE sensitivity to 0νββ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0\nu \beta \beta $$\end{document} decay

We report a study of the CUORE sensitivity to neutrinoless double beta (0νββ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0\nu \beta \beta $$\end{document}) decay. We used a Bayesian analysis based on a toy Monte Carlo (MC) approach to extract the exclusion sensitivity to the 0νββ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0\nu \beta \beta $$\end{document} decay half-life (T1/20ν\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{1/2}^{\,0\nu }$$\end{document}) at 90%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$90\%$$\end{document} credibility interval (CI) – i.e. the interval containing the true value of T1/20ν\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{1/2}^{\,0\nu }$$\end{document} with 90%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$90\%$$\end{document} probability – and the 3σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3~\sigma $$\end{document} 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%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$90\%$$\end{document} CI exclusion sensitivity of 2·1025\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\cdot 10^{25}$$\end{document} year with 3 months, and 9·1025\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$9\cdot 10^{25}$$\end{document} year with 5 years of live time. Under the same conditions, the discovery sensitivity after 3 months and 5 years will be 7·1024\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$7\cdot 10^{24}$$\end{document} year and 4·1025\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4\cdot 10^{25}$$\end{document} year, respectively.


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 [1,2]. This decay is currently being investigated with a variety of double beta decaying isotopes. A recent review can be found in Ref. [3]. The cryogenic underground observatory for rare events (CUORE) [4][5][6] is an experiment searching for 0νββ decay in 130 Te. It is located at the Laboratori Nazionali del Gran Sasso of INFN, Italy. In CUORE, 988 TeO 2 crystals with natural 130 Te isotopic abundance and a 750 g average mass are operated simultaneously as source and bolometric detector for the decay. In this way, the 0νββ decay signature is a peak at the Q-value of the reaction (Q ββ , 2527.518 keV for 130 Te [7][8][9]). Bolometric crystals are characterized by an excellent energy resolution (∼0.2% Full Width at Half Maximum, FWHM) and a very low background at Q ββ , which is expected to be at the 10 -2 cts/(keV·kg·yr) level in CUORE [10].
The current best limit on 0νββ decay in 130 Te comes from a combined analysis of the CUORE-0 [11,12] and Cuoricino a e-mail: cuore-spokesperson@lngs.infn.it b Deceased c Presently at: INFN-Laboratori Nazionali di Frascati, Frascati, Rome, 00044, Italy data [13,14]. With a total exposure of 29.6 kg·year, a limit of T 0ν 1/2 > 4.0 · 10 24 year (90% CI) is obtained [15] for the 0νββ decay half life, T 0ν 1/2 . 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νββ 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νββ decay analysis results.
In this work, we dedicate our attention to those factors which could strongly affect the sensitivity, such as the background index (B I ) and the energy resolution at Q ββ . In CUORE, the crystals in the outer part of the array are expected to show a higher B I than those in the middle [10]. Considering this and following the strategy already implemented by the Gerda Collaboration [16,17], we show how the division of the data into subsets with different B I could improve the sensitivity.
The reported results are obtained by means of a Bayesian analysis performed with the Bayesian analysis toolkit (BAT) [18]. 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 [19], this provides the full sensitivity probability distribution and not only its median value.
In Sect. 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 Sect. 4. Finally, we present the results in Sect. 5.

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 ββ , except for the presence of a 60 Co summation peak at 2505.7 keV. Therefore, H 0 is implemented as a flat background distribution plus a Gaussian describing the 60 Co peak. In CUORE-0, this peak was found to be centered at an energy 1.9 ± 0.7 keV higher than that tabulated in literature [15]. This effect, present also in Cuoricino [14], is a feature of all gamma summation peaks. Hence, we will consider the 60 Co 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νββ decay. The latter is modeled as a Gaussian peak at Q ββ . -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νββ decay analysis [15], and is bounded by the possible presence of peaks from 214 Bi at 2447.7 keV and 208 Tl X-ray escape at ∼2585 keV [15]. 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 ∼5 keV FWHM and using a 1 keV bin width, the ±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. -0ν is the parameter describing the 0νββ decay rate for H 1 : θ is the list of nuisance parameters describing the background processes in both H 0 and H 1 ; -is the parameter space for the parameters θ .

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νββ decay analysis, we are interested in the measurement of 0ν for the hypothesis H 1 . The probability distribution for the parameter set ( 0ν , θ ) is: The numerator contains the conditional probability P(E 0ν , θ , H 1 ) of finding the measured data E given the model H 1 for a set of parameters ( 0ν , θ ), times the prior probability π 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 0ν 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 E given the hypothesis H 1 and all possible parameter combinations, P(E|H 1 ). The posterior probability distribution for 0ν is obtained via marginalization, i.e. integrating P 0ν , θ E, H 1 over all nuisance parameters θ : 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 [20]. 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: where λ di is the expectation value for the bin di. The best-fit is defined as the set of parameter values ( 0ν , θ ) for which the likelihood is at its global maximum. In the practical case, we perform the maximization on the log-likelihood 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 λ di . As mentioned above, we parametrize H 0 with a flat distribution over the considered energy range, i.e.
[2470; 2570] keV: plus a Gaussian distribution for the 60 Co peak: 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: where E min di and E max di are the left and right margins of the energy bin di, respectively. Considering bins of size δ E di and expressing M bkg di as function of the background index B I d , of the total mass m d and of the measurement live time t d , we obtain: Similarly, the expectation value for the 60 Co distribution on the bin di is: where M Co d is the total number of 60 Co events for the subset d and can be redefined as function of the 60 Co event rate, R Co d : The total expectation value λ di for H 0 is then: In the case of H 1 an additional expectation value for 0νββ decay is required: The number of 0νββ decay events in the subset d is: where N A is the Avogadro number, m a and f 130 are the molar mass and the isotopic abundance of 130 Te and ε tot is the total efficiency, i.e. the product of the containment efficiency ε MC (obtained with MC simulations) and the instrumental efficiency ε instr .

Exclusion sensitivity
We compute the exclusion sensitivity by means of the 90% CI limit. This is defined as the value of T 0ν 1/2 corresponding to the 90% quantile of the posterior 0ν distribution: An example of posterior probability for 0ν and the relative 90% CI limit is shown in Fig. 1, top. Flat prior distributions are used for all parameters, as described in Sect. 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 0ν 1/2 is to be interpreted as the value above which, given the current knowledge, the true value of T 0ν 1/2 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  Top Marginalized probability distribution for 0ν relative to one toy-MC experiment. The 90% CI limit on T 0ν 1/2 is indicated. Bottom distribution of the 90% CI limits on T 0ν 1/2 for 10 5 toy-MC experiments. The red vertical line corresponds to the median sensitivity, while the three colors depict the 1, 2 and 3 σ quantiles of the distribution. Both plots are obtained with 1 year live time, a 10 -2 cts/(keV·kg·yr) B I and 5 keV energy resolution T 0ν 1/2 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 background-only model, H 0 . We then fit spectra with the background-plussignal model, H 1 , and obtain the T 0ν 1/2 (90% C I ) distribution (Fig. 1, bottom). Its medianT 0ν 1/2 (90% C I ) is referred as the median sensitivity. For a real experiment, the experimental T 0ν 1/2 limit is expected to be above/belowT 0ν 1/2 (90% C I ) with 50% probability. Alternatively, one can consider the mode of the distribution, which corresponds to the most probable T 0ν The exact procedure for the computation of the exclusion sensitivity is the following: ; -we extract the 90% CI limit on T 0ν 1/2 ; -we repeat the algorithm for N toy-MC experiments, and build the distribution of T 0ν 1/2 (90% C I ).

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-plussignal models. A method for the calculation of the Bayesian discovery sensitivity was introduced in Ref. [21]. 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: The application of the Bayes theorem to the models H 0 and H 1 yields: In this case, the numerator contains the probability of measuring the data E given the model H : 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 E given either the model H 0 or H 1 : At this point we need to define a criterion for claiming the discovery of new physics. Our choice is to quote the 3 σ (median) discovery sensitivity, i.e. the value of T 0ν 1/2 for which the posterior probability of the background-only model H 0 given the data is smaller than 0.0027 in 50% of the possible experiments. In other words: 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 0ν 1/2 ; -we fit the spectrum with both H 0 and H 1 ; -we compute P(H 0 |E); -we repeat the procedure for N toy-MC spectra using the same T 0ν 1/2 ; -we repeat the routine with different values of T 0ν 1/2 until the condition of Eq. 20 is satisfied. The iteration is implemented using the bisection method until a 5 · 10 -5 precision is obtained on the median P(H 0 |E).

Experimental parameters
The fit parameters of the H 1 model are B I , R Co and 0ν , while only the first two are present for H 0 . If the data are divided in subsets, different B I and R Co fit parameter are considered for each subset. On the contrary, the inverse 0νββ 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 [22][23][24][25][26][27][28][29]. 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 back- ground model [30]. We developed a full MC simulation of CUORE [10], 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 B I at Q ββ 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 α particles emitted by U and Th surface contaminations of the copper structure holding the crystals. Additionally, we consider a 60 Co contribution normalized to the 90% C.L. limit from the screening measurement. In this sense, the effect of a 60 Co background on the CUORE sensitivity is to be held as an upper limit. Given the 60 Co 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 60 Co half life (5.27 year), and set the start of data taking to January 2017.
The parameter values used for the production of the toy-MC are reported in Table 1. The quoted uncertainty on the B I comes from the CUORE MC simulations [10]. We produce the toy-MC spectra using the best-fit value of the B I . In a second time, we repeat the analysis after increasing and decreasing the B I 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 B I . Namely, being the background at Q ββ dominated by surface α contamination of the copper structure, the crystals facing a larger copper surface are expected to have a larger B I . 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 B I , as reported in Table 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 0ν 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 0ν , would yield a stronger limit. A dif-ferent prior choice based on the real characteristic of the experimental spectra might be more appropriate for B I 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 B I and 60 Co rate uncertainties obtained by the background budget as the σ of a Gaussian prior. No significant difference was found on the sensitivity band because the Poisson fluctuations of the generated number of background and 60 Co events are dominant for the extraction of the 0ν posterior probability distribution. Table 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νββ decay and 60 Co Gaussian peaks. We perform the analysis with a FWHM of 5 and 10 keV, corresponding to a σ of 2.12 and 4.25 keV, respectively. Regarding the efficiency, while in the toy-MC production the B I and R Co are multiplied by the instrumental efficiency, 1 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% [15]. We evaluate the exclusion and discovery sensitivities for different live times, with t ranging from 0.1 to 5 year and using logarithmically increasing values:

Fit procedure
We perform the analysis with the software BAT v1.1.0-DEV [21], which internally uses CUBA [31] v4.2 for the integration of multi-dimensional probabilities and the Metropolis-Hastings algorithm [32] 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, B I 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.

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 B I values and all live times, the FWHM of 5 keV yields a ∼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 Table 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 ±3σ region around Q ββ . The results in Table 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 [19] 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 [33], hence the limit on T 0ν 1/2 is systematically weaker. 2 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 2 See the discussion of the pulls for a more detailed explanation. FWHM the difference is ∼6% after 3 months and 4% after 5 years. One remark has to be done concerning the values reported in [19]: there we gave a 90% C.I. exclusion sensitivity of 9.3 · 10 25 year with 5 year of live time. This is ∼5% higher than the result presented here and is explained by the use of a different total efficiency: 87.4% in [19] and 81.3% in this work.
We then extract the exclusion sensitivity after dividing the crystals in 4 subsets, as described in Sect. 3. The median exclusion sensitivity values after 3 months and 5 years of data collection with one and 4 subsets are reported in Table 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 B I is not worth the effort. This conclusion is not always true, and strongly relies on the exposure and B I of the considered subsets. As an example, we repeated a toy analysis assuming a B I of 10 -2 cts/(keV·kg·yr), and with two subsets of equal exposure and B I 0.5 · 10 -2 cts/(keV·kg·yr) and 1.5 · 10 -2 cts/(keV·kg·yr), respectively. In this case, the division of the data in to two subsets yields a ∼10% improvement after 5 year 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 ∼3 keV to ∼20 keV FWHM [34]. 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νββ 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 [15]. Assuming an average frequency of one calibration per month, the total number of energy spectra would be ∼6·10 4 . Assuming a different but stationary B I for each crystal, and using the same 60 Co rate for all crystals, the fit model would have ∼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  They are computed for each live time value separately as described in Fig. 1, bottom. We also show the sensitivity computed as in [19] in dark green. The horizontal green line at 4 · 10 24 year corresponds to the limit obtained with CUORE-0 and Cuoricino [15]  use of different algorithms, e.g. nested sampling [35,36], 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 60 Co peak. In the first scenario we set the B I to zero, and do the same for the 60 Co 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 60 Co event leaks in the ±3σ region around Q ββ even after 5 year of measurement. As a consequence, the 90% CI limits are distributed on a very narrow band, and the median sensitivity reaches 1.2 · 10 27 year after 5 year of data collection. On the contrary, if we assume a 10 keV FWHM, some 60 Co events fall in the 0νββ 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 B I at 1.02·10 -2 cts/(keV·kg·yr), but set the 60 Co rate to zero. In both cases, the difference with respect to the standard scenario is below 1%. We can conclude that the 60 Co peak with an initial rate of 0.428 cts/(kg·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 λ bkg di and λ 0ν 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 60 Co 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 60 Co events, respectively, and with M bkg and M Co the corresponding number of reconstructed events, the pulls are defined as: where σ 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 ± 0.4 and r Co = 0.1 ± 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 (B I , R Co and 0ν ) 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 ∼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%).

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νββ 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 σ 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 σ discovery sensitivity superior to the limit obtained from the combined analysis of Cuore-0 and Cuoricino data [15] after less than one month of operation, and reaches 3.7 · 10 25 year with 5 year 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.

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 B I 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 B I 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νββ decay and 60 Co peaks can be used, and the systematics of the energy reconstruction can be included.