Impact of an eV-mass sterile neutrino on the neutrinoless double-beta decays: a Bayesian analysis

To quantitatively assess the impact of an eV-mass sterile neutrino on the neutrinoless double-beta ($0\nu \beta \beta$) decays, we calculate the posterior probability distribution of the relevant effective neutrino mass $|m^\prime_{ee}|$ in the (3+1)$\nu$ mixing scenario, following the Bayesian statistical approach. The latest global-fit analysis of neutrino oscillation data, the cosmological bound on the sum of three active neutrino masses from {\it Planck}, and the constraints from current $0\nu\beta\beta$ 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 $\overline{|m^{}_{ee}|} = 3.37\times 10^{-3}~{\rm eV}$ (or $7.71\times 10^{-3}~{\rm eV}$) in the standard 3$\nu$ mixing scenario to $\overline{|m^{\prime}_{ee}|}=2.54\times 10^{-2}~{\rm eV}$ (or $2.56\times 10^{-2}~{\rm eV}$) in the (3+1)$\nu$ 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 $0\nu\beta\beta$ decay experiment with a sensitivity to $|m^{}_{ee}| \approx \mathcal{O}(10^{-2}_{})~{\rm eV}$ will be able to set a very stringent constraint on the sterile neutrino mass and the active-sterile mixing angle.


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 double-beta (0νββ) 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 0νββ decays A(Z, N ) → A(Z + 2, N − 2) + 2e − could take place in some even-even nuclei, namely, both the proton number Z and the neutron number N for the nuclear isotope A(Z, N ) are even. Assuming the exchange of light Majorana neutrinos to be responsible for the 0νββ decays, one can find that the half-life of the relevant nuclear isotope is given by [5] (T 0ν where G 0ν is the phase-space factor, M 0ν is the nuclear matrix element (NME), and m e is the electron mass. In Eq. (1), the effective neutrino mass |m ee | collects the contributions from light Majorana neutrinos involved in the 0νββ decays.
In the standard three neutrino (3ν) mixing scenario, the effective neutrino mass is defined as |m ee | ≡ |m 1 U 2 e1 + m 2 U 2 e2 + m 3 U 2 e3 |, where the absolute neutrino masses m i and the lepton flavor mixing matrix elements U ei (for i = 1, 2, 3) appear. When the conventional parametrization of the flavor mixing matrix U is adopted [3], i.e., U e1 = cos θ 13 cos θ 12 e iρ/2 , U e2 = cos θ 13 sin θ 12 and U e3 = sin θ 13 e iσ/2 , we have m ee ≡ m 1 cos 2 θ 13 cos 2 θ 12 e iρ + m 2 cos 2 θ 13 sin 2 θ 12 + m 3 sin 2 θ 13 e iσ , where {θ 12 , θ 13 } are two of three neutrino mixing angles, and {ρ, σ} are the Majorana-type CPviolating phases. Note that m 2 is nonzero no matter whether the normal neutrino mass ordering (NO) with m 1 < m 2 < m 3 or the inverted neutrino mass ordering (IO) with m 3 < m 1 < m 2 is considered. Therefore, such a parametrization is favorable in the discussions about the limiting case of m 1 → 0 (for NO) or m 3 → 0 (for IO), for which one of two Majorana-type CP violating phases just disappears together with the lightest neutrino mass. However, if the eV-mass sterile neutrino indeed exists as a solution to the anomalies in the short-baseline neutrino experiments [11][12][13][14][15][16][17][18], it will contribute as well to the 0νββ decays. In this case, the effective neutrino mass is given by |m ee | ≡ |m 1 | with m 4 being the mass of the sterile neutrino and V ei (for i = 1, 2, 3, 4) being the first-row 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 where m ee takes the same form as in Eq. (2), θ 14 is the active-sterile neutrino mixing angle, and ω is the additional Majorana-type CP-violating phase. Using the best-fit values ∆m 2 41 ≡ m 2 4 − m 2 1 = 1.7 eV 2 and sin 2 θ 14 = 0.019 from the global-fit analysis of the short-baseline neutrino oscillation data [19,20], one can find that the contribution from the sterile neutrino |m 4 sin 2 θ 14 | ≈ 2.5×10 −2 eV can be comparable to that from active neutrinos |m ee | 0.1 eV, which is constrained by the cosmological observations [21] and current 0νββ decay experiments [22][23][24][25][26][27].
With a ton-scale target mass, the future 0νββ experiments will be able to probe |m ee | to the O(10 −2 ) eV level [28], covering the whole range of |m ee | in the IO case. However, in the NO case, the effective neutrino mass can be as small as |m ee | ≈ (1.6 · · · 3.6) × 10 −3 eV when the lightest neutrino mass m 1 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 global-fit analysis of neutrino oscillation data [32][33][34] does show a preference for the NO at the 3σ confidence level (C.L.), it is worrisome that |m ee | may be out of the reach of the next generation 0νββ decay experiments. To quantitatively assess how likely |m ee | is small, the authors of Refs. [28,35] have carried out a Bayesian analysis and obtained the posterior distribution of |m ee |, given the neutrino oscillation data, current experimental upper bounds on |m ee | 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 eV-mass sterile neutrino on the effective neutrino mass |m ee | 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 |m ee | in this work by using the global-fit 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 global-fit analysis of neutrino oscillation data [19,33], the cosmological observations [21] and the existing 0νββ decay experiments [22][23][24][25]. Then, the posterior distribution of the standard effective neutrino mass |m ee | and that of |m ee | are presented in Section 3. Two-dimensional posterior probability densities in the |m ee |-m L plane and those in the |m ee |-ρ plane have also been given, where m L denotes the lightest neutrino mass. Finally, we make some concluding remarks in Section 4.

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] where Θ denotes the set of model parameters, D stands for the available experimental data, and {H i } are the hypotheses or models with i being the model index. Here L(D|Θ, H i ) is the likelihood of the data D, assuming the model H i with the parameters Θ, π(Θ, H i ) is the prior distribution of Θ, and Z i is the evidence. The evidence Z i is given by which measures the compatibility of the model with the data, and N is just the dimension of the parameter space. The hypotheses relevant for our analysis are H NO for the NO and H IO for the IO in the 3ν or (3+1)ν mixing scenario. The model parameters in the (3+1)ν mixing scenario include: (i) the involved neutrino oscillation parameters {sin 2 θ 13 , sin 2 θ 12 , sin 2 θ 14 , ∆m 2 sol , ∆m 2 atm , ∆m 2 41 }, where ∆m 2 sol ≡ m 2 2 − m 2 1 and ∆m 2 atm ≡ m 2 3 − (m 2 2 + m 2 1 )/2 are two mass-squared differences of ordinary neutrinos; (ii) the lightest neutrino mass m L , which is m 1 for H NO and m 3 for H IO ; (iii) the Majorana-type CP-violating phases {ρ, σ, ω}; (iv) the phase-space factor and the nuclear matrix element {G 0ν , |M 0ν |} for the 0νββ decays. The overall likelihood function can be constructed as L = L 3ν × L cosmo × L 0νββ × L sterile , and the details of the individual likelihood function are summarized as follows.
• L 3ν : the likelihood function of the 3ν mixing parameters {sin 2 θ 13 , sin 2 θ 12 , ∆m 2 sol , ∆m 2 atm }. Given the ∆χ 2 function from the global-fit analysis in Ref. [33], we can fix the likelihood function L 3ν = exp(−∆χ 2 /2), where ∆χ 2 is defined as with Θ i running over {sin 2 θ 13 , sin 2 θ 12 , ∆m 2 sol , ∆m 2 atm }, Θ bf i the corresponding best-fit value from the global analysis, and σ i the symmetrized 1σ error. See Table. 1 of Ref. [33] for more details about the global-fit results of neutrino oscillation data. To be explicit, we list the best-fit values and the corresponding symmetrized 1σ errors as below for H NO ; and sin 2 θ 12 = (3.03 ± 0.14) × 10 −1 , ∆m 2 sol = (7.34 ± 0.16) × 10 −5 eV 2 , for H IO . The latest neutrino oscillation data favor the NO over the IO at the 3σ level, i.e., the difference between the minima of χ 2 in these two cases is ∆χ 2 min ≡ χ IO min − χ NO min ≈ 9. The preference for the NO arises mainly from two different data sets. First, the excess of ν e -like events in the multi-GeV energy range in Super-Kamiokande atmospheric neutrino data can be accommodated by the resonant enhancement of the oscillation probability in the ν µ → ν e channel, leading to ∆χ 2 min ≈ 4. Second, the running long-baseline accelerator experiments T2K and NOνA prefer the value of θ 13 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 ∆χ 2 min ≈ 4 to the mass ordering discrimination. To be conservative, we will take ∆χ 2 min = 4 as the preference for the NO over the IO from neutrino oscillation data.
• L cosmo : the likelihood function for the cosmological observations on the sum of three neutrino masses Σ ≡ m 1 + m 2 + m 3 . After combining several different sets of cosmological data (P lanck TT, TE, EE+lowE+lensing+BAO), the Planck Collaboration has recently updated the upper limit on the sum of neutrino masses as Σ < 0.12 eV at the 95% C.L. [21]. We obtain the likelihood information by making use of the Markov chain file available from the The likelihood function L cosmo for the sum of three neutrino masses Σ ≡ m 1 + m 2 + m 3 from cosmological observations, which has been derived by combining the P lanck TT, TE, EE + lowE + lensing + BAO data sets [21].
Planck Legacy Archive (PLA) 1 . 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.
• L 0νββ : the likelihood function derived from the experimental constraints on the effective neutrino mass |m ee | or |m ee | due to the existing searches for 0νββ decays. For simplicity, we implement the likelihood function available from Refs. [25,35] in our analysis. Although both L 0νββ and L cosmo contain the information about the absolute scale of neutrino masses, the constraint on |m ee | from the 0νββ decays suffers from a large theoretical uncertainty in the prediction for the NME. For instance, the tightest bound comes from the KamLAND-Zen experiment [23], namely, |m ee | (61 · · · 165) meV. Given further uncertainties from the mixing parameters and the unknown Majorana CP-violating phases, the 0νββ decays are not so informative about the absolute scale of neutrino masses when compared to the cosmological observations.
• L sterile : the likelihood function encoding the global-fit analysis of sterile neutrino mass and mixing parameters {θ 14 , ∆m 2 41 }. In practice, we determine the likelihood function as L sterile = exp[−∆χ 2 sterile (θ 14 , ∆m 2 41 )/2] by using the ∆χ 2 distribution in Fig. 9 of Ref. [19]. The result of the so-called 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 low-energy ν e -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 mass-squared differences and mixing angles {sin 2 θ 13 , sin 2 θ 12 , sin 2 θ 14 , ∆m 2 sol , ∆m 2 atm , ∆m 2 41 } 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 CP-violating phases are completely unknown, so it is reasonable to adopt the flat priors in the range of [0 · · · 2π). 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 phase-space factor G 0ν and the NME |M 0ν | are needed when we try to translate the experimental constraint on the half-life into that on the effective neutrino mass. The phase-space factors for different nuclear isotopes have been computed in Refs. [5,50,51], and we use the central values from Ref. [51], e.g., G 0ν ( 76 Ge) = 6.15 × 10 −15 yr −1 , G 0ν ( 130 Te) = 3.70 × 10 −14 yr −1 and G 0ν ( 136 Xe) = 3.79 × 10 −14 yr −1 , which have been obtained with the axial vector coupling constant g A = 1.27. We assume that G 0ν can be described by the Gaussian distribution with the aforementioned central value and a relative error of 7%. 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 σ NME ≡ i (|M i 0ν | − |M 0ν |) 2 /n NME , where |M i 0ν | is the NME value of the ith model, |M 0ν | is the averaged NME value of all models, and n NME is the total number of models. Using the tabulated NME values in Ref. [44], we find that |M 0ν |( 76 Ge, 130 Te, 136 Xe) = (4.88, 3.94, 2.73) and σ NME ( 76 Ge, 130 Te, 136 Xe) = (1.14, 0.90, 0.80). Then the Gaussian distribution with the central value |M 0ν | and the standard deviation σ NME is assumed for each nuclear isotope.
• For the prior of the lightest neutrino mass m L , a more careful study should be performed.
Four kinds of prior distributions for m L are usually considered: (i) a logarithmic prior on m L with an adjustable lower cutoff that we choose to be 10 −4 eV; (ii) a logarithmic prior on Σ with a natural lower cutoff at 0.06 eV for NO or at 0.1 eV for IO, as required by neutrino oscillation experiments; (iii) a flat prior on m L ; (iv) a flat prior on Σ. The prior probability distributions have been plotted with respect to log 10 (m L /eV) in the left panel of Fig. 2, where one can see that the flat priors on m L (gray solid curve) and Σ (gray dashed curve) lead to nearly the same distribution. After incorporating the experimental limits from Planck 2018 and the 0νββ decays, as shown in the right panel of Fig. 2, we observe that the logarithmic prior on m L (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 m L and the logarithmic prior on Σ, both of which are scale invariant. Since the posterior distribution of m L 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 0νββ decays, for which the likelihood functions are given by L 3ν , L cosmo and L 0νββ , respectively. The preference odds for NO over IO can be represented by the Bayes factor, i.e., B ≡ Z NO /Z IO . With the help of Eq. (5), one can calculate the evidences for NO and IO and thus their ratio. The dependence of B on the choice of the m L prior distribution is found to be very weak. Given identical prior information on both mass orderings, we consider only the cosmological observations L cosmo and obtain the logarithm of the Bayes factor as log(B cosmo ) ≈ 0.85 2 , corresponding to B cosmo ≈ 2.34, which is in concordance with the results from Refs. [49,[52][53][54]. If only the 0νββ decay experiments are considered, then we get log(B 0νββ ) ≈ 0.2. A combination of the cosmological observations and 0νββ decay data leads to log(B cosmo+0νββ ) ≈ 1.1. Regarding the three-flavor neutrino oscillation data, if we take the conservative choice of ∆χ 2 min ≈ 4 for two neutrino mass orderings, which has been used to construct L 3ν , the logarithm of the Bayes factor turns out to be log(B 3ν ) = 2. Combining L cosmo , L 0νββ and L 3ν together, one can find the total Bayes factor B tot ≈ 22. As we have mentioned before, the global-fit analysis of all the neutrino oscillation data gives rise to a 3σ preference for the NO, corresponding to B 3ν ≈ 90. 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 B tot ≈ 270, showing a strong evidence for the NO according to the Jeffreys scale [55]. The addition of L sterile into the analysis does not alter the above conclusions, since the short-baseline neutrino oscillation experiments are insensitive to the mass ordering of three ordinary neutrinos.

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 |m ee | and |m ee | 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 |m ee |-m L plane for the standard 3ν mixing scenario (the upper row) or in the |m ee |-m L 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 p. Then, the posterior probability P can be calculated by definition as the product of p and the area A of a small region, in which the point is located. For instance, take a small square in the |m ee |-m L plane, and its area is thus given by dA ≡ d [log 10 (|m ee |/eV)] × d [log 10 (m L /eV)]. Notice that the total posterior probability is normalized to one for each plot. Several comments on the numerical results in Fig. 3 are helpful.
1. In the upper-left panel, the posterior distribution in the |m ee |-m L plane is shown for the standard 3ν mixing scenario, where the logarithmic prior on m L is assumed. The results for the logarithmic prior on Σ are plotted in the upper-right panel. In both panels, the thin dot-dashed (or dashed) curves indicate the boundaries of the effective neutrino mass |m ee | in the IO (or NO) case, where the best-fit values of neutrino mixing angles and mass-squared differences are input. Moreover, the current limit (taken from Ref. [23] for the tightest one) on or the future sensitivity (of a ton-scale 0νββ decay experiment like nEXO [52]) to |m ee | 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 m L in the assumption of a logarithmic prior on m L , as already emphasized in the previous section.
2. An urgent question is how likely |m ee | is vanishingly small in the NO case, which has been quantitatively addressed in Refs. [28,35]. In order to draw a prior-independent conclusion from the posterior distributions, we treat the scenarios with different values of m L as different models. For each fixed m L , the posterior distribution of |m ee | can be derived with the help of the likelihood L 3ν . Then, one can calculate the probability for the true value of the effective neutrino mass to be above a certain |m ee |. The probability contours are plotted as the blue curves in Fig. 3, where several representative values, i.e., 68%, 95%, 99% and 99.7%, are shown. It is evident that the probability for |m ee | to be vanishingly small, e.g., |m ee | < 10 −4 eV, is tiny (less than 0.3%). This conclusion is independent of the priors on m L , as it should be. In particular, the probability for |m ee | > 10 −3 eV is larger than 95% even when m L is located in the regime where the destructive cancellation caused by the unknown Majorana CP phases occurs.
3. 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 eV-mass sterile neutrino shifts the effective neutrino mass to higher values. As the future ton-scale 0νββ decay experiments are able to explore the effective neutrino mass to the level of O(10 −2 ) eV, the inclusion of the sterile neutrino can raise the effective mass to the level that is within the reach of the next-generation experiments even for a very small m L . If the sensitivity at the O(10 −2 ) eV level is achieved, more than 99.7% of the region of |m ee | can be covered for m L 10 −2 eV. When m L 10 −2 eV, the chance for |m ee | to fall into the cancellation region increases. However, even in this case, at least 95% of the |m ee | range can be probed. Therefore, in the statistical sense, it is quite promising to check the (3+1)ν mixing scenario with an eV-mass sterile neutrino in the future 0νββ decay experiments.
In Fig. 4, we present the posterior distributions in the |m ee |-ρ (the upper row) or |m ee |-ρ plane (the lower row) by marginalizing over the lightest neutrino mass m L instead of the Majorana CP phase ρ. The notations and conventions are the same as those in Fig. 3. The area in the |m ee |-ρ plane is defined as dA ≡ d [log 10 (|m ee |/eV)] × d [ρ/rad] in the 3ν 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 |m ee | or |m ee |. These contours become dependent on the m L priors, because the prior information of m L 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 3ν mixing scenario, the fine structure around ρ ≈ π due to the cancellation can be observed. Therefore, it seems difficult to determine the Majorana CP phase ρ if |m ee | 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 0νββ decays, it is interesting to see the posterior distribution of |m ee | or |m ee |, which can be obtained by marginalizing over both m L and ρ. The final results can be found in Fig. 5. For the standard 3ν case in the left panel, if we choose the logarithmic prior on m L for NO (red solid curve), a large fraction (about 92%) of the probable range of |m ee | is unreachable for the future ton-scale 0νββ decay experiments. With a logarithmic prior on Σ (blue solid curve), the next generation experiments can cover about 41% of the range. As we have observed before, adding an eV-mass sterile neutrino can greatly enhance the probability of the effective neutrino mass |m ee | to larger values. The future 0νββ decay experiments with a sensitivity to the effective neutrino mass of O(10 −2 ) eV can cover around 99.4% (97.4%) of the posterior space for the logarithmic prior on m L (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 |m ee | = 3.37 × 10 −3 eV (or 7.71 × 10 −3 eV) in the standard 3ν mixing scenario to |m ee | = 2.54 × 10 −2 eV (or 2.56 × 10 −2 eV) in the (3+1)ν mixing scenario, with the logarithmic prior on m L (or on Σ). Therefore, a null signal from the future 0νββ decay experiments will be able to set a very stringent constraint on the sterile neutrino mass and mixing angle.

Concluding Remarks
In this short note, we have carried out a Bayesian analysis of the effective neutrino mass in the 0νββ decays in both the standard 3ν mixing scenario and the (3+1)ν mixing scenario. With the latest experimental information, including the global-fit analysis of neutrino oscillation data, the cosmological observations from the Planck satellite and the current limits from the 0νββ decay experiments, the posterior probability distributions of the effective neutrino mass |m ee | in the standard 3ν mixing scenario and |m ee | 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 eV-mass sterile neutrino slightly mixing with ordinary neutrinos is likely to enhance the effective neutrino mass to the level of O(10 −2 ) eV, which is within the reach of the next generation 0νββ 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 ton-scale 0νββ 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 0νββ decay experiments to the sterile neutrino mass and mixing angle deserves a dedicated study, which will be left for the upcoming works.