Majorana neutrino mass constraints in the landscape of nuclear matrix elements

We discuss up-to-date constraints on the Majorana neutrino mass $m_{\beta\beta}$ from neutrinoless double beta decay ($0\nu\beta\beta$) searches in experiments using different isotopes: KamLAND-Zen and EXO ($^{136}$Xe), GERDA and MAJORANA ($^{76}$Ge) and CUORE ($^{130}$Te). Best fits and upper bounds on $m_{\beta\beta}$ are explored in the general landscape of nuclear matrix elements (NME), as well as for specific NME values obtained in representative nuclear models. By approximating the likelihood of $0\nu\beta\beta$ signals through quadratic forms, the analysis of separate and combined isotope data becomes exceedingly simple, and allows to clarify various aspects of multi-isotope data combinations. In particular, we analyze the relative impact of different data in setting upper bounds on $m_{\beta\beta}$, as well as the conditions leading to nonzero $m_{\beta\beta}$ at best fit, for variable values of the NMEs. Detailed results on $m_{\beta\beta}$ from various combinations of data are reported in graphical and numerical form. Implications for future $0\nu\beta\beta$ data analyses and NME calculations are briefly discussed.

We discuss up-to-date constraints on the Majorana neutrino mass m ββ from neutrinoless double beta decay (0νββ) searches in experiments using different isotopes: KamLAND-Zen and EXO ( 136 Xe), GERDA and MAJORANA ( 76 Ge) and CUORE ( 130 Te). Best fits and upper bounds on m ββ are explored in the general landscape of nuclear matrix elements (NME), as well as for specific NME values obtained in representative nuclear models. By approximating the likelihood of 0νββ signals through quadratic forms, the analysis of separate and combined isotope data becomes exceedingly simple, and allows to clarify various aspects of multi-isotope data combinations. In particular, we analyze the relative impact of different data in setting upper bounds on m ββ , as well as the conditions leading to nonzero m ββ at best fit, for variable values of the NMEs. Detailed results on m ββ from various combinations of data are reported in graphical and numerical form. Implications for future 0νββ data analyses and NME calculations are briefly discussed.

I. INTRODUCTION
The process of neutrinoless double beta decay (0νββ), expected to occur for some candidate isotopes (Z, A) if neutrinos are Majorana fermions, may be interpreted as a miniature event of leptonic matter creation or "Little Bang," whose discovery would have profound implications for particle and nuclear physics and for cosmology [1].
In the standard three-neutrino paradigm [2], the process would be mediated by three Majorana neutrino mass states ν i (i = 1, 2, 3) mixed with the three known flavor states ν α (α = e, µ, τ ) via a unitary mixing matrix U αi , parametrized in terms of three mixing angles (θ 12 θ 13 , θ 23 ), one (Dirac) phase δ, and two (Majorana) phases ϕ 1,2 . The relevant particle physics parameter is the effective Majorana neutrino mass m ββ , defined as and related to the observable 0νββ decay half-life T i in each isotope i = (Z, A) via where G i is the phase space, and M i is the nuclear matrix element (NME) for the decay. It is useful to contrast the Majorana ν mass m ββ with the sum of neutrino masses that, being a source of gravity, can produce observable cosmological effects [2]. Figure 1 shows the regions allowed in the (Σ, m ββ ) plane at the 2σ level (∆χ 2 = 4) by a global analysis of neutrino oscillation data [3], for masses m i either in normal ordering (NO, m 1,2 < m 3 ) or in inverted ordering (IO, m 3 < m 1,2 ). For a given value of Σ, the vertical spread of m ββ is mostly due to the unknown relative phases of the U ei matrix elements in Eq. (2). Current cosmological data provide typical upper bounds on Σ at the level of O(100) meV, that are more easily accommodated in NO than in IO [3]; see also the overview of recent constraints on Σ and their impact on ν mass ordering in [4]. Several 0νββ decay searches are also exploring the O(100) meV range for m ββ [1,2]; in particular, the latest constraints from KamLAND-Zen [5] ( 136 Xe) can plunge into the region m ββ ∼ few × 10 meV for favorable values of the NME. 1 Other very sensitive 0νββ searches, all probing half-lives T i > 10 25 y at 90% C.L., have been performed by the experiments EXO [6] ( 136 Xe), GERDA [7] and MAJORANA [8] ( 76 Ge), and CUORE [9] ( 130 Te). 1 The range m ββ 16-49 meV, spanning the leftmost edge of the IO region in Fig. 1, is often (but improperly) dubbed in 0νββ jargon as "IO region," despite being compatible with both IO and NO (as well as the quasi-degenerate region at higher m ββ ). The misnomer may originate from plots of m ββ versus the lightest ν mass m l , where two elongated stripes for IO and NO appear in log scale as m l → 0, see e.g. [1,5]. However, the asymptotic separation of such stripes has no physical relevance, since m l is not directly measurable and cannot be resolved with an accuracy better than an observable such as Σ. Projecting away m l (as in Fig. 1) makes the point clear. Building upon previous work [3], we discuss in detail how to combine current (Xe, Ge, Te) data for given NME values. 2 The approach allows a do-it-yourself 0νββ global analysis in terms of χ 2 functions with (up to) quadratic dependence on the signal strength 1/T i , which is a good approximation to recent results [3,10,11]. In particular, we describe how to derive m ββ constraints at a given confidence level, using both separate and combined (Xe, Ge, Te) data, for generic values of the nuclear matrix elements (the "NME landscape"), as well as for representative NME values from different nuclear models. Our approach clarifies interesting aspects of the 0νββ data analysis, such as the relative importance of each isotope in determining (non)zero best fits and upper bounds for m ββ .
The paper is structured as follows: In Section II we describe the ingredients of our analysis in terms of notation, parametrization of experimental results for (Xe, Ge, Te), and associated NME's. In Section III we discuss the main results of the analysis in terms of m ββ constraints, by considering two qualitatively different situations: (1) cases where, a priori, m ββ = 0 is preferred, and (2) more general cases where the best fit may be at m ββ > 0. Upper bounds on m ββ are explored both graphically and numerically in the NME landscape, by using separate and combined (Xe, Ge, Te) data. In Sec. IV we summarize our results and comment on further applications and perspectives.

II. INGREDIENTS OF THE ANALYSIS
In this Section we introduce the notation, the experimental results and their parametrization, the landscape of NME and the phase space related to the three isotopes Xe, Ge and Te.

A. Notation and units
Following [3], we introduce the inverse half-life that represents, up to a constant factor, the observable decay rate or signal strength in each i = (Z, A) isotope. Equation (3) reads then To keep the notation compact, we absorb in G i terms as 1/m 2 e and g 4 A (where g A = 1.276 [12] is the bare value of the axial-vector coupling), that are factorized out in other conventions. In particular, we can make contact with the notation of [1], where 1/T = G 01 g 4 A M 2 0ν m 2 ββ /m 2 e , by identifying G = G 01 g 4 A /m 2 e and M = M 0ν for each isotope i. We also follow [1] by taking the M i as positive real numbers, referred to the bare value of g A (unless otherwise noticed). Qualitative effects of the so-called quenching of g A in nuclear matter [13] are separately commented below.

B. Experimental inputs and parametrizations
In principle, the 0νββ data analysis would be straightforward, if likelihood profiles were provided for the signal strength S i (or for T i ) in each experiment, e.g., in terms of a function χ 2 i (S i ). Barring error correlations among independent experiments, one should sum up the χ 2 i functions, express the S i in terms of m ββ via Eq. (6) for a given set of M i , and map the resulting best fits and allowed regions for m ββ . In practice, experimental papers often focus on a single point of the likelihood profile (e.g., the T i bound at 90% C.L.), whereas its shape has to be derived from supplementary information.
A useful empirical fact, first noted in [10] and further elaborated in [3], is that the functions χ 2 i (S i ) are often well approximated by (up to) quadratic forms in S i ; see also the recent results in [11]. Such forms cover 0νββ decay searches ranging from zero background (with poissonian, linear dependence on S i ) to large backgrounds (with gaussian, quadratic dependence on S i ) [3]. In particular, we have checked that the quadratic approximation works very well also for the latest KamLAND-Zen data [5,14], up to 3σ level at least.
Each experimental result is thus parametrized through a function ∆χ 2 i (S i ) of the form: where the offset c i is set by the condition that the minimum value ∆χ 2 i = 0 is reached within the physical region S i ≥ 0, namely, For a i > 0 the ∆χ 2 i functions are parabolic, with a vertex placed at either S i = 0 (null result, b i = 0), or S i < 0 (negative fluctuation in the unphysical region, b i > 0), or S i > 0 (physical signal or positive fluctuation, b i < 0). In the latter case, the offset c i guarantees ∆χ 2 i = 0 at S i = −b i /2a i . For a i = 0, the ∆χ 2 i functions are linear. For the same isotope, the results of independent experiments are combined by summing their ∆χ 2 i 's, and readjusting the total offset as per Eq. (12). In all cases, 90% C.L. bounds on the half-life (T 90 = 1/S 90 ) are obtained by imposing ∆χ 2 i (S 90 ) = 2.706. Table I, updated from [3] with the inclusion of the latest KamLAND-Zen results [5], reports the coefficients of the parametrization in Eq. (11) and the T 90 bounds for the most sensitive current experiments (T 90 > 0.1), as well as for combinations of experiments using the same isotope. For later purposes, we also consider hypothetical CUORE results for an exactly null signal, denoted as CUORE * (or Te * ). i in terms of the signal strength Si = 1/Ti. The first two columns report the isotope and the names of the experiments or their combination. The next three columns report our evaluation of the coefficients (ai, bi, ci) for the various experiments (upper five rows) and for their combinations in the same isotope (lower rows). The bottom row refers to the case of CUORE sensitivity for null result (tagged by * ). The sixth column reports our 90% C.L. (∆χ 2 = 2.706) half-life limits T90, to be compared with the experimentally quoted one in the seventh column (as taken from the reference in the eighth column, when applicabile).    Table I in graphical form; the left and right panels refer, respectively, to separate experiments and to same-isotope combinations (Xe, Ge, Te). A few remarks about these graphs and the numerics are in order. The case of linear ∆χ 2 functions applies to current GERDA and MAJORANA results, whose combination (denoted as Ge) sets a bound T 90 = 2.07 stronger than for GERDA alone (T 90 = 1.8). All the other experiments are characterized by parabolic functions. KamLAND-Zen and EXO report, respectively, a negative and a positive fluctuation, that partly cancel in their combination (denoted as Xe). As a result, the Xe bound T 90 = 2.26 is slightly weaker (T 90 = 2.26) than for KamLAND-Zen alone (T 90 = 2.3). Note that, for both the Ge and Xe combinations, it is ∆χ 2 i = 0 at S i = 0. Results for the Te isotope depend on the single CUORE experiment, which shows a positive fluctuation at the level of 0.64σ = [∆χ 2 i (0)] 1/2 . As anticipated we consider, besides the real Te results, also hypothetical Te * results, where this fluctuation is canceled by setting b i = 0 (and thus also c i = 0). The half-time limit for Te * (T 90 = 0.301) is in reasonable agreement with the median sensitivity quoted by the CUORE experiment for null result (T 90 = 0.28). In Sec. III, the combination of Xe, Ge, and Te * results (all with ∆χ 2 i = 0 at S i = 0) will provide a simple starting point, before discussing the full combination of Xe, Ge and Te constraints on m ββ .

C. Landscape of nuclear matrix elements
In order to study the combination of 0νββ results in full generality, we consider unconstrained values of the nuclear matrix elements M Xe , M Ge , M Te in the numerical range M i ∈ [0. 2,20]. Within this landscape, we also consider representative M i values from four different approaches to nuclear modeling, including the nuclear shell model (SM), the quasi-particle random phase approximation (QRPA), the energy-density functional theory (EDF), and the interacting boson model (IBM). The NME values are taken from a recent compilation of results [15][16][17][18][19][20][21][22][23][24][25][26] as reported in [1] for the bare value of g A (see Tab. I therein), and are listed in Table II for the sake of completeness. Figure 3 shows the NME landscape in each of the three planes charted by pairs (M i , M j ), together with the representative M i values reported in Table II, which refer to the bare g A . The issue of the effective g A value to be used in nuclear matter, either bare or quenched by a factor q (g A → q g A with q < 1), is largely debated and model-dependent [1,13]. For NMEs dominated by the axial-vector (Gamow-Teller) component (as it is often the case), the leading quenching effect would amount to rescaling the product G i M 2 i by a factor q 4 ∼ g 4 A , that can be assumed to operate on As a representative quenching effect one may consider, e.g., the typical case q g A 1, namely, q ∼ 1/g A , leading to an approximate rescaling M i → M i /g 2 A , as shown in each panel of Fig. 3 by an arrow (applicable to any marked point). Stronger quenching would be associated to longer arrows (not shown). On the other hand, quenching effects would be weaker for sizeable NME vector (Fermi) components, not scaling with q 2 . Moreover, some M i calculations may exhibit a milder dependence on qg A for different reasons. In some QRPA calculations, e.g., 2νββ data are used to adjust the particle-particle parameter g pp , partly trading the effect of quenching g A from its bare value to unity [27]. In the same approach, large statistical covariances are observed among the M i values for different isotopes, inducing noticeable effects on m ββ constraints [27], as recently discussed in [3]. The marked points in Fig. 3 also seem to suggest an overall positive correlation (possibly enhanced by quenching effects) but, since they do not represent a statistical distribution, their covariances (if any) will be ignored. Summarizing, the NME landscape in Fig. 3 is meant to cover a wide and continuous range of M i values, either unquenched or arbitrarily quenched. Graphical results will be shown for unconstrained M i values in this landscape. Marked points in Fig. 3 are meant to represent typical unquenched M i (central values) as taken from the literature (see Table II), while the arrows provide visual guidance for typical quenching effects (q g A ∼ 1).  Table II from

D. Phase space
The last ingredient is represented by the phase space G i for 0νββ decay in i = Xe, Ge and Te, that we take from the calculation in [26]. In our notation and units: The phase space uncertainties [28] are much smaller than those related to 0νββ data and are not considered herein.

III. CONSTRAINTS ON THE MAJORANA NEUTRINO MASS
The previously discussed functions ∆χ 2 i (S i ) can be recast in terms of quadratic functions of m 2 ββ through Eq. (6): where the offset γ i is set by The best-fit value of m ββ is set by: Constraints on m ββ from two or more isotopes are obtained by summing the corresponding ∆χ 2 i , and by adjusting the offset so that it obeys Eq. (17), namely: ∆χ 2 = αm 4 ββ + βm 2 ββ + γ, where α = i α i , β = i β i , and γ = 0 for β ≥ 0 (γ = β 2 /4α otherwise). Bounds on m ββ at a given confidence level are obtained by solving where, e.g., ∆ CL = 2.706, 4 and 9 for limits at 90% C.L., 2σ and 3σ, respectively. Two qualitatively different cases arise from current results: (a) for Xe, Ge and Te * , either separately or in combination, there is no offset γ and the ∆χ 2 function is zeroed at m ββ = 0; (b) for Te results, characterized by γ i > 0 (positive fluctuation), the combination with Xe or Ge (or both) may lead to γ > 0, implying a nonzero Majorana neutrino mass at best fit: m ββ = [−β/2α] 1/2 . We discuss separately these two cases below.
A. Combination of Xe, Ge, Te * constraints In this section we consider current constraints from Xe, Ge, and from the pseudo-experiment Te * (corresponding to null signal in CUORE). In this case the analysis is straightforward, since the best fit is m ββ = 0 for all isotopes and their combinations (independently of the NME). The analysis including real Te data will bring forward cases with m ββ > 0 at best fit, as discussed in the next Section.  We consider both separate and combined Xe, Ge and Te * constraints, as obtained by summing up the corresponding functions defined in Table III, Table II. Concerning constraints from single isotopes, in most cases Xe sets the strongest bound, followed by weaker ones from Ge and Te * . However, for the cases numbered as 9 and 10 (QRPA), the bounds from Xe, Ge and Te * are comparable to each other, and for case 6 (QRPA) the Ge bound actually prevails over the Xe (and Te * ) bound. Notice that such a hierarchy of m ββ constraints may change at different confidence levels, since the S i bounds scale up at different rates (see Fig. 2). In Table IV the combination of pairs of constraints improves appreciably upon each separate constraint; the relative balance in each pair is highlighted below. Finally, the total combination Xe+Ge+Te * provides even stronger bounds on m ββ , that range from a minimum of 38.5 meV (case 12, EDF) to a maximum of 120.4 meV (case 9, QRPA) at 2σ. Figure 4 shows isolines of the 2σ bounds on m ββ , as derived by combining any two pairs among Xe, Ge and Te * , for unconstrained values of the NME. In each panel, single-isotope bounds are asimptotically recovered along each axis, for vanishing matrix element on the other axis. When both matrix elements are sizeable, the joint bound improves upon separate ones. In particular, at each marked point, the bounds in Table IV are recovered for the corresponding NME and pair of isotopes.
For each (M x , M y ) panel and (x, y) isotope pair, the condition for the dominance of one isotope constraint over the other is easily derived. The two isotopes contribute equally to ∆ CL when ∆χ 2 x = ∆ CL /2 = ∆χ 2 y . The solutions to these equations read M x m ββ = ξ x and M y m ββ = ξ y , where ξ x,y are positive numbers. For ∆ CL = 4 it is M Te /M Ge = 1.154, M Ge /M Xe = 2.488, and M Te /M Xe = 2.871, shown as a dashed line in each panel of Fig. 4. Along the dashed line, the two isotopes contribute with equal strength to the 2σ upper bound; above the dashed line, the y-axis isotope dominates over the x-axis one, and viceversa. In this way one gets a graphical interpretation of the hierarchy of bounds for different NME, that was inferred from numerical inspection of Table IV.

B. Combination of Xe, Ge, Te constraints
In this section we consider real Te data, as opposed to the previous cases including Te* pseudo-data. The slight preference of Te data from CUORE for a nonzero signal (as compared with Xe, Ge and Te * , see Fig. 2) brings forward new features of multi-isotope data constraints, although still at embryonic stages.
In general one may expect that, for relatively small values of M Te (with respect to M Xe and M Ge ), the Xe+Ge results will dominate over Te, keeping the best fit at m ββ = 0. However, for increasing M Te , the Te results will eventually prevail and set m ββ > 0 at ∆χ 2 = 0, affecting also upper bounds at some value ∆ CL .
This situation anticipates what could happen with future and more accurate 0νββ data: their combination may (or may not) be consistent with some indications for nonzero m ββ , depending on both the data and the NME values. A future preference for m ββ > 0 might even lead to lower bounds on m ββ , either separately or in combination, depending in part on (un)favorable values of the NME. Eventually, precise multi-isotope data might even test specific NME's by selecting allowed ratios M x /M y [29,30] namely, slanted allowed stripes in the NME landscape of Fig. 3. 3 In our approach, the occurrence of m ββ > 0 at best fit is simply signaled, for a single isotope, by a coefficient β i < 0 in the ∆χ 2 function (currently occurring only for Te) and, for any combination of multi-isotope data, by a negative coefficient β = i β i < 0. The best-fit value of m ββ is then m ββ = (−β/2α) 1/2 with α = i α i , and its specific value depends on the NME's. For β < 0, the offset must be taken as γ = β 2 /4α. In all cases, upper bounds at a chosen confidence levels are set by ∆χ 2 = ∆ CL . Table V shows the numerical results from current Xe, Ge and Te data, regarding the 2σ limits (lower part) and the best-fit values (lower half) of m ββ . In the upper half, the rows corresponding the Xe, Ge and Ge+Xe are unchanged with respect to Table IV, but are repeated for completeness.
Let us first comment on the m ββ best fits in the lower half of Table V. Of course, the Te row displays nonzero results, with m ββ scaling as 1/M Te . In almost all cases including Te (combined with Xe or Ge or both), the positive contributions to β = i β i from i = Xe and Ge are never erased by the negative contribution from Te, and m ββ = 0 is preferred. Only for cases 9 and 10 (QRPA), it turns out that the large ratio M Te 2.6M Xe makes Te prevail over Xe in the corresponding Xe+Te combinations, which show nonzero best fits.  Concerning the upper half of Table V, the combination of Te with Xe (or Ge) does not necessarily improve upon the separate 2σ bounds. Roughly speaking, when the best-fit value of m ββ in Te is comparable or larger than the upper bound from Xe (Ge) alone, the joint bound from Xe+Te (Ge+Te) is weakened, as a result of the slight tension between the two isotopic data. A slight weakening also occurs whenever when Te is added to Xe+Ge in the global combination Xe+Ge+Te. 4 The 2σ bounds on m ββ compiled in Table V The lowest (most optimistic) edge of this range would significantly cut from above the IO and NO allowed regions in Fig. 1, setting also an upper limit on Σ at the level of ∼ 400 meV. Figure 5 shows isolines of the best-fit value of m ββ in the left and right panels charted by (Xe, Te) and (Ge, Te), respectively. In the left panel, the condition i β i > 0 for m ββ > 0 implies M Te /M Xe > 2.42, satisfied only by two marked points (corresponding to the QRPA cases 9 and 10 in Table II)   locate the best-fit values of m ββ as a function of the NME, for each isotope pair. The (in)consistency of the best fits in different pairs will provide interesting clues about the interpretation of data in terms of light Majorana neutrinos. Figure 6 shows isolines of the m ββ upper bounds at 2σ, in the same planes of Fig. 5. In comparison with the lower panels of Fig. 4, a slight weakening of the bounds can be appreciated. Note that, in the presence of subregions where m ββ > 0 at best fit, the offset γ depends on information coming from both isotopes, whose χ 2 contributions cannot be separated in the combination. The condition of equal contributions to ∆ CL cannot be defined in general terms, and the dashed lines of Fig. 4 are thus absent in Fig. 6.
We conclude this Section by discussing the constraints on m ββ at various confidence levels, as derived from the global combination of current (Xe+Ge+Te) data, using the representative NME values in Table II. Since the best fit is m ββ = 0 in all Xe+Ge+Te cases (see Tab. V), only upper bounds need to be quoted. Table VI reports the 90% C.L., 2σ, and 3σ upper limits on m ββ (in meV). Figure 7 shows the N σ = (∆χ 2 ) 1/2 bounds as continuous functions of m ββ . Qualitatively, the strongest limits in are obtained using NMEs from the EDF and IBM models, followed by QRPA and SM cases in mixed order. These results can be generalized to any other choice of NME calculations, using the information provided in this paper. In Fig. 7, one of the N σ curves labelled as QRPA shows a a markedly different (almost linear) slope, intersecting three SM curves. This peculiar curve corresponds to the lowest QRPA point in both panels of Fig. 6 (case 6 in Table II), that is characterized by a rather large value of M Ge as compared with M Xe,Te . As a result, Ge data prevail over Xe and Te in the combination, and the leading dependence is ∆χ 2 ∝ S ∝ m 2 ββ (rather than ∆χ 2 ∝ S 2 ∝ m 4 ββ ), implying a roughly linear N σ (m ββ ) function. Once more, this observation shows the importance of considering the full likelihood profile of the experimental results (e.g., in terms of S = 1/T ), rather than pointlike information (such as the 90% C.L. limit on T ).
As a final step, one could include a joint probability distribution or a ∆χ 2 penalty defined over the NME landscape (M Xe , M Ge , M Te ), and numerically minimize the total ∆χ 2 function. This exercise was performed in [3] by assuming a conservative characterization of the NME and their correlated uncertainties, derived within QRPA calculations [27]. A limit m ββ < 110 meV was obtained at 2σ [3]. By repeating the same exercise with the updated (Xe+Ge+Te) combination considered herein, we get the following marginalized bounds (in meV): m ββ < 79.5 at 90% C.L., m ββ < 99.8 at 2σ, and m ββ < 169 at 3σ. Roughly speaking, from these results and from the summary in Eq. (20) one can state that the combination of current 0νββ experiments sets 2σ upper bounds on m ββ at the level of ∼ 90-100 meV for "average" NME values, possibly lowered to ∼ 40-50 meV for favorable NME values.
A final remark is in order. In principle, one should replace the QRPA input from [27] with more general and up-to-date estimates of the NME's and their uncertainties, characterizing also the spread among different models and calculations. However, no consensus estimates exist yet for NME fiducial values and covariances, although relevant work is in progress toward this goal [1,31]. Part of the planned strategy involves benchmarking nuclear models for 0νββ decay against a variety of data, coming from related electroweak and strong interaction processes or from nuclear structure [32][33][34][35].

IV. CONCLUSIONS AND PERSPECTIVES
We have discussed an approach to the analysis of neutrinoless double beta decay experiments, in terms of ∆χ 2 profiles for the signal strength S i (inverse of the half-life T i ) in the isotopes i = Xe, Ge and Te, building upon previous work [3]. The approach becomes exceedingly simple for quadratic approximations to such profiles, implying quadratic (in)equalities in the landscape of nuclear matrix elements M i that connect the S i to the Majorana mass m ββ . For convenience, some results have been discussed in terms of pseudo data for null signal in Te (dubbed Te * ). Simple relations among the M i have been derived to gauge the relative contributions of different isotopic data in setting upper limits to m ββ (for null best fits in Xe, Ge and Te * ), and to identify the conditions leading to a preference for nonzero m ββ (for generic Xe, Ge and Te data). Using the latest available 0νββ data, as well as representative values of the NME from different models, we have discussed current constraints on m ββ at several confidence levels and in various combinations, both numerically and graphically. Global 2σ upper limits on m ββ are found in the range from 40.5 to 135 meV, depending on the NME.
The approach can be easily extended to nonstandard processes for 0νββ decay [1,36,37] by replacing the relation S i = G i M 2 i m 2 ββ with the appropriate phase space, NME and particle physics parameter characterizing the process. Also, the approach can be extended to generic ∆χ 2 (S i ) functions, with a modest price to pay in terms of numerical (rather than analytical) solutions. We invite the experimental collaborations involved in 0νββ decay searches to publicly provide such ∆χ 2 (S i ) functions or equivalent ones, as they contain much more information than the usually quoted 90% C.L. limits on T i . Indeed, the relative impact of such limits and of the resulting bounds on m ββ in a multi-isotope combination depend sensitively on the likelihood profiles of S i (or T i ), and not only on the relative size of the M i .
Our approach to the multi-isotope data analysis would be complete if one could also assign joint probability densities to the M i , whose variations could then be treated as nuisance parameters and marginalized. So far, detailed results for the NME central values and covariances, including g A quenching uncertainties, have been obtained in a specific (QRPA) model [27]. In perspective, it would be important to extend such investigations to other nuclear models, eventually reaching consensus values for the M i and for their correlated (and possibly reduced) uncertainties.
In this sense, the combined analysis of 0νββ results is proceeding through to the same steps that have characterized similar fields (e.g., solar neutrinos) in the past: from low-statistics data and theoretical models with large uncertainties, to a wealth of accurate experimental results interpreted in increasingly refined and constrained models. Our work aims at providing one methodological step along this path.