Tentative sensitivity of future 0νββ-decay experiments to neutrino masses and Majorana CP phases

In the near future, the neutrinoless double-beta (0νββ) decay experiments will hopefully reach the sensitivity of a few meV to the effective neutrino mass |mββ|. In this paper, we tentatively examine the sensitivity of future 0νββ-decay experiments to neutrino masses and Majorana CP phases by following the Bayesian statistical approach. Provided experimental setups corresponding to the experimental sensitivity of |mββ| ≃ 1 meV, the null observation of 0νββ decays in the case of normal neutrino mass ordering leads to a very competitive bound on the lightest neutrino mass m1. Namely, the 95% credible interval in the Bayesian approach turns out to be 1.6 meV ≲ m1 ≲ 7.3 meV or 0.3 meV ≲ m1 ≲ 5.6 meV when the uniform prior on m1/eV or on log10(m1/eV) is adopted. Moreover, one of two Majorana CP phases is strictly constrained, i.e., 140° ≲ ρ ≲ 220° for both scenarios of prior distributions of m1. In contrast, if a relatively worse experimental sensitivity of |mββ| ≃ 10 meV is assumed, the constraint on the lightest neutrino mass becomes accordingly 0.6 meV ≲ m1 ≲ 26 meV or 0 ≲ m1 ≲ 6.1 meV, while two Majorana CP phases will be essentially unconstrained. In the same statistical framework, the prospects for the determination of neutrino mass ordering and the discrimination between Majorana and Dirac nature of massive neutrinos in the 0νββ-decay experiments are also discussed. Given the experimental sensitivity of |mββ| ≃ 10 meV (or 1 meV), the strength of evidence to exclude the Majorana nature under the null observation of 0νββ decays is found to be inconclusive (or strong), no matter which of two priors on m1 is taken.


Introduction
The experimental observation of neutrinoless double-beta (0νββ) decays A Z N → A Z+2 N +2e − of some heavy nuclei A Z N , which possess an even atomic number Z and an even mass number A, is currently the most promising way to probe the Majorana nature of massive neutrinos and to prove the existence of lepton number violation in nature [1]. In the framework of three-flavor neutrino mixing, the 0νββ decays are mediated by three active neutrinos and the corresponding half-life of the 0νββ-decaying even-even nuclear isotope is given by [2] T 0ν where G 0ν denotes the relevant phase-space factor, M 0ν is the nuclear matrix element (NME), and m e = 0.511 MeV is the electron mass. As advocated by Particle Data Group [3], the lepton flavor mixing matrix U is usually parametrized in terms of three mixing angles {θ 12 , θ 13 , θ 23  where c ij ≡ cos θ ij and s ij ≡ sin θ ij (for ij = 12, 13, 23) and P ν = diag{1, e iφ 21 /2 , e iφ 31 /2 }. The effective neutrino mass |m ββ | for 0νββ decays appearing in eq. (1.1) reads |m ββ | ≡ 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σ , (1.3) where m i (for i = 1, 2, 3) stand for the absolute masses of three ordinary neutrinos. Out of three neutrino mixing angles only two {θ 12 , θ 13  2δ have been redefined as in ref. [4]. The other neutrino mixing angle θ 23 is irrelevant for 0νββ decays.
In the past few decades, neutrino oscillation experiments have measured with a good precision the two neutrino mixing angles {θ 12 , θ 13 }, and two independent neutrino masssquared differences ∆m 2 21 ≡ m 2 2 − m 2 1 and |∆m 2 31 | ≡ |m 2 3 − m 2 1 | [5]. In the near future, the JUNO experiment will offer an unambiguous answer to whether neutrino mass ordering is normal m 1 < m 2 < m 3 (NO) or inverted m 3 < m 1 < m 2 (IO), and improve the precisions of all four parameters {sin 2 θ 12 , sin 2 θ 13 } and {∆m 2 21 , ∆m 2 31 } to the O(0.1%) level [6][7][8][9]. Given the precision data on these parameters, the observation of 0νββ decays will be extremely important in the determination of other fundamental parameters that cannot be probed in neutrino oscillation experiments, such as the absolute scale m L of neutrino masses, i.e., the lightest neutrino mass m 1 (for NO) or m 3 (for IO) and two Majorana CP phases {ρ, σ}. In particular, the experimental constraints on the Majorana CP phases can be obtained only in the lepton-number-violating processes, among which 0νββ decays should be most feasible and promising.
In this paper, we demonstrate that it is scientifically beneficial and even indispensable to reach the meV frontier of |m ββ |, by quantitatively examining the projected sensitivities of future 0νββ-decay experiments to the absolute neutrino masses and two Majorana CP phases. The main motivation for such an investigation is three-fold. First, the upper bound on the absolute scale of neutrino masses extracted from the tritium betadecay experiments is m β < 2.3 eV at the 95% confidence level (CL) from Mainz [10], m β < 2.2 eV at the 95% CL from Troitsk [11], and m β < 1.1 eV at the 90% CL from KATRIN [12], where the effective neutrino mass m β for beta decays is defined as m β ≡ m 2 1 |U e1 | 2 + m 2 2 |U e2 | 2 + m 2 3 |U e3 | 2 1/2 with the moduli of the matrix elements of lepton flavor mixing matrix being |U e1 | = cos θ 13 cos θ 12 , |U e2 | = cos θ 13 sin θ 12 and |U e3 | = sin θ 13 in the standard parametrization. The future operation of KATRIN [13,14] and the next-generation tritium beta-decay experiment Project 8 [15] will hopefully be able to bring the upper limit down to m β 0.2 eV and m β 40 meV, respectively. On the other hand, the cosmological observations of cosmic microwave background by the Planck satellite gives the most restrictive bound on the sum of three neutrino masses Σ ≡ m 1 + m 2 + m 3 < 0.12 eV [16]. However, there is still a long way to go until the mass region of a few meV can be accessed. Second, if massive neutrinos are indeed Majorana particles, then two associated CP-violating phases {ρ, σ} are fundamental parameters in nature and must be experimentally determined. At present, the 0νββ decays are the unique feasible pathway to get close to this goal [17][18][19][20][21][22][23]. In this connection, the neutrino-antineutrino oscillations and other lepton-number-violating processes could in principle also provide some useful information about Majorana CP phases [24,25], but the observations of these processes are currently still far away from reality. Even though a number of analytical studies of the effective neutrino mass |m ββ | have been performed in the literature [17][18][19][20][21][22][23], some particular values of |m ββ | are assumed to derive the constraints on neutrino masses and Majorana CP phases. However, the effective neutrino mass |m ββ | itself is not a direct observable of 0νββ-decay experiments. A robust statistical analysis is desirable to answer JHEP03(2021)084 the following question: (i) given an experimental setup, what can we learn from a null signal after systematically taking into account the uncertainties of oscillation data, the nuclear matrix element and the phase-space factor? (ii) or conversely, to derive competitive bounds on the neutrino mass and Majorana phases, which kind of experimental setups will be required in the future? Finally, the latest global-fit analysis of neutrino oscillation data yields a 2σ hint at the normal neutrino mass ordering [26], so it is timely to investigate the physics potential of the future 0νββ-decay experiments that aim at the ultimate discovery even in the NO case. Strategically speaking, whether the target value of the effective neutrino mass |m ββ | should be set to 10 meV or 1 meV makes a significant difference.
The remaining part of this paper is structured as follows. In section 2, the sensitivities of 0νββ-decay experiments to the half-life T 0ν 1/2 and to the effective neutrino mass |m ββ | are discussed. In section 3, the sensitivities to the absolute neutrino mass scale and Majorana CP phases are examined by following the Bayesian statistics, where the physics potential of future experiments is investigated. Then we implement the Bayesian factors to discriminate between NO and IO, as well as the Dirac and Majorana nature of the massive neutrinos, in a quantitative way. Finally, we summarize our main conclusions in section 4.

Sensitivities to T 0ν 1/and |m ββ |
A number of nuclear isotopes have been found to be suitable for observing 0νββ decays [1,2]. In the present work we take the nuclear isotope 136 Xe for illustration, which has been implemented in the currently leading 0νββ-decay experiments (e.g., KamLAND-Zen [27] and EXO-200 [28,29]), and the other candidates can be studied in a similar way. It should be helpful to first establish the relation between an experimental configuration and its sensitivity to the effective neutrino mass |m ββ |. This task has already been accomplished in ref. [30], but we shall reproduce the main results in this section for completeness and for setting up our notations for later discussions. As is well known, for a given setup of the 0νββ-decay experiment, its sensitivity to the half-life T 0ν 1/2 for 0νββ decays can be derived by using the following formula [30] T 0ν 1/2 = ln 2 · where N A = 6.022 × 10 23 is the Avogadro's constant, m iso is the molar mass of the relevant nuclear isotope, ξ ≡ M iso ·t is the total exposure of the experiment with M iso being the total target mass of the decaying isotope and t being the running time of the experiment, and is the detection efficiency of the signal event. In addition, S (B) in eq. (2.1) is defined as the expected number of signal events within the region of interest (ROI) when a specified fraction q of a set of identical experiments can report a discovery of the 0νββ decay signal at the CL ≥ p [31], where the dependence on the total number of background events B ≡ b · ξ has been explicitly stated with b being the background index (in units of counts per ton·yr).
Given the expectation value µ of the total number of events, the number of counts n truly observed in the experiment statistically fluctuates according to the Poisson distribution, for which the probability distribution function (PDF) is given by PDF Poisson (n, µ) = JHEP03(2021)084 |M 0ν | 4.20 as compiled in ref. [32] have been considered. The sensitivity to T 0ν 1/2 and |m ββ | at the 3σ CL for 90% experiments has been plotted as contours in the plane of the effective exposure ξ and the effective background index b (right panel). The colored bands in both left and right panels stem from the NME uncertainty. The corresponding dashed curves in the middle of the band of the right panel is obtained by using the average NME value of |M 0ν | = 2.94. e −µ µ n /n! and the corresponding cumulative distribution function (CDF) reads The expectation value of the signal event number to set the experimental sensitivity S (B) can be figured out by solving the equations where n p is the smallest number of counts to exclude the null-signal hypothesis at the CL ≥ p, and CDF Poisson = 1 − CDF Poisson is the complementary function of the CDF.
In the extreme background-free case, any positive signal events mean a discovery, i.e. the required number of counts is always n p = 1 regardless of p, and q can be interpreted as the probability that an experiment can report a positive signal (otherwise null signal). The median sensitivity usually adopted in the literature refers to a discovery probability of 50%, which is quite reasonable when the background is large and events are Gaussian distributed.
In the background-free scenario, however, this implies that there is a probability of 50% to observe null signal; therefore a larger value of q (e.g. 90%, 95%, etc) should be adopted such that the obtained sensitivity is more solid and close to the true case.
With the above definitions in mind, we are ready to derive the sensitivities to T 0ν 1/2 and |m ββ | for any given experimental setup of the exposure and background index. For later JHEP03(2021)084 convenience, we introduce the effective exposure ξ ≡ ξ · and the effective background index b ≡ b/ such that the signal detection efficiency in eq. (2.1) is no longer present explicitly. In the left panel of figure 1, the relationship between T 0ν 1/2 and |m ββ |, as indicated in eq. (1.1), has been shown for the nuclear isotope 136 Xe. In the numerical calculations, the phase-space factor G 0ν = 3.79 × 10 −14 yr −1 with the axial vector coupling g A = 1.27 has been used [33][34][35], while the NME for the 0νββ decays of 136 Xe has been taken from table II of ref. [32], where one can find that the theoretical predictions for NME via different methods span a wide range of 1.68 |M 0ν | 4.20. The gray band in the left panel of figure 1 shows the NME uncertainties, and seven typical values of NME have been plotted as dashed lines. Therefore, depending on the NME, |m ββ | = 10 meV and |m ββ | = 1 meV correspond to the half-life of T 0ν 1/2 ∈ 3.8 × 10 27 · · · 2.5 × 10 28 yr and T 0ν 1/2 ∈ 3.8 × 10 29 · · · 2.5 × 10 30 yr, respectively. This observation can be perfectly understood by noting that the half-life T 0ν 1/2 is inversely proportional to |m ββ | 2 as in eq. (1.1). The experimental sensitivity to |m ββ | should be given together with specific values of p and q defined in eq. (2.3). For instance, the sensitivity at the 3σ CL for 90% experiments corresponds to p = 99.73% and q = 90%. In the right panel of figure 1, we have displayed the contours of the sensitivity to |m ββ | at the 3σ CL for 90% identical experiments in the plane of the effective exposure ξ and the effective background index b . Note that the discrete characteristic of the Poisson distribution becomes apparent when the background counts reach the threshold value of O(1). Similar results can also be found in figure 3 of ref. [30] with a different concerned region of ξ and b and a smoothed Poisson distribution. The dashed curves denote the contours of the sensitivity to |m ββ |, where the white colored number in the band is calculated with the average NME value |M 0ν | = 2.94. Each colored band stands for the same sensitivity to |m ββ | as indicated, and its width signifies the NME uncertainty. It should be noticed that to achieve the sensitivity of |m ββ | = 1 meV is very challenging. Even with a background-free environment (namely, b → 0), an effective exposure of at least ξ 300 ton · yr is required to reach the sensitivity of |m ββ | = 1 meV at the 3σ CL.

The Bayesian approach
The Bayesian statistics provides us a logical and practical approach to comparing different models as well as inferring the posterior probability distributions of model parameters. According to the Bayesian theorem, the posterior probability of a hypothesis in light of the experimental data D is where H i stands for the hypothesis with i being the index of different models, and P (H i ) is the prior probability for the model to be true. In addition, P (D|H i ) is identical to the so-called evidence Z i , which is the total likelihood to observe D given the hypothesis H i , and P (D) = i P (D|H i )P (H i ) can be regarded as a normalization factor that fixes JHEP03(2021)084 The model favored by the experimental data among a set of models can be selected by taking their posterior ratios, i.e., If we assume no prior preference for any models, the Bayes factor B ≡ Z i /Z j can directly reflect the odds of different models. We will adopt the Jeffreys scale [36] to interpret the Bayes factor. The posteriors in the parameter space of a specific model can also be updated in light of the experimental data. The posterior probability distribution of the model parameter set Θ can be derived according to where P (D|H i , Θ) denotes the likelihood function in the assumption that the model H i with the parameter set Θ is true, and P (Θ|H i ) is the prior probability of Θ. Here P (D|H i ) is the aforementioned evidence Z i , which can be obtained by integrating over all model parameters, We will use the MultiNest routine for the Bayesian analysis [37][38][39].

Sensitivities to m 1 , ρ and σ
The next-generation 0νββ-decay experiments aim to cover entirely the whole range of |m ββ | in the IO case. The lower boundary of |m ββ | is always lying above 10 meV, which will be taken as a representative value for the sensitivity of next-generation experiments to |m ββ |.
The target value of |m ββ | ∼ 10 meV can be hopefully achieved in a number of proposed experiments, e.g., LEGEND [40], CUPID [41], nEXO [42], JUNO Xe-LS [43] and PandaX-III [44]. Moreover, we try to explore the physics potential of the 0νββ-decay experiment with a sensitivity to |m ββ | 1 meV in the NO case. Thus two scenarios will be considered: • Setup-I with the total exposure ξ = 50 ton · 5 yr, the background index b = 1.35 ton −1 · yr −1 , and the detection efficiency = 0.634. Such a setup is inspired by the preliminary study of the future JUNO Xe-LS experiment in ref. [43]. Given this experimental setup, one can derive its projected sensitivity to half-life T 0ν 1/2 = 6.24×10 27 yr at the 3σ CL, which can be translated into the sensitivity to the effective neutrino mass |m ββ | = (7.9 · · · 19.7) meV (depending on the NME) at the same CL.
• Setup-II with the total exposure ξ = 400 ton · 5 yr, the background index b = 0 ton −1 · yr −1 , and the detection efficiency = 1. In comparison with the previous setup, the exposure is now increased by one order of magnitude, while the background is assumed to be vanishing. As it is very challenging in reality to achieve these improvements, this experimental setup may just stand for the ultimate goal of the 0νββ-decay experiments in the far future. With this ideal setup, we find that JHEP03(2021)084 the sensitivity to T 0ν 1/2 is 2.67 × 10 30 yr at the 3σ CL, or equivalently the sensitivity to |m ββ | is (0.38 · · · 0.95) meV.
For each specific experimental setup, one is able to examine its sensitivities to the fundamental parameters, such as the lightest neutrino mass m 1 and two Majorana CP phases {ρ, σ}, which is the main task in this subsection. The Bayesian approach will be implemented to derive the posterior distributions of the model parameters and to select favorable models [36,45]. The strategy for our statistical analysis is outlined as below.
First, we assume that the future experiments would have not discovered any signals of 0νββ decays, so the observed events should be ascribed solely to the background. For the background event number B and a hypothetical signal event number N 0ν , the probability to observe the number n tot of total events in the ROI is determined by the likelihood function where the Poisson distribution is adopted. Second, the prior distributions of two relevant neutrino mixing angles {sin 2 θ 12 , sin 2 θ 13 } and two neutrino mass-squared differences {∆m 2 21 , ∆m 2 31 } are taken to be flat in some ranges, which are chosen to be wide enough to cover the latest global-fit results of all neutrino oscillation data in ref. [26]. The unconstrained Majorana CP phases {ρ, σ} are uniformly distributed in the whole range [0, 360 • ). As a fundamental parameter, the lightest neutrino mass m 1 /eV or its logarithm log 10 (m 1 /eV) will be uniformly distributed in the range of m 1 /eV ∈ [10 −7 , 10] or log 10 (m 1 /eV) ∈ [−7, 1], which will be referred to as the flat or log prior on m 1 in the following discussions. Note that one may also adopt the flat or log prior on the sum of three neutrino masses Σ instead of m 1 . We have numerically checked that with a sensitivity of |m ββ | 10 meV, these two prior options of Σ lead to posteriors very similar to that with a flat prior on m 1 . In connecting the fundamental parameters to the hypothetical signal events in 0νββ-decay experiments, one must specify the phase-space factor G 0ν and the NME value |M 0ν |. In our calculations for 136 Xe, the phase-space factor G 0ν is supposed to be Gaussian distributed with the central value and 1σ error as found in ref. [35], namely G 0ν = 3.79 × 10 −14 yr −1 with an error of 0.1%, while the NME |M 0ν | is uniformly distributed in the range [1.68, 4.20] as obtained in various nuclear models [32]. Now, all the priors of model parameters in our analysis have been specified.
Third, the posterior distributions can be derived by imposing the experimental likelihood information of both the existing data and the simulated data of future 0νββ-decay experiments. To be explicit, the likelihood functions of neutrino oscillation parameters are extracted from the global-fit analysis of ref. [26]. For the likelihood of the future 0νββ-decay experiments, we will generate the Asimov data, for which the simulated event number is the same as the theoretical expectation. In general, one should follow the Feldman-Cousins approach [46] by taking the median projection of Monte Carlo simulations. In the assumption of null signals, we can simply set n tot = B in eq. (3.5). For each set of parameters in the model under test, one can predict the expected number of events N 0ν as explained in the second step and find out its associated likelihood by using eq. (3.5). It is worthwhile to JHEP03(2021)084 mention that one may generate the data based on a true event signal, and the estimation of model parameters in this case can also be studied.
To demonstrate the independent constraining power of future 0νββ-decay experiments on the absolute scale of neutrinos masses, we do not include the likelihood of other existing non-oscillation experiments in limiting the neutrino masses and Majorana CP phases. But they will be included for the discriminations between NO and IO as well as the Majorana and Dirac nature of neutrinos in section 3.3. The likelihood information about the effective neutrino mass m β in beta decays, the effective neutrino mass |m ββ | in 0νββ decays and the sum of three neutrino masses Σ ≡ m 1 + m 2 + m 3 is extracted from the existing beta-decay experiments (i.e., Mainz [10], Troitsk [11] and KATRIN [12]), 0νββ-decay experiments (including GERDA [47], KamLAND-Zen [27], EXO [29] and CUORE [48]) and the cosmological observations [16], respectively.
Assuming a null signal in the aforementioned experimental setup of 0νββ decays and following the strategy outline above, one can set limits on the absolute scale of neutrino masses as well as the Majorana CP phases. In figure 2 For the flat prior, the interval of m 1 is bounded from below because of the prior effect. For the log prior, the upper limit of the credible interval is slightly subject to the ad hoc lower bound when we set the prior. For instance, if we shift this prior bound from 10 −7 eV to 10 −4 eV, the upper limit will be changed from 6.1 meV to 9.2 meV accordingly. The upper bound can be transformed into the limit on the sum of neutrino masses Σ by using the best-fit values of mass-squared differences [26] as for flat prior on m 1 ; Σ ≡ m 1 + m 2 + m 3 < 0.067 eV , for log prior on m 1 , (3.7) which are very competitive with the cosmological bounds of Planck. We should emphasize that the 0νββ-decay experiments can provide a direct information on the absolute scale of the lightest neutrino mass instead of a bound on the sum of JHEP03(2021)084

JHEP03(2021)084
all neutrino masses Σ as in cosmology. The limits on m 1 from future cosmological surveys are not expected to be so strong by transforming from the future cosmological bounds on Σ, e.g., from Σ 0.087 eV to m 1 16 meV at 95% CL [49]. The nullsignal constraints on the Majorana CP phases ρ and σ are rather weak. The 95% credible interval of the strongest one reads ρ ∈ [28 • · · · 342 • ] when we adopt the flat prior on m 1 , while the other Majorana CP phase σ is almost unconstrained.
• With Setup-II corresponding to the sensitivity of |m ββ | ≈ 1 meV, the null-signal simulation will exclude a large fraction of the parameter space of m 1 , ρ and σ. As shown in the second row of figure 2, very informative conclusions can be made in this case. The lightest neutrino mass m 1 can be constrained into the 95% credible range m 1 ∈ [1.6 · · · 7.3] meV , for flat prior on m 1 ; m 1 ∈ [0.3 · · · 5.6] meV , for log prior on m 1 , (3.8) which are much better than other observational constraints from beta decays and cosmology in the foreseeable future. These two limits are mostly stable against a change on the model parameter from log 10 (m 1 /eV) to m 1 /eV in obtaining the credible intervals. Apparently, the lower bounds on m 1 in eq. (3.8) arise from the "well"-like structure of |m ββ | [19,20]. In addition, the constraints on the Majorana CP phases at the 95% CL turn out to be ρ ∈ [148 • · · · 212 • ] , for flat prior on m 1 ; ρ ∈ [138 • · · · 222 • ] , for log prior on m 1 , (3.9) The constraint on ρ is almost independent of the priors on m 1 , which contains only 20% of the whole range of [0 · · · 360 • ). The limits on σ are not so strong for both priors, e.g. σ ∈ [121 • · · · 239 • ] for the log prior on m 1 and basically unconstrained for the flat prior. The correlations among m 1 , ρ and σ agree well with the analytical results in ref. [22].
With those two specific setups, we have shown their constraining power on the lightest neutrino mass m 1 . In more general cases, we vary the exposure ξ with the background-free assumption, and plot the 1σ (solid red lines) and 2σ (dotted red lines) credible intervals of m 1 under the null-signal assumption in figure 3. We can notice an apparent converging behavior for two priors. This observation makes sense for the Bayesian analysis, namely, as more and more data have been collected the impact of priors will eventually fade away. For the case of the log prior on m 1 , a lower bound on m 1 appears only after the 3σ sensitivity of the setup to |m ββ | has reached around 1 meV. This result is quite meaningful, as there is only a very small fraction of the parameter space in the "well"-like structure. See, e.g., blue curves of figure 3 in ref. [50]. One can obtain a lower limit on m 1 only when the parameter space with m 1 → 0 meV is highly disfavored, which requires a sensitivity of |m ββ | 1 meV. In the right panel of figure 3, as the exposure increases, the credible intervals do not strictly shrink for the exposure between 10 2 ton · yr and 10 3 ton · yr. This effect is due to the shift of probability from m 1  structure. Beyond a critical value of the exposure, the credible intervals will gradually stop changing as the well-like structure spans a certain range. The bounds from the KATRIN projection and Planck 2018 results are transformed into those on m 1 and shown as gray horizontal lines for comparison, while the future cosmology sensitivity to m 1 at 95% CL corresponding to σ(Σ) ∼ 14 meV [49] is given as the horizontal blue line. One can clearly note the advantage of 0νββ-decay experiments in probing the absolute scale of neutrino masses when the O(meV) sensitivity is achieved [22].

NO vs. IO and Majorana vs. Dirac
Although it is reasonable to assume that the neutrino mass ordering is normal, as indicated by the latest global-fit analysis of neutrino oscillation data, and massive neutrinos are Majorana particles, as in a class of seesaw models of neutrino masses, we can use the Bayesian approach to perform a model comparison. Such a study is based on no prior preference for NO or Majorana neutrinos and maximizes the information from current and future experimental observations. In the Bayesian analysis the preference of NO over IO can be represented by the Bayes factor B N/I which is defined as the ratio of the evidence of NO to that of IO [45]. In either case of NO or IO, the fundamental parameters are the same, including two neutrino mixing angles, two neutrino mass-squared differences, two Majorana CP phases and the lightest neutrino mass. Given these parameters, one can compute the effective mass |m ββ |, which together with the phase-space factor and the NME will predict the 0νββ-decay rate. Then the likelihoods for the signal events can be calculated for a nominal 0νββ-decay experiment. Finally, the Bayes factor can be conveniently obtained by integrating the JHEP03(2021)084 to be IO in future neutrino oscillation experiments like JUNO, one can discriminate the Dirac hypothesis from the Majorana one with a very high statistical significance in the JUNO Xe-LS experiment of 0νββ decays or other similar experiments with a competitive sensitivity. For the NO case, to have an adequate evidence in favor of the Dirac hypothesis over the Majorana one, one must go beyond Setup-I. We show in figure 4 the Bayesian factor ln(B D/M ) as a function of the exposure with the background-free assumption. It can be clearly observed that to obtain a strong evidence, i.e., ln(B) = 5, the exposure should be as large as 10 3 ton · yr while keeping the background vanishing, which corresponds to the ultimate meV sensitivity of the 0νββ-decay experiment. In other words, if there is null 0νββ-decay signal at the meV frontier, we can then claim that the Majorana nature of neutrinos is excluded with a strong evidence. However, it is worthwhile to stress that the conclusions here are based on the standard mechanism of exchanging three light neutrinos, which may not apply to the 0νββ decays induced by some non-standard physics (e.g., sterile neutrinos and left-right symmetric models [2,34]).

Summary
In order to explore the physics potential of future 0νββ-decay experiments with a sensitivity of |m ββ | ≈ 1 meV, we have investigated the projected constraints on the lightest neutrino mass m 1 and the Majorana CP phases, in the assumption of a null signal. For comparison, the experimental setup for the sensitivity of |m ββ | ≈ 10 meV is also considered. Our main results and conclusions are summarized in eqs. (3.6)-(3.9), where the Bayesian approach is adopted for statistical analysis.
We believe that our analysis is very important and suggestive for setting up the future program for 0νββ-decay experiments. As already pointed out in ref. [22], if the experimental sensitivity of |m ββ | = 1 meV is ultimately realized, the determination of absolute neutrino masses and the constraints on Majorana CP phases are very promising, which cannot be reached in other types of future neutrino experiments. We have examined these issues in a quantitative way by performing the Bayesian analysis. Furthermore, the determination of neutrino mass ordering and the Majorana or Dirac nature of massive neutrinos are also studied. Certainly, to achieve all these goals, one has to make great efforts in increasing the target mass and reducing the background by two orders of magnitude compared to the present design of next-generation 0νββ-decay experiments. These technical challenges will be left for more future works [51].

JHEP03(2021)084
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.