A Global Bayesian Analysis of Neutrino Mass Data

TL;DR

This paper conducts a comprehensive Bayesian analysis of neutrino data, integrating oscillation, double beta decay, and cosmology, to assess mass ordering and future experiment prospects.

Contribution

It provides a unified Bayesian framework for neutrino data analysis, evaluating mass ordering and discovery potential with various priors and data combinations.

Findings

Normal ordering indication remains mild.

Discovery potential varies widely with parameters.

Bayesian analysis shows robustness against prior choices.

Abstract

We perform a global Bayesian analysis of currently available neutrino data, putting data from oscillation experiments, neutrinoless double beta decay ($0\nu\beta\beta$), and precision cosmology on an equal footing. We evaluate the discovery potential of future $0\nu\beta\beta$ experiments and the Bayes factor of the two possible neutrino mass ordering schemes for different prior choices. We show that the indication for normal ordering is still very mild and does not strongly depend on realistic prior assumptions or different combinations of cosmological data sets. We find a wide range for $0\nu\beta\beta$ discovery potential, depending on the absolute neutrino mass scale, mass ordering and achievable background level.

Tables

Table 1. Table 1 : Bayes factors for NO vs IO, for all combinations of cosmologies and priors.

Cosmology	Prior on $m_{\mathrm{lightest}}$	Bayes factor NO/IO
restrictive	flat	1.86
conservative	flat	1.64
restrictive	log	1.70
conservative	log	1.61

Table 2. Table 2 : NMEs as extracted from the references given. See text for details.

Method [Reference]	Value for 76Ge	Value for 136Xe
Energy Density Functional Method Rodriguez and Martinez-Pinedo (2010)	$4.60$	$4.20$
Interacting Shell Model Menendez et al. (2009)	$2.30$	$1.77$
Microscopic Interacting Boson Model Barea et al. (2013)	$5.47$	$3.36$
Proton-Neutron QRPA Suhonen and Civitarese (2010)	$5.52$	$3.35$
Self-Consistent Renormalized QRPA (Bonn potential) Meroni et al. (2013)	$5.82$	$3.36$
Self-Consistent Renormalized QRPA (Argonne V18 potential) Meroni et al. (2013)	$5.44$	$2.75$
QRPA (Bonn potential) Simkovic et al. (2013)	$5.40$	$2.38$
QRPA (Argonne V18 potential) Simkovic et al. (2013)	$5.00$	$2.11$
Deformed Self-Consistent Skyrme QRPA Mustonen and Engel (2013)	$5.53$	$1.68$

Table 3. Table 3 : Fit parameters for the cosmological likelihood factor parametrized by Eq. ( 17 ) for both data sets used in our analysis.

Level	Data set	$p_1$	$p_2$	$p_3$	$p_4$
Conservative	base_mnu_plikHM_TT_lowTEB_lensing	$2.073$	$0.401$	$0.552$	$2.514$
Restrictive	base_mnu_plikHM_TT_lowTEB_lensing_BAO	$7.385$	$0.407$	$0.207$	$2.765$

Equations

$$\mathbf{\theta }=\left(m_{\mathrm{lightest}},\Delta m_{\odot }^2,\Delta m_A^2,s_{12}^2,s_{13}^2,{\alpha }_1,{\alpha }_2,\mathcal{G}\right).$$

$${\mathcal{L}}_{\mathrm{glob}}={\mathcal{L}}_{\mathrm{cosmo}}\times {\mathcal{L}}_{0\nu \beta \beta }\times {\mathcal{L}}_{\mathrm{osc}}.$$

$$\Sigma \equiv m_1+m_2+m_3,$$

$${\mathcal{L}}_{0\nu \beta \beta ,\mathrm{i}}=\exp \left[-\left(a_{0,i}+a_{1,i}/T_{1/2}+a_{2,i}/T_{1/2}^2\right)\right],$$

$$\exp \left(-\frac{\Delta {\chi }^2\left({\theta }_i\right)}{2}\right),$$

$$\nu =\frac{N_A\ln 2}{m_{\mathrm{enr}}}\frac{Eϵ}{T_{1/2}},$$

$$P(H|D,I)=\frac{P(D|H,I)P_0(H|I)}{P(D|I)},$$

$$P(D|I)=\sum _i\int d\mathbf{\theta }P(D|M_i,\mathbf{\theta },I)P(\mathbf{\theta }|M_i,I)P(M_i|I)$$

$$\mathbf{\theta }=\left(m_{\mathrm{lightest}},\Delta m_{\odot }^2,\Delta m_A^2,s_{12}^2,s_{13}^2,{\alpha }_1,{\alpha }_2,\mathcal{G}\right)$$

$$P(\log (m_{\mathrm{lightest}}))=const.\iff P(m_{\mathrm{lightest}})\propto 1/m_{\mathrm{lightest}}.$$

$$\frac{1}{T_{1/2}^{0\nu }}=G_{0\nu }|{\mathcal{M}}_{0\nu }{|}^2{\left(\frac{|m_{ee}|}{m_e}\right)}^2,$$

$$\mathcal{P}(1/T_{1/2}^{\mathrm{Ge}})\propto \exp \left[-\frac{1}{2}\frac{(1/T_{1/2}^{\mathrm{Ge}}+1.48{)}^2}{0.461^2}\right].$$

$${\mathcal{L}}_{\mathrm{GERDA}}(T_{1/2}^{\mathrm{Ge}})\propto \exp \left[-\frac{1}{2}\frac{(1/T_{1/2}^{\mathrm{Ge}}+1.48{)}^2}{0.461^2}\right].$$

$${\mathcal{L}}_{\text{KamLAND-Zen}}(T_{1/2}^{\mathrm{Xe}})\propto \exp \left[-\frac{1}{2}\Delta {\chi }^2(T_{1/2}^{\mathrm{Xe}})\right].$$

$$\Delta {\chi }^2=\{\begin{array}{c}2\left[2.3/T_{1/2}^{\mathrm{Xe}}+1.09/(T_{1/2}^{\mathrm{Xe}}{)}^2\right]\mathrm{Phase}-\mathrm{I},\text{hfill}\text{hfill}\text{hfill}\text{hfill}\text{hfill}\\ 2\left[9.71/T_{1/2}^{\mathrm{Xe}}+28.1/(T_{1/2}^{\mathrm{Xe}}{)}^2\right]\mathrm{Phase}-II.\text{hfill}\text{hfill}\text{hfill}\text{hfill}\text{hfill}\end{array}$$

$${\mathcal{L}}_{\mathrm{EXO}}(T_{1/2}^{\mathrm{Xe}})\propto \int e^{-(\nu +\lambda )}\frac{(\nu +\lambda {)}^n}{n!}\mathcal{G}(\lambda ;{\lambda }_0=31.1,{\sigma }_{\lambda }=3.3)d\lambda ,$$

$${\mathcal{L}}_{\mathrm{EXO}}(T_{1/2}^{\mathrm{Xe}})\propto \exp \left[-\frac{1}{2}\frac{(1/T_{1/2}^{\mathrm{Xe}}-0.32{)}^2}{0.30^2}\right].$$

$${\mathcal{L}}_{\mathrm{cosmo}}\left(\Sigma \right)=p_1{\Sigma }^{p_2{p}_3-1}\exp \left[-{\left(\frac{\Sigma }{p_3}\right)}^{p_4}\right].$$

$${\mathcal{O}}_1=\frac{P(H_1|D)}{P(H_0|D)}$$

$$P(H_0|D)$$

$$P(H_1|D)$$

$${\mathcal{O}}_1=\frac{P(D|H_1)}{P(D|H_0)}{\mathcal{O}}_0$$

$${\mathcal{O}}_0=\frac{P_0(H_1)}{P_0(H_0)}$$

$$P(D|H_1)=P(n|H_1)={\int }_0^{\infty }P(n|\lambda +{\nu }^{\prime })P({\nu }^{\prime })d{\nu }^{\prime },$$

$$P(D|H_0)=P(n|H_0)=e^{-\lambda }\frac{{\lambda }^n}{n!}.$$

$$\nu =\frac{N_A\ln 2}{m_{\mathrm{enr}}}\frac{Eϵ}{T_{1/2}},$$

$$P_{\mathrm{discovery}}$$

$$I$$

$$P(n|H_1)$$

$$=$$

Full text

This is an updated version of the original paper dated 4 May 2017. The update contains a correction to an error found in our use of the PLANCK collaboration data. All relevant figures were updated but no new data were included.

††preprint: MPP-2017-79

I Introduction

Neutrino physics is one of the most attractive fields to look for new physics, and many parameters in the neutrino sector are not yet determined. A series of oscillation experiments has established the fact that at least two neutrino masses are distinct from zero, but their smallness cannot be accommodated within the Standard Model of Particle Physics. Furthermore, quantities like the Dirac phase $\delta_{\mathrm{CP}}$ describing the difference between matter and antimatter, the absolute neutrino mass scale, or the mass ordering (i.e., which neutrino is the lightest) are currently only poorly restricted Esteban et al. (2017). The even more fundamental question of whether neutrinos are of Majorana or of Dirac nature (i.e., whether they are identical to their antiparticles, or not) is also still unanswered. A global analysis, combining all relevant experimental results and using the current data as efficiently as possible, is the most suitable way to address these open issues. It additionally informs upcoming experimental choices.

The first fully comprehensive Bayesian analysis of this kind is presented in this letter. We use information from oscillation experiments, precision cosmology, and neutrinoless double beta decay ($0\nu\beta\beta$). While previous works focus mostly on one of these aspects, we adopt a fully global view. Our analysis is based on the minimal framework of three light Majorana neutrinos, the most predictive setting for neutrino physics. We use global oscillation data from the nu-fit collaboration Esteban et al. (2017), cosmological data from the Planck Legacy Archive (PLA),111Based on observations by Planck (http://www.esa.int/Planck). and data from the $0\nu\beta\beta$-experiments KamLAND-Zen Gando et al. (2016a), EXO-200 Albert et al. (2014a), and GERDA Agostini et al. (2017a). We do not include single $\beta$ decay results as they currently do not provide additional constraints.

We extract the implications for future $0\nu\beta\beta$-experiments and also address whether current data may already exhibit a tendency towards the normal ordering scheme – which is currently under intense debate Gerbino et al. (2016); Capozzi et al. (2017); Simpson et al. (2017); Schwetz et al. (2017).

Bayes’ Theorem provides the logical path from the probability of data under different hypotheses $H_{i}$, $P\left(D|H_{i}\right)$, to the probability of a hypothesis being correct given the data. It requires the explicit definition of prior probabilities. Any other approach is either incoherent (for a discussion of $p$-values as evidence, see e.g. Schervish (1996)), requires predefined error rates (Neyman frequentist approach Neyman (1977)), or contains implicit prior choices – all of which are to be disfavored. If the outcome depends strongly on the choice of prior, this is simply an indication that the data is not powerful enough to draw a conclusion. In this context, it is important to find the implications from a range of reasonable prior choices.

For our analysis we assign equal prior probabilities to the two mass orderings. For the absolute neutrino mass scale, the choice of the prior on the mass, $m_{\rm lightest}$, of the lightest neutrino is critical. At present, only an upper limit is available. We consider two choices: A prior flat in $m_{\rm lightest}$ (this will tend to favor larger values and reflects the sizes of the measured solar and atmospheric mass square differences $\Delta m^{2}_{\odot}\equiv m_{2}^{2}-m_{1}^{2}$ and $\Delta m^{2}_{A}\equiv|m_{3}^{2}-m_{1}^{2}|$) and a scale-free prior (flat in the logarithm of $m_{\rm lightest}$, which tends to favor small values and reflects the fact that the absolute neutrino mass scale is still unknown). We then evaluate the consequences of both prior choices on the posterior probabilities.

Technical details regarding priors and likelihoods are provided in the supplementary material.

II Parameter Space, Likelihoods, Priors

Parameter space

Our analysis depends on eight parameters, aggregated into the vector $\bm{\theta}$:

[TABLE]

Here, $m_{\mathrm{lightest}}$ is the smallest neutrino mass eigenvalue, while $s^{2}_{ij}$ denotes the sine squared of the mixing angle $\theta_{ij}$. The Majorana phases $\alpha_{1}$ and $\alpha_{2}$ take trivial values if neutrinos are of Dirac nature. The nuclear matrix elements (NMEs) required to calculate rates of $0\nu\beta\beta$ for different considered isotopes are condensed into $\mathcal{G}$. Note that $\bm{\theta}$ does not contain the Dirac phase $\delta_{\mathrm{CP}}$ or $s^{2}_{23}$, as they only affect neutrinos oscillations, but neither $0\nu\beta\beta$ nor cosmology. Furthermore we allow for the choice of either normal (NO: $m_{1}<m_{2}<m_{3}$) or inverted (IO: $m_{3}<m_{1}<m_{2}$) ordering.222Such an “external” parameter is called hyperparameter.

Likelihood

The global likelihood function $\mathcal{L}_{\mathrm{glob}}$ is a function that depends non-trivially on all eight parameters. Since every suite of experiment constitutes an independent data set, we can factorize the likelihood:

[TABLE]

The cosmological term $\mathcal{L}_{\mathrm{cosmo}}$ only depends on the sum of neutrino masses,

[TABLE]

and hence on $m_{\mathrm{lightest}}$, $\Delta m^{2}_{\odot}$, and $\Delta m^{2}_{A}$. We reconstructed the likelihood from the Markov chains that are publicly available on the PLA. As those chains were sampled on a prior flat in $\Sigma$, one can construct a function of $\Sigma$ that is directly proportional to $\mathcal{L}_{\mathrm{cosmo}}$ from these chains.333Note that the choice of prior used by the Planck collaboration is, in fact, unphysical, since $\Sigma$ has a minimal value $\Sigma_{\rm min}\neq 0$. As shown in Gerbino et al. (2016), using a more realistic prior, the resulting constraints become stronger, but they constrain $(\Sigma-\Sigma_{\rm min})$, making the corresponding limit on $\Sigma$ appear looser. To be agnostic about the choice of datasets used in our analysis, we chose two models: both contain the $TT$ correlation of the Cosmic Microwave Background (CMB), a joint likelihood on $TT$, $EE$, $BB$, and $TE$ correlations444Here, $T$ refers to the temperature of the CMB, while $E$ and $B$ denote the basic patterns of its polarization. for low multipoles ($2\leq l\leq 29$), as well as weak gravitational lensing data. We refer to this combination of data as conservative model. Augmenting the data by data from Baryonic Accoustic Oscillations (BAO) defines a restrictive model, which puts tighter constraints onto $\Sigma$.555We use the PLA data sets base_mnu_plikHM_TT_lowTEB_lensing (conservative) and base_mnu_plikHM_TT_lowTEB_lensing_BAO (restrictive).

The term $\mathcal{L}_{\mathrm{0\nu\beta\beta}}$ is factorized into three terms for GERDA, KamLAND-Zen, and EXO, respectively. Each term is well parametrized by the form

[TABLE]

where $T_{1/2}$ is the neutrinoless double beta decay half-life in units of $10^{25}\;\mathrm{a}$. The parameters entering this formula have either been determined directly from published data or confirmed by the experimental collaborations.

The neutrino oscillation data is adopted from the global oscillation fit provided by nu-fit.org Esteban et al. (2017). This collaboration does not provide a public likelihood function depending on all oscillation parameters. We thus approximate the likelihood using a spline-fit of the 1-dimensional projections of their $\Delta\chi^{2}$ fit onto the relevant parameters, shifting the distributions so they have a minimum at $\Delta\chi^{2}=0$, and taking the likelihood to be a product of the resulting four factors. Each such factor is given by

[TABLE]

where $\chi^{2}\left(\theta_{i}\right)$ is the one dimensional projection onto the parameter $\theta_{i}$. To account for the slight tendency of oscillation data to favor NO, we then added one offset $\Delta\chi^{2}=0.83$ Esteban et al. (2017) to the product value for IO after re-calibrating the projections. This procedure avoids counting the slight inclination towards NO each time it gets projected onto one parameter. The quoted value of $\Delta\chi^{2}$ corresponds to a Bayes factor of $1.5$ in favor of NO.666In future analyses, this procedure could be refined by using the full multiparameter likelihoods or by directly fitting the oscillation data.

Priors

Our choice of priors is intended to be informative whenever possible, and otherwise agnostic for parameters with little or no information available. For mass square differences and mixing angles, we choose flat priors extending well beyond the best-fit values for both orderings, since the likelihood functions severely constrain their values in any case. The choice of range does not influence our results. For the Majorana phases $\alpha_{1}$ and $\alpha_{2}$, we choose priors flat on the entire range $\left[0,2\pi\right]$, as we do not have any prior knowledge or theoretical motivation for an informative distribution. The prior on the NME factor $\mathcal{G}$ of a certain isotope is discrete, equally weighing the different theoretical calculations, which deviate from one another by factors of up to three. This procedure will be further explained and justified in the Appendix (cf. Fig. 4 and subsequent discussion).

We choose two very different priors for $m_{\mathrm{lightest}}$, to explicitly show how results depend on the respective assumption. Both choices span the range $\left[10^{-7}\;\mathrm{eV},0.6\;\mathrm{eV}\right]$.777The upper end of our range is motivated by the restrictive power of the cosmological data sets: Extending it to larger values does not alter the results significantly. The lower end of the range can be argued to be as low as $10^{-13}\;\mathrm{eV}$ Davidson et al. (2007), using a general theoretical argument from perturbation theory. As we are interested in implications for $0\nu\beta\beta$, there will however be a point where lowering the minimal value does not anymore alter any implications. The flat prior allocates 90% of the probability mass at values greater than $60\;\mathrm{meV}$, while the scale-free log prior allocates about 85% in a region where $m_{\mathrm{lightest}}<60\;\mathrm{meV}$.

III Results of our analysis

For all combinations of neutrino mass ordering, cosmological data set, and prior for $m_{\mathrm{lightest}}$, we used the BAT Markov chain Monte Carlo code in version v2.0.4 Schulz et al. (2021) to sample from the posterior probability density. For each scenario, we ran 4 Markov chains with a length of $1\times 10^{6}$ samples each. The results derived from the posterior samples are presented in figs. 1, 2, 3.

Posterior probabilities

For fixed mass ordering, cosmology, and prior on $m_{\mathrm{lightest}}$, we show the posterior probability density in the plane of the effective mass $|m_{ee}|$ controlling the $0\nu\beta\beta$ rate and the sum $\Sigma$ of neutrino masses in fig. 1. We also show theoretical envelopes bounding the regions in the $\Sigma$–$|m_{ee}|$ parameter space if the Majorana phases are varied within their full ranges while the relevant oscillations parameters are kept at their best fit values (dashed lines) or varied within their $3\sigma$ ranges (solid lines).

Choosing the restrictive cosmological data sets forces $\Sigma$ (and by means of oscillation data also $m_{\mathrm{lightest}}$) to be smaller, resulting in more probability mass shifted towards smaller values of $|m_{ee}|$. Another clear feature is the fact that a scale-free prior on $m_{\mathrm{lightest}}$ leaves almost no probability mass in regions where $\Sigma\gtrsim 0.1\;\mathrm{eV}$. Finally, we note that the highest probability density tends to be at the upper ends of the allowed $|m_{ee}|$ ranges for a given $\Sigma$. This reflects the flat prior assumption on the Majorana phases, which makes it highly unlikely to get a cancellation of terms leading to vanishing $|m_{ee}|$ Lindner et al. (2006); Merle and Rodejohann (2006); Maneschg et al. (2009).

Double Beta Decay

The posterior samples from our global fits can be used to infer the discovery potential for $0\nu\beta\beta$ experiments. The relevant quantity is the half-life probability distribution, examples of which are shown in fig. 2.888The full set of plots is available in the supplementary material. As expected, the probability mass at smaller half-life is always bigger in the case of IO. However, for the log prior (right panel), favoring small $m_{\mathrm{lightest}}$, the effect is much more pronounced than in the case of the flat prior (left panel).

The relationship between the signal expectation and the half-life is

[TABLE]

where $N_{A}$ is Avogadro’s number and $m_{\mathrm{enr}}$ is the molar mass of the relevant enriched isotope. $E$ is the exposure and $\epsilon$ is the efficiency to find a signal in the region of interest. The latter includes active/fiducial mass considerations as well as signal reconstruction efficiencies.

Selected results for the discovery potential for $0\nu\beta\beta$, assuming light Majorana neutrinos as the source of a signal, are presented as a function of effective exposure $E\epsilon$ for four different background levels $({\rm BG}=0.1,1,10,100\leavevmode\nobreak\ {\rm events})$ for 76Ge and three for 136Xe $({\rm BG}=1,10,100)$999Discovery is defined as a Bayes factor of at least $100$ in favor of the presence of $0\nu\beta\beta$. The details of the statistical analysis are available in the supplementary material. in fig. 3. The discovery potential can be used to judge the merits of individual experiments. In general, we see that a prior flat in $m_{\mathrm{lightest}}$ gives a much larger probability for discovery than a prior flat in $\log(m_{\mathrm{lightest}})$. The importance of keeping the background small is obvious from these plots.

Mass ordering

As noted earlier, the oscillation data give a Bayes factor of 1.5 in favor of the normal mass ordering. We have evaluated the Bayes factor that results from including the data of cosmology and $0\nu\beta\beta$ into the analysis for the different cosmological data sets used as well as for our two choices of priors on $m_{\mathrm{lightest}}$. All evidence integrals were computed using the Cuhre quadrature algorithm provided by the CUBA library Hahn (2005). The results are shown in Table 1. As can be seen, the additional data and our choices of priors do not change the conclusions significantly. Also the choice of the cut-off scale at $10^{-7}$ eV is perfectly reasonable, as lowering it does not change the numerical results significantly anymore.

IV Conclusions

In this work, we have conducted a global Bayesian analysis of neutrino mass parameters within the minimal framework of three light Majorana neutrinos, combining data from oscillation experiments, $0\nu\beta\beta$ decay, and precision cosmology. Working with one prior flat in $m_{\mathrm{lightest}}$ and another one flat in $\log\left(m_{\mathrm{lightest}}\right)$, we have investigated both extreme cases of uninformative priors for the critical parameter $m_{\mathrm{lightest}}$. Combining these priors with a more conservative cosmological data set on one hand and a more restrictive one on the other hand, we conclude that the posterior probability for NO is still very mild, even in the extreme case of a restrictive cosmology and a flat prior on $m_{\mathrm{lightest}}$. In the other cases, the slight inclination towards NO is almost entirely driven by neutrino oscillation experiments.

Furthermore we have evaluated the posterior distributions of $T_{1/2}$ in the different combinations of priors and cosmological data sets as well as for different isotopes. This allowed us to infer the discovery potential of different experimental approaches for $0\nu\beta\beta$. Depending on the neutrino mass ordering, the achievable effective exposure, and the background level, the discovery potential spans a wide range. Assuming a flat prior for $m_{\mathrm{lightest}}$ is in all cases favorable for $0\nu\beta\beta$ searches as is the inverted mass ordering.

Our approach of putting all relevant experimental insights on a consistent and equal footing helps to pave the way for future comprehensive analyses of neutrino data. Additional data can be easily accommodated and probabilities updated in this scheme.

Acknowledgments

We would like to thank Matteo Costanzi, Matteo Viel, and in particular Massimiliano Lattanzi for valuable discussions concerning cosmology. Furthermore we are grateful to Itaru Shimizu and Luciano Pandola for checking our likelihood parametrizations for KamLAND-Zen and GERDA, respectively. We thank Frederik Beaujean for advice regarding the implementation of our model in BAT. AM acknowledges partial support by the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreements No. 690575 (InvisiblesPlus RISE) and No. 674896 (Elusives ITN), as well as by the Micron Technology Foundation, Inc. MT acknowledges support by the IMPRS-EPP and by the Studienstiftung des deutschen Volkes.

Supplementary Material

This supplementary material describes how we obtained the likelihood functions, it illustrates the procedure we used to determine discovery potentials, and it also presents our main result (namely the posterior probability densities) in an alternative coordinate system, which only uses observables accessible in laboratory experiments.

Statistical Formulation

The analysis was performed using Bayes’ Theorem:

[TABLE]

where $H$ is the hypothesis to which we are assigning a probability, $D$ represents the data being used to update the probability, and $I$ is the additional information and constraints that have been used to define the probabilities. $P_{0}(H|I)$ is the prior probability assignment to $H$. The denominator can be rewritten using the Law of Total Probability as

[TABLE]

where $M_{i}$ represent different models allowed under $I$ subject to $\sum_{i}P(M_{i}|I)=1$ and $\bm{\theta}$ are the (possibly different) parameters present in the model. $P(\bm{\theta}|M_{i},I)$ is also normalized.

In our analysis, we assume that there are three light Majorana neutrinos, and have two models representing the two allowed mass ordering schemes: normal (NO: $m_{1}<m_{2}<m_{3}$) or inverted (IO: $m_{3}<m_{1}<m_{2}$). We assign equal prior probabilities to the two mass orderings: $P(M_{1}={\rm NO}|I)=P(M_{2}={\rm IO}|I)=1/2$. The parameters used are

[TABLE]

where $m_{\mathrm{lightest}}$ is the smallest neutrino mass eigenvalue and $s^{2}_{ij}$ denotes the sine squared of the mixing angle $\theta_{ij}$. The Majorana phases $\alpha_{1}$ and $\alpha_{2}$ take trivial values if neutrinos are of Dirac nature. The nuclear matrix elements (NMEs) required to calculate rates for $0\nu\beta\beta$ are collectively denoted by $\mathcal{G}$. Note that $\bm{\theta}$ does not contain the Dirac phase $\delta_{\mathrm{CP}}$ or $s^{2}_{23}$, as they only affect neutrinos oscillations but not $0\nu\beta\beta$ or cosmology. We assign flat probabilities to all parameters except for $m_{\mathrm{lightest}}$ and $\mathcal{G}$. The latter are discussed in more detail below. For the parameters with flat priors, the ranges were defined either to be the maximal range (for the phases) or to cover a sufficiently wide range so as not to affect the results. For example, in the case of the mixing angles and squared mass differences, the parameters are so strongly constrained by the oscillation data that values even remotely close to the edges are already strongly disfavored.

We consider two prior probability assignments for $m_{\mathrm{lightest}}$: A prior flat in $m_{\mathrm{lightest}}$ (this will tend to favor larger values and it reflects the sizes of the measured mass square differences $\Delta m^{2}_{\odot}\equiv m_{2}^{2}-m_{1}^{2}$ and $\Delta m^{2}_{A}\equiv|m_{3}^{2}-m_{1}^{2}|$) and a scale-free prior (flat in the logarithm of $m_{\mathrm{lightest}}$, which tends to favor small values and reflects the fact that the absolute neutrino mass scale is still unknown). We note that

[TABLE]

The prior probability for $\mathcal{G}$ was defined by assigning equal probabilities of $1/9$ to the different matrix element computations, the results of which are collected in Table 2, using the normalization corresponding to the following half-life formula Doi et al. (1985):

[TABLE]

where $G_{0\nu}$ is a phase space factor which can be easily computed for any isotope under consideration (we make use of the values from Ref. Suhonen and Civitarese (1998), which were slightly updated by Ref. Rodejohann (2011)). $\mathcal{M}_{0\nu}$ is the NME, which encodes all nuclear physics that goes into the process. Some remarks are at order concerning the values shown in Tab. 2. For definiteness, we have adopted the standard value $g_{A}=1.25$ for the axial vector coupling, which means that we had to rescale some of the NMEs Menendez et al. (2009); Simkovic et al. (2013). Also the phase space factors used Suhonen and Civitarese (1998); Rodejohann (2011) correspond to $g_{A}=1.25$. In cases where different versions of a computations are available, we have for definiteness always chosen the most optimistic result. For example, Ref. Meroni et al. (2013) reports the results for both, an intermediate size model space and a large size single particle space, the latter of which tends to yield larger values of the NME. Hence, we have decided to use the large size results. This differs from the choices e.g. made in Ref. Bhupal Dev et al. (2013), where in some cases the smaller and in others the larger value has been chosen. Nevertheless, none of the treatments is wrong in the sense that at the moment the values still suffer from nuclear physics uncertainties, and there exists no way to decide which value is closer to reality.

The NMEs are shown together with their posterior probabilities in Fig. 4. While we only show results for 76Ge and 136Xe, it is straightforward to extend the analysis to any isotope for which NMEs have been calculated. Note that the posterior distribution for the NMEs tracks the prior distribution quite closely. This is to be expected, as the data from the $0\nu\beta\beta$ experiments does not yet hint towards a signal but only gives lower bounds on the lifetimes of the isotopes. Thus, no insights about the NMEs mediating the decays can be generated, and our prior knowledge does not get updated in a relevant way. Due to this argument, we can also justify taking a prior distribution sharply centered around discrete values. Even if each calculation were equipped with an uncertainty, which could be accounted for by modeling the prior by a sum of narrow Gaussians, this would not change the results significantly. As the spread between the different calculations is presumably much larger than this uncertainty, the overall picture would stay unaltered. The choice of assigning a equal normalization of $N^{-1}_{calc}$ to each calculation just formalizes our current state of ignorance concerning the best method—a state that can only be changed by observation of $0\nu\beta\beta$ decay.

The probability of the data, $P(D|M_{i},\bm{\theta},I)$, viewed as function of the parameters with the data fixed, is known as the likelihood $\mathcal{L}(\bm{\theta})$. In the following, we describe which prescriptions were used to define these likelihoods.

Neutrino Oscillation Data

The neutrino oscillation likelihood is taken from nu-fit.org Esteban et al. (2017), of which we have used v3.0 of their neutrino oscillation global fit results. As described in the main text, we have converted the $\Delta\chi^{2}$ projections from nu-fit.org into likelihood functions, $\mathcal{L}(\theta_{i})\propto\exp\left[-\Delta\chi^{2}(\theta_{i})/2\right]$, where $\Delta\chi^{2}(\theta_{i})$ is the $\Delta\chi^{2}$ projection for the parameter $\theta_{i}$. The likelihoods are shown in Fig. 5. Note that the fit by nu-fit.org very slightly disfavors IO, by an offset $\Delta\chi^{2}=0.83$. However, while this value appears in each of the projections shown by that collaboration, we have to be careful not to penalize IO more than once in our analysis, because the value of $0.83$ also appears only once in the full $\chi^{2}$-function, before projecting out individual parameters. We have accounted for this by adding the offset only once, to avoid an artificial worsening of IO which would even increase with the number of parameters used. Note that this results in a global factor in front of the likelihood function for IO, which drops out when computing posterior probabilities. However, once we compute a Bayes Factor, both likelihoods need to share a common normalization to make sense of the values obtained, and in that case this constant prefactor does play a role, and it corresponds to the factor of $1.5$ in favor of NO, as mentioned in the main text.

Double Beta Decay data

The most strongly constraining experiments to date are the GERDA, KamLAND-Zen, and EXO experiments, and we have included likelihood terms for these experiments in our analysis. The data were analyzed and presented in different ways by the three experiments. We discuss these in turn to explain how we included these data into our Bayesian analysis. Note that $T_{1/2}$ values are always given in units of $10^{25}$ years. The functions we will present in the following are plotted in the left panel of Fig. 6 for illustration.

GERDA

The GERDA Collaboration Agostini et al. (2017b) performed both a Bayesian and a frequentist limit setting analysis on their data. We use the posterior probability density for the inverse of the half-life of 76Ge as determined by the GERDA Collaboration, which is well parametrized by the form101010The parametrizations have been cross-checked by the collaboration (L. Pandola, private communication).

[TABLE]

Since GERDA used a flat prior for $1/T^{\rm Ge}_{1/2}$, we use

[TABLE]

KamLAND-Zen

The KamLAND-Zen collaboration Gando et al. (2016b) performed frequentist analyses based on a profile likelihood test statistic. They presented their results in the form of $\Delta\chi^{2}$ as a function of $T^{\rm Xe}_{1/2}$ (Wilk’s Theorem was assumed to apply). We have interpreted these $\Delta\chi^{2}$ as likelihoods for the half-life of 136Xe by setting

[TABLE]

The Phase-I and Phase-II results were used separately, and the following parametrizations were found to represent the KamLAND-Zen results well:111111The parametrizations have been cross-checked by the collaboration (I. Shimizu, private communication).

[TABLE]

EXO

The EXO Collaboration also performed a frequentist analysis and presented confidence level limits on $T^{\rm Xe}_{1/2}$. However, the background expectation, background uncertainty, and observed numbers of events were presented in the paper Albert et al. (2014b), allowing for the construction of the likelihood as

[TABLE]

where the signal expectation $\nu$ depends on $T^{\rm Xe}_{1/2}$ as well as on the exposure and the efficiency, which are both given in Ref. Albert et al. (2014b). The background level is given by $\lambda$ and is centered on $\lambda_{0}$. The uncertainty is introduced by smearing the expectation with a Gaussian of width $\sigma_{\lambda}$. The observed number of events is $n=39$. The resulting likelihood is well parametrized by the form

[TABLE]

We note that the peak is at finite $T^{\rm Xe}_{1/2}$ since more events were observed than expected from background processes.

Cosmological Data sets

The cosmological factor $\mathcal{L}_{\rm cosmo}$ contains extensive information about cosmic structure formation and cannot be given in a simple closed form. Instead, we have used Markov chains sampling the posterior for different combinations of datasets that are publicly available on the Planck Legacy Archive (PLA).121212Based on observations by Planck (http://www.esa.int/Planck). These posterior samples were obtained from a uniform prior in $\Sigma$, such that the posterior is identical to the likelihood except for some unknown normalization that is irrelevant in our implementation. The total number of samples for each model is about $50\times 10^{3}$. The samples were binned into 50 bins to obtain a likelihood estimator, which was then fit to the following four-parameter template:

[TABLE]

This functional form is related to the generalized Gamma distribution and it is a good empirical fit to likelihood, as the right panel of Fig. 6 shows. Tab. 3 lists the fit coefficients we obtained for the two cosmological data sets used in our analysis.

Double Beta Decay: Discovery Probability Analysis

We define two hypotheses:

• $H_{0}$

There are only background processes producing the data;

• $H_{1}$

There is, additionally to the background processes, also a signal from neutrinoless double beta decay.

The full information is contained in the posterior probability for the hypothesis given the data, $P(H_{1}|D)$. If this probability is large enough, then a ‘discovery’ can be claimed. We will use the related concept of posterior odds

[TABLE]

and place requirements on this quantity to claim a discovery. Assuming the two hypotheses are exhaustive131313We exclude here the possibility that there are background sources not included in the background model expectations. We also assume that the signal would come from the exchange of a light Majorana neutrino. we have

[TABLE]

The posterior odds are then given by

[TABLE]

where

[TABLE]

are the prior odds. The factor $\frac{P(D|H_{1})}{P(D|H_{0})}$ is known as the ‘Bayes Factor’.

We consider the case of a single experiment measuring $n$ events in a given region of interest defined by the $Q$ value of the decay, the background level, and the energy resolution. We use simple Poisson statistics to evaluate the discovery potential. Small improvements are possible by making use of the background and signal spectral shapes, but the additional sensitivity gained by doing this is not significant for our discussion. We will also ignore uncertainties in the background and assume it is well known. The uncertainty in the background typically has little effect on the discovery calculation if the background is small, but it can become important if this is not the case.

We use the symbol $\lambda$ to represent the background expectation in the region of interest, and $\nu$ is the signal expectation given a half-life for the decay, the exposure, and the detection efficiency.

Given these definitions, we have

[TABLE]

where $P(\nu^{\prime})$ is the probability to have the signal expectation $\nu^{\prime}$ given the experimental conditions and prior knowledge on $T_{1/2}$. The probability of the data for $H_{0}$ is simply

[TABLE]

The relationship between the signal expectation and the half-life is

[TABLE]

where $N_{A}$ is Avogadro’s number and $m_{\rm enr}$ is the molar mass of the relevant enriched isotope. $E$ is the exposure and $\epsilon$ is the efficiency to find a signal in the region of interest. The latter includes active/fiducial mass considerations as well as signal reconstruction efficiencies.

We evaluate the probability of a claim of discovery assuming three light Majorana neutrinos by:

[TABLE]

where $I$ is the indicator function:

[TABLE]

In addition to the exposure and the signal efficiency, the discovery potential depends on the background level and on the cut value for the posterior odds. We take the prior odds to be $1$. The probability distribution for $\nu$ is given by the results of our neutrino parameter analysis, and depends on the mass orderings, on the probability distribution for the effective neutrino mass in these orderings, and on the matrix element values. We give the probability distributions for $T_{1/2}$ for 76Ge and 136Xe for the different orderings, cosmological data sets sets, and mass scale priors in Fig. 7.

We use the output of BAT to carry out the integrals:

[TABLE]

with $N$ being the number of samples of $T_{1/2}$ taken from the Markov chain. Recall that the background level $\lambda$ is fixed.

The discovery potential is presented as a function of efficiency reduced exposure, $E\epsilon$, for four different background levels $(\lambda=0.1,1,10,100)$ for Ge and three for Xe $(1,10,100)$. The results are shown in Figs. 8 and 9.141414Results for other isotopes can be requested from the authors.

The cut value for the Bayes Factor was taken to be $100$ – i.e., we would claim a discovery if $H_{1}$ was 100 times more probable than $H_{0}$. This is an arbitrary choice, and varying it by a factor $10$ is certainly possible, but it does not have strong influence on the results of this analysis. The jumps in the curve are due to the Poisson nature of the statistical fluctuations.

Posterior probability for lab-based experiments

In the main text, we have shown the heat map of posterior probability density in the $\Sigma$–$|m_{ee}|$ plane, i.e., in terms of the observables that are probed by cosmology and by $0\nu\beta\beta$. This choice is inspired by these two observables being the most promising both to identify the neutrino mass ordering and to get information on the absolute neutrino mass scale. In the theoretical literature, one often finds plots of $m_{\rm lightest}$ vs. $|m_{ee}|$, see e.g. the literature on how to probe certain neutrino mass models by $0\nu\beta\beta$ King et al. (2013); Gehrlein et al. (2016, 2015); Agostini et al. (2016). However, this choice is problematic in our case because not only is $m_{\rm lightest}$ only a theoretical parameter rather than an actual observable, but it also exhibits a very strong dependence on the prior, which is poorly constrained.

When considering laboratory-based experiments, there is a long history of attempts to measure the neutrino mass in single beta decay experiments. These types of experiments are looking for deviations in the kinematic endpoint of the resulting electron spectrum, and what they constrain (at least with realistic resolution) is the so-called effective electron neutrino mass $m_{\beta}$, defined by151515Note that, in fact, there is no mass that could be associated with an electron neutrino, as this is a not a mass eigenstate. Thus, the $m_{\beta}$ defined in Eq. (24) is really nothing more than an effective quantity, however, an observable one.

[TABLE]

Constraints on $m_{\beta}$ have mainly been derived from tritium decay, by both MAINZ Kraus et al. (2005) and Troitsk Lobashev et al. (1999), however also future proposals focus on tritium161616Alternative isotopes, such as rhenium as used by MARE Ferri et al. (2015), do not seem to be competitive. (such as Project 8 Monreal and Formaggio (2009), which aims at detecting cyclotron radiation from relativistic tritium ions). The most prominent on-going effort is provided by KATRIN Osipowicz et al. (2001), which is data taking and whose sensitivity on $m_{\beta}$ is of $\mathcal{O}(1\leavevmode\nobreak\ {\rm eV})$, with quoted values being dependent on the statistics used (e.g., Ref. Host et al. (2007) predicts a sensitivity of 0.20 eV based on frequentist methods and a sensitivity of 0.17 eV when using Bayesian statistics with a flat prior $m_{\beta}^{2}$).

The heat maps for our posterior probability densities for all combinations used of orderings and cosmological limits are displayed in Fig. 10. The underlying data are equivalent to those shown in the main text. As can be seen, the generic tendencies are basically the same as for the $\Sigma$–$|m_{ee}|$ plot. However, even if the ultimate sensitivity limit of KATRIN is reached, it will not add a significant constraint on the parameter space for the scenarios investigated. This is our justification for neglecting data from single beta decay experiments.

Conclusion

The probability density heat maps and $0\nu\beta\beta$ results are based on currently available data. They should be updated as new data becomes available.