A multi-isotope $0\nu2\beta$ bolometric experiment

There are valuable arguments to perform neutrinoless double beta ($0\nu2\beta$) decay experiments with several nuclei: the uncertainty of nuclear-matrix-ele\-ment calculations; the possibility to test these calculations by using the ratio of the measured lifetimes; the unpredictability of possible breakthroughs in the detection technique; the difficulty to foresee background in $0\nu2\beta$ decay searches; the limited amount of isotopically enriched materials. We propose therefore approaches to estimate the Majorana neutrino mass by combining experimental data collected with different $0\nu2\beta$ decay candidates. In particular, we apply our methods to a next-generation experiment based on scintillating and Cherenkov-radiation bolometers. Current results indicate that this technology can effectively study up to four different isotopes simultaneously ($^{82}$Se, $^{100}$Mo, $^{116}$Cd and $^{130}$Te), embedded in detectors which share the same concepts and environment. We show that the combined information on the Majorana neutrino mass extracted from a multi-candidate bolometric experiment is competitive with that achievable with a single isotope, once that the cryogenic experimental volume is fixed. The remarkable conceptual and technical advantages of a multi-isotope investigation are discussed. This approach can be naturally applied to the proposed CUPID project, follow-up of the CUORE experiment that is currently taking data in the Gran Sasso underground laboratory.


Introduction
Neutrinoless double beta (0ν2β) decay is considered as one of the most promising ways to investigate the properties of neutrino and weak interactions, to test lepton number conservation, and in general to probe the Standard Model (SM) of particle physics [1,2,3,4,5,6,7]. Despite the almost seventy-year-long history of 0ν2β decay search, this elusive process is still not observed: the most sensitive experiments provide limits on the effective Majorana neutrino mass of the electron neutrino on the level of lim m ν ∼ 0.1 − 1 eV (see the status of the 0ν2β decay search in the reviews [2,3,4,5,6,7,8,9,10,11,12] and the recent experimental results [13,14,15,16,17,18]).
The exploration of the inverted neutrino mass hierarchy ( m ν ∼ 0.02 − 0.05 eV) is a challenge of the experiments in preparation or R&D stage. Such a sensitivity can be reached in searches able to detect an extremely low activity with half-life in the range T 1/2 ∼ 10 26 − 10 27 yr for the most promising candidates. However, taking into account the possible quenching of the axial vector coupling constant g A [19,20,21,22], the experimental program should foresee the possibility to go towards an-order-of-magnitude higher half-life sensitivities, of the order of T 1/2 ∼ 10 27 − 10 28 yr. The accomplishment of such an ambitious plan requires the construction of detectors containing a large number of 2β active nuclei (10 27 −10 28 nuclei, ∼ 10 3 −10 4 moles of isotope of interest), with an as-high-as-possible detection efficiency and an extremely low (ideally zero) radioactive background, and able to distinguish the effect searched for. The last requirement implies, in particular, an as-high-as-possible energy resolution -to identify the effect unambiguously and to avoid the presence of a tail from the SM-allowed two-neutrino double beta (2ν2β) decay in the region of the 0ν2β peak (region of interest, ROI) -and time resolution -to avoid random coincidences of events from 2ν2β decays and other components of the background [23].
The technique of scintillating and Cherenkov-radiation low-temperature bolometers [24,25,26,27,28,29] is one of the most suitable detection method for the next generation 2β experiments, thanks to the high detection efficiency (70%−90%, depending on the crystal composition and volume) and the excellent energy resolution (a few keV). In addition, this technique suppresses the background caused by α radioactive contamination to a very low level thanks to the excellent particle discrimination. Following this approach, the CUPID group of interest [30] proposes to perform a tonne-scale bolometric 0ν2β decay search with a sensitivity high enough to probe the inverted hierarchy of the neutrino mass and even go towards the normal hierarchy region. This experiment would be a follow-up to the CUORE experiment with particle identification for background rejection.
There are both theoretical and experimental arguments to develop 0ν2β decay experiments with several nuclei. First of all, the ambiguity of the theoretical calculations of nuclear matrix elements (NME) of the 0ν2β decay remains significant, not to say that there are no experimental tests of the calculations so far. Therefore, no one can safely suggest the nucleus with the highest decay probability. Furthermore, data on 0ν2β decay rate in a few nuclei will be definitely useful to adjust theoretical calculations of the NME since the ratio of the calculated lifetimes is very sensitive to different models [31]. From the experimental point of view, taking into account the extremely low decay probability and the possibility that spurious effects mimick the signal searched for, an observation of 0ν2β decay in different nuclei will be requested in case of positive evidence in one of the candidates. Indeed, there have already been many false claims of double beta decay in the history of the experimental search for this process [32].
Therefore, there is a strong call for experiments with several nuclei candidates. At the same time one can extract competitive estimations of the Majorana neutrino mass by combining data of experiments with several different isotopes. In this paper we demonstrate that the sensitivity to the Majorana neutrino mass that can be achieved in a cryogenic experiment with several nuclei is similar to the sensitivity of an experiment utilizing only one of the candidates and exploiting the same experimental volume. Of course, the same argument can be applied to other physics parameters related to different mechanisms inducing 0ν2β decay, although deepening this aspect goes beyond the scope of this paper.

Experimental concepts
On the basis of the above considerations, we propose a next generation 0ν2β decay bolometric experiment studying several isotopes simultaneously. To complete the general discussion of the previous Section, we stress that there are also practical, but not less important, reasons to implement such cryogenic 0ν2β experiment: the production of a large amount of a given isotope (at the hundreds of kg scale) looks rather questionable, while the production, e.g., of one quarter of the same quantity for four isotopes looks more realistic, given the existing enrichment techniques and capabilities. We remark that different sets of centrifugues are used for different nuclei in the centrifugation enrichment technique (presently the only one that provides a high enough throughput). Enrichment in different isotopes can therefore be performed in parallel. A similar argument applies to crystal production, which is mandatory for the bolometric approach. Crystals containing different candidates can be grown in parallel, since each type of crystal requires specific technologies and is typically produced in a specialized laboratory or company. A multi-isotope experiment can therefore provide a remarkable reduction of the time required to build the detector, especially in the case, discussed below, in which the single module structure is the same for each type of crystal.
We will consider the candidates 82 Se, 100 Mo, 116 Cd and 130 Te, which can be studied with the scintillating (the first three ones) and Cherenkov-radiation (the last one) bolometric technique with excellent prospects. Although our considerations apply to a generic nextgeneration cryogenic experiment, the most natural implementation of this multi-isotope approach is in the framework of the CUPID proposal. Therefore, in the following we will refer to the four-isotope option as CUPID-4, while we will designate CUPID-1 the conventional one-isotope version preliminarily discussed in Ref. [30].
As for the detector technology, we will stick to the basic configuration of the single module of CUORE [14,15], consisting of a large crystal (125 cm 3 in case of CUORE) coupled to a neutron transmutation doped (NTD) Ge thermistor, working as a temperature sensor, and to a Si-doped heater for detector response stabilization. Each crystal is hold by small PTFE elements inside a copper frame. The read-out electronics operates at room temperature and is based on the well-known Si-JFET technology. CUORE is taking data successfully [15] and there is no reason to change substantially this experimental configuration.
The bolometric light detector, which is not present in CUORE but will be necessary to reject the α background in future bolometric searches, will consist of a thin Ge disk (50 mm diameter, 0.15 mm thickness) with an NTD Ge thermistor as well. It will be operated in the Neganov-Luke mode in order to achieve a high enough sensitivity to reject α background in the TeO 2 crystals -used in the 130 Te section -by exploiting Cherenkov radiation [28,33,34] and the two-neutrino double beta decay random coincidences in the 100 Mo section [35]. The same Si-JFET-based read-out electronics adopted for the main crystal, and used now in CUORE, can be employed. The same type of light detectors can be coupled to all the four crystal types, but the Neganov-Luke option is not strictly necessary in the 82 Se and 116 Cd cases. Of course, other technologies for the light detector may be viable: a vibrant R&D activity is ongoing on this subject [29,30].
We will assume an array of 1300 scintillating and Cherenkov-radiation bolometers (compatible with the experimental volume of the current CUORE cryostat or of a similar dedicated infrastructure) based on cylindrical ultra-radiopure crystals -5 cm diameter and 5 cm height -enriched to 100% in the isotope of interest. In case of the CUPID-1 option, the crystals will consist of only one of the following compounds: zinc selenide (Zn 82 Se), lithium molybdate (Li 2 100 MoO 4 ), cadmium tungstate ( 116 CdWO 4 ), or tellurium oxide ( 130 TeO 2 ). As for the CUPID-4 option, we will assume that the 1300 elements will be shared among crystals of each type. We stress that a remarkable advantage of the bolometric technique for a multi-isotope experiment is that detectors are based exactly on the same technology, are operated in the same set-up and are exposed to the same background sources.
All the four aforementioned compounds have been successfully tested as scintillating or Cherenkov-radiation bolometers, using crystals grown from enriched materials. It was shown as well that the Neganov-Luke light detectors have the requested performance [34]. The energy resolution of the detectors at the energies Q 2β were estimated from existing experimental data concerning ZnSe [36], Li 2 MoO 4 [37], CdWO 4 [38], and TeO 2 [39]. The detection efficiencies for the 0ν2β decay effect were taken from Ref. [27] (ZnSe and TeO 2 ), while for Li 2 MoO 4 and CdWO 4 detectors the efficiencies were estimated by using the GEANT4 simulation package.
3 Approaches to estimate neutrino mass from several 0ν2β experiments A value (limit) on the effective Majorana neutrino mass can be derived from the experimental data collected by an experiment following the CUPID-4 scheme. Attempts to estimate a limit on the effective Majorana neutrino mass combining data of the most sensitive experiments are already reported in the literature [40,41]. We propose here other approaches driven also by the quite similar performance of the cryogenic detectors foreseen in CUPID.

Sum of counts in the region of interest
Let us consider the mass mechanism for light neutrino exchange, which leads to the following expression for the half-life T 0ν 1/2 of the nuclei candidates to the 0ν2β decay: where G 0ν is the phase-space factor, M 0ν is the NME, g A is the axial vector coupling constant, m ν is the Majorana neutrino mass, and m e is the mass of electron. The number of detected 0ν2β events S i in an i-th 0ν2β experiment with N i nuclei, with a detection efficiency ε i , over a live time t i is inversely proportional to the nucleus half-life: The sum of 0ν2β events acquired by several detectors containing different 2β isotopes can be expressed as following: Therefore, the Majorana neutrino mass can be derived from several experiments by using the following formula: S (lim S) can be obtained from the analysis of the experimental data accumulated with several detectors as following (assuming a large enough number of counts): However, the method of sum of counts has a drawback in the case of different detectors performance: if all the S i values are close to 0, the biggest contribution will be given by the measurement with the highest background. Thus, the neutrino mass limit could be defined by the worst experimental case, which is not logical.
Nevertheless, as it will be shown in Sections 4.3 and 4.4, the method is quite effective in case of experiments with similar mass of the isotope of interest, background, efficiency and energy resolution, which is expected to be the case of CUPID-4.
In Section 4.4 we will consider "zero background" CUPID-like experiments. Approximation of data with mainly bins with zero counts is not robust enough. Thus, one can calculate the sum of counts in selected energy intervals to estimate experimental sensitivities to the Majorana neutrino mass by using only the given background and no true signal.

Weighted averages and errors of neutrino mass square
The following method allows us to suppress the impact of measurements with high backgrounds. A value (limit) of neutrino mass can be derived from several experiments with different nuclei by using weighted averages and errors of square of the neutrino mass obtained in different experiments as recommended by the Particle Data Group [42] (assuming once again a large enough number of counts): where In this approach negative values of the square values of the neutrino mass m ν 2 i can appear due to possible negative values of the peaks area. Finally, a combined limit on the neutrino mass can be calculated as the square root of the lim m ν 2 , obtained from the values of m ν 2 and σ m ν 2 . This method reduces the impact of experiments with high background on the combined neutrino mass limit. This approach can be applied to extract a combined limit even from experiments with rather different characteristics (e.g., bolometers, HPGe detectors, Xe-filled scintillators, etc).
Below we will show how the discussed methods can be applied to the case of CUPID-like experiments.

Input data
The various contributions to the background in a large scintillating and Cherenkov-radiation bolometric experiment and in the CUORE experiment are discussed respectively in Ref. [27] and Ref. [43]. According to the discussions in these works, we will assume that the external contributions to the background (external γ radiation, muons, neutrons) can be reduced to below 10 −4 counts/yr/keV/kg by proper material selection and shielding. A simililarly low contribution will be provided by the bulk crystal radioactivity. The dominant source of background will be then related to contamination of close structures, typically the surface of the crystals and of the holders. Presently, α particles generated at the surface of the crystal holders are indeed the dominant background source in CUORE, at the level of 10 −2 counts/yr/keV/kg [15,43]. Switching on particle identification capability by the detection of light emitted by the crystals, the residual level of background is estimated by the simulations based on the CUORE background model and reported in Ref. [44]. The remaining background events are determined by β's and  Table 1. The masses of the detectors and the numbers of nuclei should be divided by 4 to obtain the corresponding values for CUPID-4. We assume 10 yr of data taking in both configurations. We have not tried to optimize the sharing among the different isotopes, as this work intends just to show in general the potential of a multi-isotope bolometric approach. It is likely that, after the conclusion of the current CUPID R&D phase, we will find out that the optimum crystal size and shape as well as the technical performances of the various detectors depend on the compound type: therefore, a sharing based on an equal number of crystals could not be the best choice. CUPID-2 or CUPID-3 options could also result more convenient. The choice of the isotopes can also depend on the future results of CUORE: a larger statistics will lead to a better definition of the background model and this can influence the 0ν2β candidate selection, given the different values of the Q 2β 's.

Sensitivity of CUPID-1 experiments
First we estimate the sensitivity to the 0ν2β decay of 82 Se, 100 Mo, 116 Cd and 130 Te in the CUPID-1 option, assuming the simulated background reported in Ref. [44]. 1 Ten thousands of randomly generated energy 1 Strictly speaking, the background counting rate in Ref. [44] was estimated for TeO 2 crystals. We use the simulation re- spectra were fitted by the maximum-likelihood method 2 to estimate the average errors of the peak area σ(S i ) for each nucleus. Examples of the fits are presented in Fig. 1 (energy spectra with S i ≈ 0 and σ(S i ) equal to the average values for each detector were chosen for the presentation). The values of lim S i were estimated by using the Feldman-Cousins procedure [45] from the σ(S i ) values assuming S i = 0 counts. A slightly higher sensitivity was achieved by using energy spectra with bins in units of energy resolution as it is shown in Fig.  2. 3 This improvement can be explained by an increase of statistics thanks to the binning in units of energy resolution. The results are presented in Table 2. 4 The neutrino mass limits here and below have been calculated by using the matrix elements from Ref. [46], the phase space factors from Ref. [47] and the value of the axial vector coupling constant g A = 1.27. We would like to stress that the choice of the matrix elements from Ref. [46] does not mean that we give preference to these calculations, but it is used only to illustrate our approach. In Section 4.5 we will give a picture of the CUPID-1 and CUPID-4 sensitivity to neutrino mass using other calculations too. In the right column of Fig. 1 positive signals for the effective neutrino mass m ν = 0.05 eV are shown (we have used the same NME calculations to sults also for the other compounds (taking into account of course the different Q 2β 's) assuming that the background level does not depend significantly on the detector material, as expected if the contamination is mainly located in the holders. 2 We have applied the likelihood fit taking into account that approximation by the least squares method is not robust enough in the case of low statistics data, since it ignores zero entry bins. estimate the numbers of counts as for the limit estimations above).

Sensitivity of CUPID-4 and combined limits on the neutrino mass
The sensitivity of CUPID-4 was estimated by using both the sum of counts in the ROI and the weighted averages methods.   Table 1 for details). The figures at the left column show the "no effect observed" case, while those at the right column illustrate the effects for a neutrino mass of 0.05 eV calculated with the matrix elements from Ref. [46], the phase space factors from Ref. [47] and the value of the axial vector coupling constant g A = 1.27. Values of the excluded lim S and peaks areas S are in counts.  Counts/keV Fig. 2 The same spectra as in Fig. 1 are reported, but binning is units of energy resolution (1 bin = FWHM(keV)/5). This approach leads to a slightly better sensitivity, as discussed in the text.
Values of S and σ(S) can be estimated by using formula (5) and by a fit of the sum spectrum of several detectors. It should be noted that one can sum spectra measured with different nuclei (and therefore with different Q 2β values). For this purpose the spectra can be binned equally (e.g., 1 keV per channel, as in Fig. 1) or in units of energy resolution (as in Fig. 2). They can then be shifted in energy so that the decay energies of all the nuclei are in the same channel. Binning in units of energy resolution improves statistics and looks a preferable approach in case of different energy resolutions. Below we will use energy spectra binned in units of energy resolution.
The generated energy spectra of the four sets of detectors and the sum spectrum are shown in Fig. 3. The results of the fits are presented in Fig. 3 and Table 3.
An estimation by using formula (5) gives lim m ν = 0.020 eV. Then the sum spectrum was fitted by the sum of four Gaussian functions with energy resolution FWHM=5 channels. This approach gives almost the same result: lim m ν = 0.021 eV. The excluded peaks for the four detectors and for the sum spectrum are shown in Fig. 3.
The best limit lim m ν = 0.018 eV was obtained with the second method by using the values of weighted errors of squares of neutrino mass σ m ν 2 listed in Table 3.
Energy spectra generated for a Majorana neutrino mass of 0.05 eV and their fits are given in the right column of Fig. 3. One can see that an effect for m ν = 0.05 eV will be reliably observed both by CUPID-1 and CUPID-4 with the background levels extracted from Ref. [44]. However, only the simultaneous observation in more than one isotope -with excesses of counts in the correct energy positions -would give overwhelming evidence for the occurrence of 0ν2β decay, providing a cross check of the results in the moment itself of the discovery.

Sensitivity of "zero-background" experiments
We will now assume that the background will be suppressed down to the level of 4×10 −6 counts/yr/keV/kg. This value is very challenging and we use it here just to simulate the case of a "zero background experiment" both in the CUPID-1 and CUPID-4 options. Generated energy spectra for Zn 82 Se, Li 2 100 MoO 4 , 116 CdWO 4 and 130 TeO 2 detectors for the CUPID-1 option are presented in Fig. 4, while the energy spectra for the CUPID-4 configuration are shown in Fig. 5. A fit of such statistically poor data is questionable. Thus, we have estimated the sensitivity of this thought experiments in accordance with the procedure reported in Ref. [45] for an expected background and no true signal (see Table XII in Ref. [45]). The estimations of Table 3 Sensitivity of the CUPID-4 experiment with Zn 82 Se, Li 2 100 MoO 4 , 116 CdWO 4 and 130 TeO 2 detectors assuming the level of background simulated in Ref. [44] (denoted as I BG , counts/yr/keV/kg). The errors σ m ν 2 and the neutrino mass limits were calculated by using the matrix elements from Ref. [46], the phase space factors from Ref. [47] and the value of the axial vector coupling constant g A = 1.27.
Isotope Approach to estimate  [44]. The figures at the left column show the "no effect observed" case, while those at the right column illustrate effects for a neutrino mass of 0.05 eV calculated with the matrix elements from Ref. [46], the phase space factors from Ref. [47] and the value of the axial vector coupling constant g A = 1.27. Values of excluded lim S and peaks areas S are in counts.
the experimental sensitivity for the CUPID-1 option are reported in Table 4. The expected background N BG for each detector and for the sum spectrum was evaluated in ≈ ±2σ energy intervals (the intervals are shown in Fig. 4 and 5). The estimations of the CUPID-4 experimental sensitivity are presented in Table 5.  Table 4) are shown. The figures at the right column illustrate effects for a neutrino mass of 0.02 eV calculated with the matrix elements from Ref. [46], the phase space factors from Ref. [47] and the value of the axial vector coupling constant g A = 1.27. Values of excluded lim S and peaks areas S are in counts.
The energy spectra for the Majorana neutrino mass m ν = 0.02 eV are presented in the right columns of Figs. 4 and 5. The effect can be detected in both CUPID configurations. show the "no effect observed" case. The energy intervals where the numbers of background counts N BG were estimated (the N BG values are given in Table 5) are shown. The figures at the right column illustrate effects for a neutrino mass of 0.02 eV calculated with the matrix elements from Ref. [46], the phase space factors from Ref. [47] and the value of the axial vector coupling constant g A = 1.27. Values of excluded lim S and peaks areas S are in counts.

Dependence on nuclear matrix elements calculations
The dependence of the CUPID thought experiment sensitivity to the effective Majorana neutrino mass on the NME calculations is shown in Fig. 6. The experimental conditions are the same as those reported in Tables 4 and 5. The values of lim m ν were estimated by using the phase space factors from Ref. [47], the axial vector coupling constant g A = 1.27, and the NME discussed in a recent review [49] (however, we have considered only the works where all the "CUPID nuclei" -82 Se, 100 Mo, 116 Cd and 130 Te -were calculated): the different versions of the quasiparticle random-phase approximation from the Tübingen (QRPA Tu) [50] and Jyväskylä (QRPA Jy) [51]; the microscopic interacting boson model IBM-2 [46]; and energy density functional theory (EDF), non-relativistic (NR-EDF) [52] and relativistic (R-EDF) [53].

Conclusions
Two approaches to estimate the Majorana neutrino mass from 0ν2β decay experiments with several nuclei are proposed. The former method uses a sum of counts in the region of interest, while the latter one estimates the Majorana neutrino mass from weighted averages and errors of neutrino mass squares. The former approach is especially effective in case of similar detector performance and can be used even in the "zero background" conditions with a very poor statistics, while the latter approach allows one to estimate the neutrino mass also from experiments with completely different detectors.
By applying these methods, we have shown that the sensitivity to the neutrino mass that can be achieved in a single bolometric experiment with several nuclei (e.g., with 82 Se, 100   counts/yr/keV/kg. The values of lim m ν were estimated by using the phase space factors from Ref. [47], the value of the axial vector coupling constant g A = 1.27, and the NME calculated in the framework of different nuclear models: QRPA Tu [50]; QRPA Jy [51]; IBM-2 [46]; NR-EDF [52]; and R-EDF [53]. based on only one of these nuclei and contained in the same cryostat volume. Actually, Tables 2 and 4, as well as Fig. 6, show that the pure 130 Te option gives better results. We have to consider however that predictions of background are less safe for this nucleus, since this is the only one for which the signal is below the end point of the bulk of the natural γ radioactivity (2615 keV). A possible way to mitigate the effects of a residual background component from external γs is to place the TeO 2 crystals in the central part of the array so that the other crystals act as an active shield for the 130 Te section.
It should be stressed that a similar data analysis can be applied not only to the Majorana neutrino mass assessment in the framework of the mass mechanism involving light neutrino exchange, but also to the estimation of physics parameters related to other possible extensions of the SM that predict 0ν2β decay.
Since a multi-isotope bolometric experiment is not less sensitive than a single-isotope one to the physics parameters that rule 0ν2β decay (once that background and NME uncertainties are taken into account), we believe that it represents by far a better option, for a number of fundamental and practical reasons we have extensively discussed in the previous sections. This multi-isotope approach can be naturally implemented in the framework of CUPID, which can be based in principle on very similar detectors containing different candidates and exploiting well-established technologies. We emphasize that a competitive single experiment combining results from several candidates with similar sensitivities represents a new original idea in the 0ν2β decay search strategy, which currently can be effectively implemented only in the bolometric approach.
Finally, it should be noted that a multi-isotope experiment could allow for precise measurements of the 2ν2β decay half-lives, which are important nuclear data also for the development of a model aiming at accurate calculations of the neutrinoless mode of the transition.

Acknowledgements
The authors thank Dr. V.V. Kobychev for the Monte-Carlo simulations of 0ν2β decay detection efficiencies of Li 2 MoO 4 and CdWO 4 detectors. The researches were supported in part by the joint scientific project "Development of Cd-based scintillating bolometers to search for neutrinoless double-beta decay of 116 Cd" in the framework of the PICS (Program of International Cooperation in Science) of CNRS in years 2016-2018. F.A.D. gratefully acknowledges support from the Scientific High-Level Scholarship Commission of the French Embassy in Kyiv (Ukraine) and the "Jean d'Alembert" Grants program of the University of Paris-Saclay. F.A.D. and V.I.T. were supported in part by the IDEATE International Associated Laboratory (LIA) and by the project "Investigation of neutrino and weak interaction in double beta decay of 100 Mo" in the framework of the Programme "Dnipro" based on Ukraine-France Agreement on Cultural, Scientific and Technological Cooperation.