Abstract
To quantitatively assess the impact of an eVmass sterile neutrino on the neutrinoless doublebeta () decays, we calculate the posterior probability distribution of the relevant effective neutrino mass in the (3+1) mixing scenario, following the Bayesian statistical approach. The latest globalfit analysis of neutrino oscillation data, the cosmological bound on the sum of three active neutrino masses from Planck, and the constraints from current decay experiments are taken into account in our calculations. Based on the resultant posterior distributions, we find that the average value of the effective neutrino mass is shifted from (or ) in the standard 3 mixing scenario to (or ) in the (3+1) mixing scenario, with the logarithmically uniform prior on the lightest neutrino mass (or on the sum of three active neutrino masses). Therefore, a null signal from the future decay experiment with a sensitivity to will be able to set a very stringent constraint on the sterile neutrino mass and the activesterile mixing angle.
Impact of an eVmass sterile neutrino on the neutrinoless doublebeta decays: a Bayesian analysis
Guoyuan Huang ^{*}^{*}*Email: huanggy@ihep.ac.cn, Shun Zhou ^{†}^{†}†Email: zhoush@ihep.ac.cn
Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China
School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
1 Introduction
Whether massive neutrinos are Majorana or Dirac particles is one of the most important problems in particle physics [1, 2, 3]. Quite a number of neutrinoless doublebeta () decay experiments are devoted to answering this question [4, 5, 6, 7, 8, 9, 10]. If massive neutrinos are Majorana particles and thus lepton number violation exists in nature, then the decays could take place in some eveneven nuclei, namely, both the proton number and the neutron number for the nuclear isotope are even. Assuming the exchange of light Majorana neutrinos to be responsible for the decays, one can find that the halflife of the relevant nuclear isotope is given by [5]
(1) 
where is the phasespace factor, is the nuclear matrix element (NME), and is the electron mass. In Eq. (1), the effective neutrino mass collects the contributions from light Majorana neutrinos involved in the decays.
In the standard three neutrino () mixing scenario, the effective neutrino mass is defined as , where the absolute neutrino masses and the lepton flavor mixing matrix elements (for ) appear. When the conventional parametrization of the flavor mixing matrix is adopted [3], i.e., , and , we have
(2) 
where are two of three neutrino mixing angles, and are the Majoranatype CPviolating phases. Note that is nonzero no matter whether the normal neutrino mass ordering (NO) with or the inverted neutrino mass ordering (IO) with is considered. Therefore, such a parametrization is favorable in the discussions about the limiting case of (for NO) or (for IO), for which one of two Majoranatype CP violating phases just disappears together with the lightest neutrino mass.
However, if the eVmass sterile neutrino indeed exists as a solution to the anomalies in the shortbaseline neutrino experiments [11, 12, 13, 14, 15, 16, 17, 18], it will contribute as well to the decays. In this case, the effective neutrino mass is given by with being the mass of the sterile neutrino and (for ) being the firstrow elements of the mixing matrix in the (3+1) mixing scenario. Adopting the standard parametrization of the mixing matrix, one can express the effective neutrino mass as
(3) 
where takes the same form as in Eq. (2), is the activesterile neutrino mixing angle, and is the additional Majoranatype CPviolating phase. Using the bestfit values and from the globalfit analysis of the shortbaseline neutrino oscillation data [19, 20], one can find that the contribution from the sterile neutrino can be comparable to that from active neutrinos , which is constrained by the cosmological observations [21] and current decay experiments [22, 23, 24, 25, 26, 27].
With a tonscale target mass, the future experiments will be able to probe to the level [28], covering the whole range of in the IO case. However, in the NO case, the effective neutrino mass can be as small as when the lightest neutrino mass is vanishing, or even vanishing in the contrived region of parameter space when the cancellation among the contributions from different neutrino mass eigenstates occurs [29, 30, 31]. Moreover, the latest globalfit analysis of neutrino oscillation data [32, 33, 34] does show a preference for the NO at the confidence level (C.L.), it is worrisome that may be out of the reach of the next generation decay experiments. To quantitatively assess how likely is small, the authors of Refs. [28, 35] have carried out a Bayesian analysis and obtained the posterior distribution of , given the neutrino oscillation data, current experimental upper bounds on and the cosmological bound on the sum of three neutrino masses. For the earlier relevant works, see Refs. [36, 37, 38]. Although the impact of an eVmass sterile neutrino on the effective neutrino mass has been considered in Refs. [39, 40, 41, 42, 43, 44, 45, 46, 47], a statistical assessment is still lacking. Therefore, we are motivated to perform a Bayesian analysis of in this work by using the globalfit results of neutrino oscillation data and other available information on the absolute neutrino masses.
The rest of the present paper is organized as follows. In Section 2, we describe the necessary information for the Bayesian analysis. The prior information can be extracted from the globalfit analysis of neutrino oscillation data [33, 19], the cosmological observations [21] and the existing decay experiments [22, 23, 24, 25]. Then, the posterior distribution of the standard effective neutrino mass and that of are presented in Section 3. Twodimensional posterior probability densities in the  plane and those in the  plane have also been given, where denotes the lightest neutrino mass. Finally, we make some concluding remarks in Section 4.
2 The Bayesian Analysis
The Bayesian analysis provides us with a reasonable statistical framework to update the probability distribution of model parameters in light of the new experimental data. The posterior distribution of model parameters can be obtained according to the Bayesian theorem [48]
(4) 
where denotes the set of model parameters, stands for the available experimental data, and are the hypotheses or models with being the model index. Here is the likelihood of the data , assuming the model with the parameters , is the prior distribution of , and is the evidence. The evidence is given by
(5) 
which measures the compatibility of the model with the data, and is just the dimension of the parameter space. The hypotheses relevant for our analysis are for the NO and for the IO in the or (3+1) mixing scenario. The model parameters in the (3+1) mixing scenario include: (i) the involved neutrino oscillation parameters , where and are two masssquared differences of ordinary neutrinos; (ii) the lightest neutrino mass , which is for and for ; (iii) the Majoranatype CPviolating phases ; (iv) the phasespace factor and the nuclear matrix element for the decays. The overall likelihood function can be constructed as , and the details of the individual likelihood function are summarized as follows.

: the likelihood function of the mixing parameters . Given the function from the globalfit analysis in Ref. [33], we can fix the likelihood function , where is defined as
(6) with running over , the corresponding bestfit value from the global analysis, and the symmetrized error. See Table. 1 of Ref. [33] for more details about the globalfit results of neutrino oscillation data. To be explicit, we list the bestfit values and the corresponding symmetrized errors as below
(7) for ; and
(8) for . The latest neutrino oscillation data favor the NO over the IO at the level, i.e., the difference between the minima of in these two cases is . The preference for the NO arises mainly from two different data sets. First, the excess of like events in the multiGeV energy range in SuperKamiokande atmospheric neutrino data can be accommodated by the resonant enhancement of the oscillation probability in the channel, leading to . Second, the running longbaseline accelerator experiments T2K and NOA prefer the value of that is slightly larger than the precisely measured value from reactor neutrino experiments. Such a tension between accelerator and reactor neutrino experiments will be relieved in the NO case, contributing another to the mass ordering discrimination. To be conservative, we will take as the preference for the NO over the IO from neutrino oscillation data.

: the likelihood function for the cosmological observations on the sum of three neutrino masses . After combining several different sets of cosmological data (), the Planck Collaboration has recently updated the upper limit on the sum of neutrino masses as at the C.L. [21]. We obtain the likelihood information by making use of the Markov chain file available from the Planck Legacy Archive (PLA) ^{1}^{1}1This is based on the observations with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada.. The likelihood function of is produced and shown in Fig. 1 by marginalizing over the other cosmological parameters. Although the sampling file given by PLA has assumed a degenerate mass spectrum of neutrinos, a more solid analysis with the realistic neutrino mass spectrum should not change the result much [49]. For this reason, the likelihood shown in Fig. 1 will be used in the following discussions.

: the likelihood function derived from the experimental constraints on the effective neutrino mass or due to the existing searches for decays. For simplicity, we implement the likelihood function available from Refs. [35, 25] in our analysis. Although both and contain the information about the absolute scale of neutrino masses, the constraint on from the decays suffers from a large theoretical uncertainty in the prediction for the NME. For instance, the tightest bound comes from the KamLANDZen experiment [23], namely, . Given further uncertainties from the mixing parameters and the unknown Majorana CPviolating phases, the decays are not so informative about the absolute scale of neutrino masses when compared to the cosmological observations.

: the likelihood function encoding the globalfit analysis of sterile neutrino mass and mixing parameters . In practice, we determine the likelihood function as by using the distribution in Fig. 9 of Ref. [19]. The result of the socalled pragmatic 3+1 global fit âPrGlo17â will be utilized [19], where the tension between appearance and disappearance oscillation data can be somewhat relaxed by ignoring the excess of lowenergy like events observed in the MiniBooNE experiment.
After having the likelihood functions constructed from various experimental observations, we need to make clear the prior probability distributions of the model parameters, which reflect our knowledge about them prior to any experimental data. First, neutrino masssquared differences and mixing angles are assumed to be uniformly distributed in their allowed ranges that are wide enough to cover their global fit results. Since the oscillation data are rather informative, different choices of prior distributions of these parameters do not have much impact on the final posterior distributions. Second, the Majorana CPviolating phases are completely unknown, so it is reasonable to adopt the flat priors in the range of . In addition, we have to mention that the prior distributions for the following relevant parameters are by no means unique but will be incorporated into our calculations for practical purposes.

As indicated in Eq. (1), the phasespace factor and the NME are needed when we try to translate the experimental constraint on the halflife into that on the effective neutrino mass. The phasespace factors for different nuclear isotopes have been computed in Refs. [5, 50, 51], and we use the central values from Ref. [51], e.g., , and , which have been obtained with the axial vector coupling constant . We assume that can be described by the Gaussian distribution with the aforementioned central value and a relative error of . On the other hand, the NME for a specific nuclear isotope encoding the information about the nuclear structure has been theoretically calculated in a variety of nuclear models. The differences among these calculations can be treated as the theoretical uncertainty. We define this uncertainty as , where is the NME value of the th model, is the averaged NME value of all models, and is the total number of models. Using the tabulated NME values in Ref. [44], we find that and . Then the Gaussian distribution with the central value and the standard deviation is assumed for each nuclear isotope.

For the prior of the lightest neutrino mass , a more careful study should be performed. Four kinds of prior distributions for are usually considered: (i) a logarithmic prior on with an adjustable lower cutoff that we choose to be ; (ii) a logarithmic prior on with a natural lower cutoff at for NO or at for IO, as required by neutrino oscillation experiments; (iii) a flat prior on ; (iv) a flat prior on . The prior probability distributions have been plotted with respect to in the left panel of Fig. 2, where one can see that the flat priors on (gray solid curve) and (gray dashed curve) lead to nearly the same distribution. After incorporating the experimental limits from Planck 2018 and the decays, as shown in the right panel of Fig. 2, we observe that the logarithmic prior on (red solid curve) gives rise to a posterior distribution that is very different from those in the other scenarios. This is because a large weight has been given to very small neutrino masses in the former case. In the following discussions, we focus only on two different prior distributions, i.e., the logarithmic prior on and the logarithmic prior on , both of which are scale invariant. Since the posterior distribution of with logarithmic prior on is very similar to those with two flat priors, the posterior distribution of the effective neutrino mass in the former case should also be roughly applicable to those in the latter two cases.
Finally, we make some comments on the current experimental hint on neutrino mass ordering by combining the data sets of neutrino oscillation experiments, cosmological observations and the decays, for which the likelihood functions are given by , and , respectively. The preference odds for NO over IO can be represented by the Bayes factor, i.e., . With the help of Eq. (5), one can calculate the evidences for NO and IO and thus their ratio. The dependence of on the choice of the prior distribution is found to be very weak. Given identical prior information on both mass orderings, we consider only the cosmological observations and obtain the logarithm of the Bayes factor as ^{2}^{2}2Note that the subscript of the Bayes factor herein refers to the cosmological data that have been used in the calculations, and likewise for the Bayes factors from other data sets and their combinations., corresponding to , which is in concordance with the results from Refs. [49, 52, 53, 54]. If only the decay experiments are considered, then we get . A combination of the cosmological observations and decay data leads to . Regarding the threeflavor neutrino oscillation data, if we take the conservative choice of for two neutrino mass orderings, which has been used to construct , the logarithm of the Bayes factor turns out to be . Combining , and together, one can find the total Bayes factor . As we have mentioned before, the globalfit analysis of all the neutrino oscillation data gives rise to a preference for the NO, corresponding to . If such a stronger preference for the NO is implemented instead of the conservative one, the total Bayes factor from all the data sets becomes , showing a strong evidence for the NO according to the Jeffreys scale [55]. The addition of into the analysis does not alter the above conclusions, since the shortbaseline neutrino oscillation experiments are insensitive to the mass ordering of three ordinary neutrinos.
3 Posterior Distributions
After specifying the likelihood functions for the relevant experimental data and fixing the prior probability distributions of model parameters in the previous section, we are ready to compute the posterior distributions of the derived parameters and by using Eq. (4). In fact, the posterior probability distribution in Eq. (4) for the model parameters is calculated via the Monte Carlo sampling, which has been done with the help of the MultiNest routine [56, 57, 58]. In Fig. 3, we present the posterior sampling distributions in the  plane for the standard 3 mixing scenario (the upper row) or in the  plane for the (3+1) mixing scenario (the lower row). The scattered points stand for the sampling data, and one can read off the corresponding posterior probabilities from their colors. Now we explain how to practically do so. For a given point, one can first look at the color legend and find the value of its posterior density, which is denoted as . Then, the posterior probability can be calculated by definition as the product of and the area of a small region, in which the point is located. For instance, take a small square in the  plane, and its area is thus given by . Notice that the total posterior probability is normalized to one for each plot. Several comments on the numerical results in Fig. 3 are helpful.

In the upperleft panel, the posterior distribution in the  plane is shown for the standard 3 mixing scenario, where the logarithmic prior on is assumed. The results for the logarithmic prior on are plotted in the upperright panel. In both panels, the thin dotdashed (or dashed) curves indicate the boundaries of the effective neutrino mass in the IO (or NO) case, where the bestfit values of neutrino mixing angles and masssquared differences are input. Moreover, the current limit (taken from Ref. [23] for the tightest one) on or the future sensitivity (of a tonscale decay experiment like nEXO [52]) to is represented by three horizontal dotted lines. The wide range between the upper and lower lines can be ascribed to the NME uncertainty. Comparing the distributions in the left and right panels, one can observe that a larger weight has been given to smaller values of in the assumption of a logarithmic prior on , as already emphasized in the previous section.

An urgent question is how likely is vanishingly small in the NO case, which has been quantitatively addressed in Refs. [28, 35]. In order to draw a priorindependent conclusion from the posterior distributions, we treat the scenarios with different values of as different models. For each fixed , the posterior distribution of can be derived with the help of the likelihood . Then, one can calculate the probability for the true value of the effective neutrino mass to be above a certain . The probability contours are plotted as the blue curves in Fig. 3, where several representative values, i.e., , , and , are shown. It is evident that the probability for to be vanishingly small, e.g., , is tiny (less than ). This conclusion is independent of the priors on , as it should be. In particular, the probability for is larger than even when is located in the regime where the destructive cancellation caused by the unknown Majorana CP phases occurs.

In the two panels in the lower row of Fig. 3, the posterior probability distributions in the (3+1) mixing scenario have been presented, where the notations and conventions for the curves are the same as those in the plots in the upper row. It is straightforward to observe that the presence of the eVmass sterile neutrino shifts the effective neutrino mass to higher values. As the future tonscale decay experiments are able to explore the effective neutrino mass to the level of , the inclusion of the sterile neutrino can raise the effective mass to the level that is within the reach of the nextgeneration experiments even for a very small . If the sensitivity at the level is achieved, more than of the region of can be covered for . When , the chance for to fall into the cancellation region increases. However, even in this case, at least of the range can be probed. Therefore, in the statistical sense, it is quite promising to check the (3+1) mixing scenario with an eVmass sterile neutrino in the future decay experiments.
In Fig. 4, we present the posterior distributions in the  (the upper row) or  plane (the lower row) by marginalizing over the lightest neutrino mass instead of the Majorana CP phase . The notations and conventions are the same as those in Fig. 3. The area in the  plane is defined as in the mixing scenario, and likewise for the (3+1) mixing scenario. Now the blue solid curves in Fig. 4 stand for the contours of the probability for the effective neutrino mass to be above a certain or . These contours become dependent on the priors, because the prior information of has been integrated into the posterior distribution. It is worthwhile to notice that the dependence of posterior distributions on is very weak for the (3+1) mixing scenario. In the mixing scenario, the fine structure around due to the cancellation can be observed. Therefore, it seems difficult to determine the Majorana CP phase if takes the value far away from that in the cancellation region.
As the effective neutrino mass can be directly extracted from the experimental data on decays, it is interesting to see the posterior distribution of or , which can be obtained by marginalizing over both and . The final results can be found in Fig. 5. For the standard case in the left panel, if we choose the logarithmic prior on for NO (red solid curve), a large fraction (about ) of the probable range of is unreachable for the future tonscale decay experiments. With a logarithmic prior on (blue solid curve), the next generation experiments can cover about of the range. As we have observed before, adding an eVmass sterile neutrino can greatly enhance the probability of the effective neutrino mass to larger values. The future decay experiments with a sensitivity to the effective neutrino mass of can cover around () of the posterior space for the logarithmic prior on (the logarithmic prior on ) in the (3+1) mixing scenario. According to the posterior distributions in Fig. 5, we find that the average value of the effective neutrino mass is shifted from (or ) in the standard mixing scenario to (or ) in the (3+1) mixing scenario, with the logarithmic prior on (or on ). Therefore, a null signal from the future decay experiments will be able to set a very stringent constraint on the sterile neutrino mass and mixing angle.
4 Concluding Remarks
In this short note, we have carried out a Bayesian analysis of the effective neutrino mass in the decays in both the standard mixing scenario and the (3+1) mixing scenario. With the latest experimental information, including the globalfit analysis of neutrino oscillation data, the cosmological observations from the Planck satellite and the current limits from the decay experiments, the posterior probability distributions of the effective neutrino mass in the standard mixing scenario and in the (3+1) mixing scenario have been updated.
Our main results of the posterior distributions have been summarized in Fig. 3 and Fig. 5. Adding an eVmass sterile neutrino slightly mixing with ordinary neutrinos is likely to enhance the effective neutrino mass to the level of , which is within the reach of the next generation decay experiments, regardless of the prior information on the absolute mass scale of ordinary neutrinos. In other words, if a null signal is observed in future tonscale decay experiments, we can place very strong limits on the parameter space of the (3+1) mixing scenario, assuming that massive neutrinos are of Majorana nature. The sensitivity of future decay experiments to the sterile neutrino mass and mixing angle deserves a dedicated study, which will be left for the upcoming works.
Acknowledgements
This work was supported in part by the National Natural Science Foundation of China under grant No. 11775232 and No. 11835013, and by the CAS Center for Excellence in Particle Physics.
References
 [1] E. Majorana, “Teoria simmetrica dellâelettrone e del positrone,” Nuovo Cim. 14, 171 (1937).
 [2] G. Racah, “On the symmetry of particle and antiparticle,” Nuovo Cim. 14, 322 (1937).
 [3] M. Tanabashi et al. [Particle Data Group], “Review of Particle Physics,” Phys. Rev. D 98, no. 3, 030001 (2018).
 [4] W. H. Furry, “On transition probabilities in double betadisintegration,” Phys. Rev. 56, 1184 (1939).
 [5] W. Rodejohann, “Neutrinoless Double Beta Decay and Particle Physics,” Int. J. Mod. Phys. E 20, 1833 (2011) [arXiv:1106.1334].
 [6] S. M. Bilenky and C. Giunti, “Neutrinoless doublebeta decay: A brief review,” Mod. Phys. Lett. A 27, 1230015 (2012) [arXiv:1203.5250].
 [7] W. Rodejohann, “Neutrinoless double beta decay and neutrino physics,” J. Phys. G 39, 124008 (2012) [arXiv:1206.2560].
 [8] S. M. Bilenky and C. Giunti, “Neutrinoless DoubleBeta Decay: a Probe of Physics Beyond the Standard Model,” Int. J. Mod. Phys. A 30, no. 04n05, 1530001 (2015) [arXiv:1411.4791].
 [9] H. Päs and W. Rodejohann, “Neutrinoless Double Beta Decay,” New J. Phys. 17, no. 11, 115010 (2015) [arXiv:1507.00170].
 [10] S. Dell’Oro, S. Marcocci, M. Viel and F. Vissani, “Neutrinoless double beta decay: 2015 review,” Adv. High Energy Phys. 2016, 2162659 (2016) [arXiv:1601.07512].
 [11] C. Giunti and T. Lasserre, “eVscale Sterile Neutrinos,” arXiv:1901.08330.
 [12] A. AguilarArevalo et al. [LSND Collaboration], “Evidence for neutrino oscillations from the observation of antineutrino(electron) appearance in a antineutrino(muon) beam,” Phys. Rev. D 64, 112007 (2001) [hepex/0104049].
 [13] A. A. AguilarArevalo et al. [MiniBooNE Collaboration], “Unexplained Excess of ElectronLike Events From a 1GeV Neutrino Beam,” Phys. Rev. Lett. 102, 101802 (2009).
 [14] A. A. AguilarArevalo et al. [MiniBooNE Collaboration], “Significant Excess of ElectronLike Events in the MiniBooNE ShortBaseline Neutrino Experiment,” arXiv:1805.12028.
 [15] C. Giunti and M. Laveder, “Statistical Significance of the Gallium Anomaly,” Phys. Rev. C 83, 065504 (2011) [arXiv:1006.3244].
 [16] G. Mention, M. Fechner, T. Lasserre, T. A. Mueller, D. Lhuillier, M. Cribier and A. Letourneau, “The Reactor Antineutrino Anomaly,” Phys. Rev. D 83, 073006 (2011) [arXiv:1101.2755].
 [17] J. N. Abdurashitov et al. [SAGE Collaboration], “Measurement of the solar neutrino capture rate with gallium metal. III: Results for the 2002–2007 datataking period,” Phys. Rev. C 80, 015807 (2009) [arXiv:0901.2200].
 [18] F. Kaether, W. Hampel, G. Heusser, J. Kiko and T. Kirsten, “Reanalysis of the GALLEX solar neutrino flux and source experiments,” Phys. Lett. B 685, 47 (2010) [arXiv:1001.2731].
 [19] S. Gariazzo, C. Giunti, M. Laveder and Y. F. Li, “Updated Global 3+1 Analysis of ShortBaseLine Neutrino Oscillations,” JHEP 1706, 135 (2017) [arXiv:1703.00860].
 [20] M. Dentler, Ã. HernÃ¡ndezCabezudo, J. Kopp, P. A. N. Machado, M. Maltoni, I. MartinezSoler and T. Schwetz, “Updated Global Analysis of Neutrino Oscillations in the Presence of eVScale Sterile Neutrinos,” JHEP 1808, 010 (2018) [arXiv:1803.10661].
 [21] N. Aghanim et al. [Planck Collaboration], “Planck 2018 results. VI. Cosmological parameters,” arXiv:1807.06209.
 [22] J. B. Albert et al. [EXO200 Collaboration], “Search for Majorana neutrinos with the first two years of EXO200 data,” Nature 510, 229 (2014) [arXiv:1402.6956].
 [23] A. Gando et al. [KamLANDZen Collaboration], “Search for Majorana Neutrinos near the Inverted Mass Hierarchy Region with KamLANDZen,” Phys. Rev. Lett. 117, no. 8, 082503 (2016) Addendum: [Phys. Rev. Lett. 117, no. 10, 109903 (2016)] [arXiv:1605.02889].
 [24] M. Agostini et al., “Backgroundfree search for neutrinoless double decay of Ge with GERDA,” Nature 544, 47 (2017) [arXiv:1703.00570].
 [25] C. Alduino et al. [CUORE Collaboration], “First Results from CUORE: A Search for Lepton Number Violation via Decay of Te,” Phys. Rev. Lett. 120, no. 13, 132501 (2018) [arXiv:1710.07988].
 [26] C. E. Aalseth et al. [Majorana Collaboration], “Search for Neutrinoless DoubleÎ² Decay in Ge with the Majorana Demonstrator,” Phys. Rev. Lett. 120, no. 13, 132502 (2018) [arXiv:1710.11608].
 [27] M. Agostini et al. [GERDA Collaboration], “Improved Limit on Neutrinoless Double Decay of Ge from GERDA Phase II,” Phys. Rev. Lett. 120, no. 13, 132503 (2018) [arXiv:1803.11100].
 [28] M. Agostini, G. Benato and J. Detwiler, “Discovery probability of nextgeneration neutrinoless double Î² decay experiments,” Phys. Rev. D 96, no. 5, 053001 (2017) [arXiv:1705.02996].
 [29] Z. z. Xing, “Vanishing effective mass of the neutrinoless double beta decay?,” Phys. Rev. D 68, 053002 (2003) [hepph/0305195].
 [30] Z. z. Xing, Z. h. Zhao and Y. L. Zhou, “How to interpret a discovery or null result of the decay,” Eur. Phys. J. C 75, no. 9, 423 (2015) [arXiv:1504.05820].
 [31] Z. z. Xing and Z. h. Zhao, “The effective neutrino mass of neutrinoless doublebeta decays: how possible to fall into a well,” Eur. Phys. J. C 77, no. 3, 192 (2017) [arXiv:1612.08538].
 [32] P. F. de Salas, D. V. Forero, C. A. Ternes, M. Tortola and J. W. F. Valle, “Status of neutrino oscillations 2018: 3 hint for normal mass ordering and improved CP sensitivity,” Phys. Lett. B 782, 633 (2018) [arXiv:1708.01186].
 [33] F. Capozzi, E. Lisi, A. Marrone and A. Palazzo, “Current unknowns in the three neutrino framework,” Prog. Part. Nucl. Phys. 102, 48 (2018) [arXiv:1804.09678].
 [34] I. Esteban, M. C. GonzalezGarcia, A. HernandezCabezudo, M. Maltoni and T. Schwetz, “Global analysis of threeflavour neutrino oscillations: synergies and tensions in the determination of , and the mass ordering,” arXiv:1811.05487.
 [35] A. Caldwell, A. Merle, O. Schulz and M. Totzauer, “Global Bayesian analysis of neutrino mass data,” Phys. Rev. D 96, no. 7, 073001 (2017) [arXiv:1705.01945].
 [36] G. Benato, “Effective Majorana Mass and Neutrinoless Double Beta Decay,” Eur. Phys. J. C 75, no. 11, 563 (2015) [arXiv:1510.01089].
 [37] J. Zhang and S. Zhou, “Determination of neutrino mass ordering in future Gebased neutrinoless doublebeta decay experiments,” Phys. Rev. D 93, no. 1, 016008 (2016) [arXiv:1508.05472].
 [38] S. F. Ge and M. Lindner, “Extracting Majorana properties from strong bounds on neutrinoless double beta decay,” Phys. Rev. D 95, no. 3, 033003 (2017) [arXiv:1608.01618].
 [39] S. Goswami and W. Rodejohann, “Constraining mass spectra with sterile neutrinos from neutrinoless double beta decay, tritium beta decay and cosmology,” Phys. Rev. D 73, 113003 (2006) [hepph/0512234].
 [40] S. Goswami and W. Rodejohann, “MiniBooNE results and neutrino schemes with 2 sterile neutrinos: Possible mass orderings and observables related to neutrino masses,” JHEP 0710, 073 (2007) [arXiv:0706.1462].
 [41] J. Barry, W. Rodejohann and H. Zhang, “Light Sterile Neutrinos: Models and Phenomenology,” JHEP 1107, 091 (2011) [arXiv:1105.3911].
 [42] Y. F. Li and S. s. Liu, “Vanishing effective mass of the neutrinoless double beta decay including light sterile neutrinos,” Phys. Lett. B 706, 406 (2012) [arXiv:1110.5795].
 [43] I. Girardi, A. Meroni and S. T. Petcov, “Neutrinoless Double Beta Decay in the Presence of Light Sterile Neutrinos,” JHEP 1311, 146 (2013) [arXiv:1308.5802].
 [44] P. Guzowski, L. Barnes, J. Evans, G. Karagiorgi, N. McCabe and S. SoldnerRembold, “Combined limit on the neutrino mass from neutrinoless doubleÎ² decay and constraints on sterile Majorana neutrinos,” Phys. Rev. D 92, no. 1, 012002 (2015) [arXiv:1504.03600].
 [45] C. Giunti and E. M. Zavanin, “Predictions for Neutrinoless DoubleBeta Decay in the 3+1 Sterile Neutrino Scenario,” JHEP 1507, 171 (2015) [arXiv:1505.00978].
 [46] S. F. Ge, W. Rodejohann and K. Zuber, “Halflife Expectations for Neutrinoless Double Beta Decay in Standard and NonStandard Scenarios,” Phys. Rev. D 96, no. 5, 055019 (2017) [arXiv:1707.07904].
 [47] J. H. Liu and S. Zhou, “Another look at the impact of an eVmass sterile neutrino on the effective neutrino mass of neutrinoless doublebeta decays,” Int. J. Mod. Phys. A 33, no. 02, 1850014 (2018) [arXiv:1710.10359].
 [48] D. Silva and J. Skilling, “Data Analysis: A Bayesian Tutorial,” Oxford University Press, Oxford, 2006.
 [49] S. Hannestad and T. Schwetz, “Cosmology and the neutrino mass ordering,” JCAP 1611, no. 11, 035 (2016) [arXiv:1606.04691].
 [50] J. Suhonen and O. Civitarese, “Weakinteraction and nuclearstructure aspects of nuclear double beta decay,” Phys. Rept. 300, 123 (1998).
 [51] J. Kotila and F. Iachello, “Phase space factors for double decay,” Phys. Rev. C 85, 034316 (2012) [arXiv:1209.5722].
 [52] M. Gerbino, M. Lattanzi, O. Mena and K. Freese, “A novel approach to quantifying the sensitivity of current and future cosmological datasets to the neutrino mass ordering through Bayesian hierarchical modeling,” Phys. Lett. B 775, 239 (2017) [arXiv:1611.07847].
 [53] S. Vagnozzi, E. Giusarma, O. Mena, K. Freese, M. Gerbino, S. Ho and M. Lattanzi, “Unveiling secrets with cosmological data: neutrino masses and mass hierarchy,” Phys. Rev. D 96, no. 12, 123503 (2017) [arXiv:1701.08172].
 [54] F. Capozzi, E. Di Valentino, E. Lisi, A. Marrone, A. Melchiorri and A. Palazzo, “Global constraints on absolute neutrino masses and their ordering,” Phys. Rev. D 95, no. 9, 096014 (2017) [arXiv:1703.04471].
 [55] R. Trotta, “Bayes in the sky: Bayesian inference and model selection in cosmology,” Contemp. Phys. 49, 71 (2008) [arXiv:0803.4089].
 [56] F. Feroz and M. P. Hobson, “Multimodal nested sampling: an efficient and robust alternative to MCMC methods for astronomical data analysis,” Mon. Not. Roy. Astron. Soc. 384, 449 (2008) [arXiv:0704.3704].
 [57] F. Feroz, M. P. Hobson and M. Bridges, “MultiNest: an efficient and robust Bayesian inference tool for cosmology and particle physics,” Mon. Not. Roy. Astron. Soc. 398, 1601 (2009) [arXiv:0809.3437].
 [58] F. Feroz, M. P. Hobson, E. Cameron and A. N. Pettitt, “Importance Nested Sampling and the MultiNest Algorithm,” arXiv:1306.2144.