Combined Constraints on Majorana Masses from Neutrinoless Double Beta Decay Experiments

Combined bounds on the Majorana neutrino mass for light and heavy neutrino exchange mechanisms are derived from current neutrinoless double beta decay (0{\nu}\b{eta}\b{eta}) search results for a variety of nuclear matrix element (NME) models. The approach requires self-consistency of a given model to predict NMEs across different isotopes. The derived bounds are notably stronger than those from any single experiment and show less model-to-model variation, highlighting the advantages of using multiple isotopes in such searches. Projections indicate that the combination of near-term experiments should be able to probe well into the inverted mass hierarchy region. A method to visually represent 0{\nu}\b{eta}\b{eta} experimental results is also suggested to more transparently compare across different isotopes and explicitly track model dependencies.

The origin of neutrino mass is one of the central puzzles in particle physics today. It is intimately connected to the question of whether neutrinos are Dirac or Majorana particles, with fundamental implications for both particle physics and cosmology. Owing to the small scale of neutrino masses, the search for neutrinoless double beta decay (0νββ) is the only known experimental approach that can be practically used to address this question.
However, the translation of experimental half-life bounds from different 0νββ isotopes into Majorana mass limits is complicated by the fact that the nuclear matrix elements (NMEs) for the transition currently have a high degree of uncertainty. There are a variety of different theoretical approaches and parameter choices (e.g. values of coupling constants and correlation functions) that result in predictions varying by a factor of two or more, often with different relative values between isotopes [1]. One particularly important parameter choice relates to the question of a potential "quenching" of the effective axial vector coupling constant, g A , based on discrepancies observed in the rates of single beta decays relative to calculations [2]. Recent ab initio calculations of some beta decay isotopes suggest that the equivalent suppression of the NMEs naturally arises as the result of meson currents and other higher order nucleon interactions [3]. Such interactions are expected to play less of a role for the higher momentum transfers related to 0νββ, but it is widely believed that some suppression will result, though the details remain uncertain. This situation significantly complicates the extraction of robust bounds on Majorana neutrino masses as well as the comparison and combination of experimental results between different isotopes.
The approach taken here is to explicitly treat each NME formalism and set of parameter choices as a separate, self-consistent model that makes linked predictions between different 0νββ isotopes. As such, different experimental results for a given model can then be compared and combined to yield improved mass limits. This approach is similar to that of [4], but here explicitly setting bounds on the model space for light and heavy neutrino exchange mechanisms and making use of likelihood functions based on recent experimental results and NME calculations. In what follows, NME calculations for an unquenched value of g A = 1.27 will be used throughout, however the first order dependence of g A on the halflife (which actually resides in the isotope-specific phase space factor) will be factored out for visibility and to retain model-independence of the isotopic correction.

I. EXPERIMENTAL LIKELIHOOD FUNCTIONS
Four experimental results will be considered in this analysis: CUORE [5], EXO-200 [6], GERDA Phase II [7] and KamLAND-Zen [8,9]. The combination of results will make use of the corresponding likelihood functions.
For CUORE, the likelihood function was directly provided in their paper in terms of a posterior probability with uniform prior as a function of assumed 0νββ decay rate, which has been extracted and parameterized by a polynomial in log(likelihood) over the observable range.
The GERDA Phase II experiment observed zero counts in the region of interest, making the likelihood function trivially a Poisson distribution with n = 0, which is just the exponential e −(µ S +µ B ) , where µ S and µ B are the expected mean number of signal and background events, respectively. However, for relative likelihoods, µ B can be ignored, leaving just e −µ S . EXO-200 results are summarized in their 2019 paper, which gives the numbers of events observed within a ±2σ window for phases 1 and 2, along with the expected background levels and their uncertainties. A likelihood function was therefore constructed based on a Poisson distribution for each phase, taking into account slightly different 0νββ detection efficiencies, and convolving these with a Gaussian for the background uncertainties. A com- bined likelihood was then constructed for a given signal hypothesis, accounting for the exposures of each dataset.
The derived likelihood function for KamLAND-Zen was based on their published observations for the numbers of events seen in each 50 keV energy bin near the signal region compared with their best fit background model. The potential signal fraction in each bin was calculated from a normal distribution centered on the 136 Xe endpoint, with a width commensurate with their stated energy resolution. For simplicity, a Poisson probability distribution was used for each bin assuming that the background model was perfectly constrained by data outside the region of interest. This errs on the side of providing slightly tighter constraints. Separate likelihood constructions were made for Phase I [8], Phase II Period 1 and Phase II Period 2 [9]. Once more, a combined likelihood was then constructed for a given signal hypothesis, taking into account the exposures of each dataset. Figure 1 shows the resulting values of −2 log L R versus the assumed 0νββ signal rate, where L R is the likelihood ratio with respect to the maximum for each experiment.

II. CONSTRUCTION OF HALF-LIFE BOUNDS
To specifically constrain the phase space of possible models, 90% CI Bayesian bounds with a prior that is uniform in the (positive) counting rate will be used here.
For counting experiments searching for rare events, this tends to provide conservative limits that are more robust to statistical background fluctuations than many frequentist constructions while, conveniently, also yielding a relatively close correspondence with numerical values characteristic of frequentist coverage [10]. For 0νββ, such a choice corresponds to a prior that is uniform in m 2 ββ for the light neutrino exchange (LNE) mechanism.
This, indeed, can be seen to yield more conservative upper bounds than, for example, a prior that is uniform in either m ββ or log m ββ . For the heavy neutrino exchange (HNE) mechanism, where one sets lower bounds on the mass scale of the heavy neutrino, a prior uniform in rate corresponds to a prior uniform in M −2 ββ , which is, again, a more conservative prior choice for a lower bound than one that is uniform in either M ββ or log M ββ .
90% CI upper bounds on the 0νββ decay rate, Γ 0ν , have therefore been derived by integrating the posterior probability (which, for a uniform prior, is equivalent to the normalized L R values inferred from Fig. 1) as a function of the event rate from zero until the point where 90% of the distribution is retained. The 90% CI upper bound on the half-life is then log(2)/Γ 0ν 90% . The Bayesian bounds thus produced exactly match those published by CUORE and GERDA. The derived half-life bounds for EXO-200 are slightly more restrictive in appearance than their published frequentist limit (4.3 × 10 25 yr versus 3.5 × 10 25 yr), which might be partly due to a more detailed treatment within their analysis window. On the other hand, the numerical value of the KamLAND-Zen bound is roughly a factor of ∼2 lower than their published frequentist limit (4.9 × 10 25 yr versus 1.07 × 10 26 yr). In both cases, the Bayesian bounds are closer to the expected median experimental sensitivities.

III. MODEL-BY-MODEL CONSTRAINTS
The 0νββ half-life, T 0ν 1/2 , can be related to the effective Majorana neutrino masses (as defined by the PMNS mixing matrix U ) according to [11]: where M 0ν L is the NME for the LNE transition; M 0ν H is the NME for the HNE transition; m ββ is the effective Majorana light neutrino mass ≡ l U 2 el m l ; M ββ is the effective heavy neutrino mass ≡ l U 2 el /M l −1 ; and the factorized form of the isotopic phase space factor, G 0ν = g 4 A G 0ν , is used to separate the dependence on the axial-vector coupling, g A . For simplicity, it will be assumed that we are in a region where one of these two terms in dominant, allowing bounds to be separated placed on LNE and HNE mechanisms.
For a given value of NME, the likelihoods of Fig. 1 can therefore be translated into functions of m ββ or M ββ . The likelihoods can then be combined and used to yield 90% CI bounds on the mass scales for LNE and HNE, as previously described, by integrating the posterior probability as a function of m 2 ββ or M −2 ββ , respectively.  list typifies the span of variation between model predictions. Figure 2 shows examples of the interplay between the experimental likelihoods for four of those NME models in the case of LNE. Derived upper bounds are indicated by filled circles both for the combined likelihood as well as for those of individual experiments. Table II summarizes the derived combined constraints on both LNE and HNE for each NME model.

IV. VISUAL REPRESENTATION FOR TRANSPARENT COMPARISONS
In order to make the nature of the model dependencies as transparent as possible when comparing experimental sensitivities, and to facilitate the relevant comparison across different isotopes, Figures 3 and 4 (for LNE and HNE, respectively) plot the product of T 0ν 1/2 × G 0ν Comb. 90% CI Comb. 90% CI NME Model Upper Bound Lower Bound (gA=1.27) on m ββ (meV) on M ββ (GeV) IBM2 [12] 70.5 7.6 × 10 7 CDFT [13] 62 1.1 × 10 8 QRPA-FFS [14] 144 7.9 × 10 7 QRPA-JY [15] 82 1.6 × 10 8 QRPA-TU [16,17] 92.5 1.25 × 10 8 ISM-TK [18] 107.5 7.6 × 10 7 QRPA-NC [19] 100 ISM-INFN [20] 120 CGM [21] 57 for the different experimental results along the x-axis. This unitless "TG sensitivity" provides a useful modelindependent, first order isotope correction to indicate the rough physics reach of a given observation before the uncertain NME values are taken into account. A given experimental bound results in a vertical line at the corresponding TG value, with the half-life bound also given at the base of each line. The two dashed vertical lines are here used to indicate the difference between the Bayesian bounds and published frequentist values for EXO-200 and KamLAND-Zen. Different matrix element models are represented by the different symbols indicated, each allowing a model-dependent translation to an effective neutrino mass on the y-axis. Different isotopes are color coded and the diagonal dotted black lines indicate how the mass sensitivity scales with TG sensitivity. The region approximately corresponding to the inverted neutrino mass hierarchy in Fig. 3 is indicated by the horizontal blue band (which also includes part of the normal mass hierarchy for degenerate neutrinos). The combined bounds derived in this work are indicated on the right side of each plot. These bounds are notably more restrictive than those from any single experiment.

V. PROJECTED COMBINED BOUNDS FOR NEAR-TERM EXPERIMENTS
A similar approach has been taken to project the physics reach for the combination of running and near-term experiments: CUORE, KamLAND-Zen 800, LEGEND-200 and SNO+ I. Five years of live time with the full experiment was assumed for each.
For CUORE and LEGEND-200, a single-bin Poisson model spanning ±1.5σ about the endpoint was used for the likelihood, assuming similar signal efficiencies as current versions of the two technologies. In the case of CUORE, a background index of 1.38 × 10 −2 counts/keVkg(T eO 2 )-yr was used [5] and, for LEGEND-200, a value of 2 × 10 −4 counts/keV-kg( 76 Ge)-yr was assumed [22].
For KamLAND-Zen and SNO+, a Poisson-based likelihood with multiple energy bins in the region of interest was used. Background contributions from 8 B solar neutrinos and 2νββ "spill-over" due to energy resolution were directly calculated. In the case of SNO+, additional internal and external backgrounds were added to produce a model consistent with reference [23]. For KamLAND-Zen, it was assumed that U/Th backgrounds could be reduced to negligible levels by improvements to analysis and electronics and that the fiducial volume could be extended to a radius of 1.65m so as to achieve their sensitivity goal of T 1/2 > 5 × 10 26 yrs [23].
The resulting 90% CI sensitivities assuming no observed excess for both individual and combined likelihoods (including current measurements) are shown in Fig. 5. Again, projected combined bounds are signifi- cantly more restrictive than those from any single measurement and suggests that the combination of near-term experiments will be able to probe well into the region corresponding to the inverted neutrino mass hierarchy.

VI. CONCLUSION
Combined constraints on Majorana masses have been derived based on the results of recent 0νββ experiments by demanding self-consistency of any given NME model (defined as the combination of formalism and parameter choice) in its predictions across different isotopes. Parameterized or derived likelihood functions for results from CUORE, EXO-200, GERDA and KamLAND-Zen were used to arrive at 90% CI bounds for LNE and HNE mechanisms using a conservative choice of prior. The derived combined constraints are notably more restrictive than those from any one experiment and also have less variation across different NME models, highlighting the complementarity of these different approaches and the advantages of using different isotopes in the search for 0νββ. While uncertainties on NME values remain, projections for near-term experiments suggest that their combination will be able to probe well into the IH region for the NME models considered here. A method to visually display 0νββ results has also been suggested to allow for the transparent comparison of sensitivities and model-dependencies across different isotopes.
The author would like to thank Frank Deppisch and Josh Klein for helpful discussions. This work is supported by the Science and Technology Facilities Council of the United Kingdom.