Inference of neutrino nature and Majorana CP phases from $0\nu\beta\beta$ decays with inverted mass ordering

Whether the neutrino mass ordering is normal or inverted remains an experimentally open issue in neutrino physics. The knowledge of neutrino mass ordering has great importance for neutrinoless double-beta ($ 0\nu\beta\beta$) decay experiments, which can establish the nature of massive neutrinos, i.e., whether they are Dirac or Majorana fermions. Recently, the KamLAND-Zen 800 measurement has reached for the first time the parameter space of the inverted ordering with a vanishing lightest neutrino mass. By assuming the inverted ordering, we attempt to derive the physical information of the neutrino nature and Majorana CP phases from a negative or positive observation of $ 0\nu\beta\beta$ decays in the near future. Moreover, the possibility of extracting the nuclear matrix element in the case of a positive observation is also examined. To avoid the ambiguity from unknown priors of neutrino masses, we adopt the maximum likelihood method instead of the Bayesian approach usually considered in previous works.


I. INTRODUCTION
The nature of massive neutrinos, whether they are Majorana or Dirac fermions, has been a longstanding quest in neutrino physics [1]. At present, the most feasible process that can uncover the neutrino nature is the neutrinoless double-beta (0νββ) decay (for recent reviews, see Refs. [2][3][4][5][6]). Such a process takes place through violating the lepton number by two units [7], i.e., (Z, A) → (Z + 2, A) + 2e − , where Z and A represent the atomic and mass numbers of a nucleus, respectively. For a given 0νββ decay process, the half-life of the isotope of concern can be calculated by [8] 1 T 0ν where G 0ν denotes the kinematic phase-space factor, M 0ν is the nuclear matrix element (NME), and M e ≈ 0.51 MeV is the electron mass. In the standard threeflavor paradigm, the half-life of 0νββ decays is controlled by the effective Majorana neutrino mass, 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 m i (for i = 1, 2, 3) stand for the absolute masses of three neutrinos, {θ 12 , θ 13 } are the leptonic mixing angles with the standard parametrization, and {ρ, σ} are the Majorana CP phases. As the overall CP phase of m ee is non-physical, it is a convenient choice to assign two Majorana CP phases to the lightest and the heaviest neutrino states [9][10][11], such that only one phase parameter will survive when the lightest neutrino mass is vanishing.
In order to establish the relation between the observable T 0ν 1/2 and the quantity of our interest m ee , it is necessary to have a robust theoretical calculation of NMEs (see, e.g., Refs. [34][35][36] for recent reviews). Interestingly, it was pointed out recently that a short-range effect at the leading order can potentially increase the NME by 40% [37][38][39][40]. Recent developments, especially those in the ab initio approach [41] and the lattice-QCD calculation [42], might allow the nuclear physicists to have a more robust evaluation of nuclear matrix elements with accessible errors. With that being achieved, one is able to derive much more information for the neutrino nature and CP phases from the 0νββ decay data. In the pessimistic case, the anarchic status of NMEs based on various phenomenological evaluations will persist, if these advances encounter possible computational bottleneck when adapting to the heavy nuclei. One may, in such a case, go the other way around to constrain NMEs with 0νββ decay data, in the assumption that our understanding of the neutrino mass pattern is further improved in the future.
Depending on whether neutrino masses are normal ordering (NO) or inverted ordering (IO), the predictions of |m ee | differ significantly. In fact, it has been recognized that this feature allows us to distinguish NO and IO based on the 0νββ measurement alone [15]. However, the determination of mass ordering, with a high significance, is already achievable by upcoming neutrino oscillation experiments such as JUNO [43], DUNE [44] and KM3NeT (ORCA) [45] in the foreseeable future. If the mass order- The effective neutrino mass |m ee | as a function of the lightest neutrino mass (upper two panels) or the Majorana CP phase (lower two panels). The darker blue (for NO) and red (for IO) regions signify the 3σ uncertainties (χ 2 osc = 9) caused by neutrino oscillation parameters, i.e., θ 13 , θ 12 , ∆m 2 sol and ∆m 2 atm , and the lighter blue and red regions stem from the unknown Majorana CP phases ρ and σ. The left two panels are calculated by adopting the NuFIT5.1 global-fit results, while the right two panels give the prospects of 3σ uncertainties after JUNO running for 6 years. For comparison, the KamLAND-Zen limit on |m ee | and projections of several future experiments (SNO+ Phase II [31], LEGEND [32] and nEXO [33]) at 90% confidence level (C. L.) are shown as the horizontal black and orange lines (with the most optimistic NME). For the lower two panels, the lightest neutrino mass m L is taken to be vanishing.
In the upper two panels of Fig. 1, we give |m ee | as a function of the lightest neutrino mass m L for NO (in blue) and IO (in red). The optimistic limit of KamLAND-Zen, |m ee | < 36 meV at 90% C. L., has touched for the first time the IO prediction with a vanishing m L [71]. The next-generation proposals with a sensitivity of |m ee | ∼ O(10 meV) at 90% C. L., such as CUPID (8.4-14 meV) [72], JUNO-ββ (5-12 meV) [73], KamLAND2-Zen (< 20 meV) [74], LEGEND 1000 (8.5 − 19.4 meV) [32], nEXO (5.7 − 17.7 meV) [33], NEXT 1T [75], PandaX-III 1k (20 − 46 meV) [76] and SNO+ Phase II (19 − 46 meV) [31], are expected to cover the full parameter space of IO (with the most optimistic NME), while a large fraction of NO parameter space remains untouched. On the other side, the cosmological observations of Planck [77] already pushed m L to the region where NO and IO predictions start to separate. In the lower panels of Fig. 1, we present |m ee | with respect to the Majorana CP phases. One can see from Eq. (2) that |m ee | depends on only one CP phase when the lightest neutrino mass is vanishing, i.e., m 1 → 0 (m 3 → 0) for NO (IO). This specific example shows visually how a measurement of |m ee | is able to reflect the information of CP phases. Note that in the later numerical analysis, we do not restrict the lightest neutrino mass to be vanishing.
The errors on neutrino oscillation parameters that enter as inputs in |m ee |, are still large even though their measurements have entered the precision era. In this regard, JUNO can not only pin down the crucial neutrino mass ordering but also provide a precise measurement of {θ 12 , ∆m 2 sol , ∆m 2 atm } after 6 years of running (see Table I). The induced 3σ uncertainty in |m ee | is shown as the darker blue and red regions in the left (or right) panel of Fig. 1, using the current global-fit (or JUNO 6-year) results of neutrino oscillation parameters. One can see from the right panel a significant improvement in the predicted regions of |m ee |. The rest of the uncertainty in |m ee | (in lighter blue and red) is attributed to the completely unknown Majorana CP phases. Similar observations have also been made in Refs. [23,24,78].
In the assumption of NO, a meV sensitivity to |m ee | is required to extract possible information of Majorana CP phases [14,18,19,[22][23][24][25], which is, however, far beyond the capability of the next-generation experiments in the near future. One would wish that IO is true from a practical point of view, as the full parameter space of IO will be covered with a 10 meV sensitivity in the not-too-distant future. Standing in this paradoxical situation, we ask ourselves what we can infer quantitatively if the next-generation 0νββ decay experiments come out with positive or null signals, assuming JUNO data point towards IO along with improved measurements of oscillation parameters.
We structure the article as follows. In Sec. II, we describe the statistical framework for our purpose. Our quantitative results are presented in Sec. III, including the implications for the nature of neutrinos (Sec. III A), the Majorana CP phases (Sec. III B) in the case of Majorana neutrinos as well as the NME (Sec. III C). Finally, we conclude in Sec. IV.

II. STATISTICAL FRAMEWORK
While it is straightforward to establish the one-to-one analytical correspondence between T 0ν 1/2 and |m ee | for any specified NMEs, a more quantitative analysis can only be made within a given statistical framework. In the present work, we use the frequentist approach to test different models and to estimate model parameters.
In contrast, previous works were mainly built on the Bayesian approach [15,16,19,20,24,25,53,56], which, however, severely relies on the unknown priors of model parameters. The prior of a parameter in Bayesian statistics should be assigned according to its physical origin. The most pronounced ambiguity for 0νββ decays comes from the prior of neutrino masses. Indeed, the origin of neutrino masses is not yet known, and different ad hoc choices, usually depending on one's belief, can yield very distinct results [15,20,25,56,63,79]. Moreover, there is also ambiguity in the priors of angles and phases. For instance, assigning a flat prior to ρ or sin ρ will produce different interpretations. The absence of a trustworthy prior choice leads to the considerations of objective (uninformative) priors when evaluating a specific parameter of concern [63,80]. However, it is not a trivial task to find such priors, and the priors objective for the estimation of one parameter might be biased for another. To avoid the ambiguity of the prior choice, we hence adopt the maximum likelihood method in line with frequentist's taste.
A simplified statistical analysis just counts the number of ββ events within the region of interest (ROI) without fitting the whole electron spectrum. This should be sufficient for our discussion at the sensitivity level. For a nominal 0νββ decay experiment with a total observed event number N exp tot , the likelihood of a hypothesis H can be evaluated by assuming the Poisson fluctuation, Here, N th 0ν is the theoretical expectation for the 0νββ decay event number of the hypothesis H, and B is the background expectation. It is convenient to work on the basis of log-likelihood −2 ln L 0ν which can be identified as the chi-square χ 2 if the data sample is large and approximately Gaussian distributed. On top of the log-likelihood from 0νββ decays, we further add the information about neutrino masses and mixing angles from neutrino oscillation experiments, beta decays and cosmological observations, namely The cosmological likelihood L cosmo is obtained by analyzing the Markov chain file from the Planck Legacy Archive corresponding to the dataset Planck TT, TE, EE + lowE + lensing + BAO, which has set an upper bound on the sum of neutrino masses Σm i < 0.12 eV [77]. The betadecay likelihood is simply taken from Ref. [96]. Note that the neutrino mass information for our choice is dominated by the cosmological observations. The likelihood of oscillation parameters is constructed by noting the equivalence χ 2 osc = −2 ln L osc assuming Gaussian distributions. Our neutrino oscillation chi-square χ 2 osc is taken 3.03 (4.5 %) 3.12 (5.1 %) 3.04 (4.1 %) 0.5 % · · · · · · sin 2 θ13/10 −2 2.23 (3.1 %) 2.23 (3.0 %) 2.24 (2.8 %) 12% 7 % · · · sin 2 θ23/10 −1 5.69 (5.5 %) 5.78 (5.0 %) 5.70 (5.9 %) · · · 2% 3% δ/π 1.52 (9.0 %) 1.58 (9.0 %) 1.54 (9.0 %) · · · 13 % · · · TABLE I. The global fits and the prospects of neutrino oscillations parameters for IO. The first three columns summarize the results of three global-fit groups in the form of "best fit (1σ precision in %)" [68][69][70]. The last three columns give the prospective precision of future neutrino oscillation experiments JUNO [81], DUNE [44], and HK [82], respectively. as where Θ i osc include {sin 2 θ 13 , sin 2 θ 12 , ∆m 2 sol , ∆m 2 atm }, Θ bf,i osc are the best-fit values, and σ i are the 1σ symmetrized errors. Our best-fit values and the error of sin 2 θ 13 are taken from NuFIT 5.1 [70] for IO, whereas the errors of {sin 2 θ 12 , ∆m 2 sol , ∆m 2 atm } are taken to be the projections after 6-year running of JUNO [81]. Note that by using Eq. (5) we have assumed the oscillation parameters to be independent in their experimental determination. In principle, their possible correlations (though mild in most cases) should be taken into account. Our results in this work would hence represent a conservative assessment, and the inclusion of correlations will strengthen more or less our major conclusions.
For completion, we summarize in Table. I the current global fits and the prospects of neutrino oscillation parameters, in the assumption of IO. The first three columns report the best-fit values along with their uncertainties (in parenthesis) from three major global-fit groups. The uncertainty of each parameter is represented by the 1σ fractional accuracy, defined as 1/6 of the 3σ range divided by the best-fit value. The next three columns show the prospective precision of the upcoming neutrino oscillations experiments, i.e., JUNO [81], DUNE [44], and HK [82], respectively.
Two different models can be compared with the help of likelihood ratio tests, as can be found in the textbook.  The essence is to compute L max (H 0 )/L max (H 1 ) for two models H 0 and H 1 , where the likelihood functions are maximized within each model by adjusting the model parameters. A decent choice of the test statistic is simply given by where the approximation of log-likelihood to a chi-square statistic holds when the sample size is large. This is not true in general for 0νββ decay searches, especially for those proposals with very low background indices, such as LEGEND, for which an O(1) signal event number can suggest a discovery. In such a case, the Monte-Carlo simulation with pseudo-experiments should be invoked to obtain the actual distribution of the test statistic −2∆ ln L for the given model. The p-value to favor a model against another one can be obtained by counting the probability that the test statistic exceeds the experimentally observed −2∆ ln L obs .
To be definitive, we fix the 0νββ-decaying isotope as 136 Xe. The theoretical expectation of events N th 0ν severely replies on the input of NME of the decaying isotope. While the ab initial approach is on its way, the existing phenomenological models available give rather diverse NME predictions, including, e.g., the nuclear shell model (NSM), the quasi-particle random phase approximation (QRPA), the energy-density functional theory (EDF), and the interacting boson model (IBM). Different evaluations [83][84][85][86][87][88][89][90][91][92][93][94][95] are summarized in Table. II following Ref. [2], which can differ by a factor more than four. For our numerical analysis, we treat the 136 Xe NME in two ways: • The NME uniformly (no preference in likelihood) distributes in the range |M 0ν | ∈ (1.11, 4.77). • The NME takes a definitive value, e.g., |M 0ν | = 3 with a negligible uncertainty, which can represent the ultimate goal of future theoretical evaluations. As we will see quantitatively later, a precise knowledge of NME plays a vital role when we extract the neutrino information from the 0νββ decay data.
The above considerations basically set up our statistical framework to compare models and to estimate model parameters. For the convenience of later discussions, we The experimental configuration is fixed as ξ = 5 ton · yr and b = 1 ton −1 · yr −1 . In the left panel, the NME is assumed to be completely known and fixed as |M 0ν | = 3, while in the right panel, the NME takes a flat likelihood in the range of |M 0ν | ∈ (1.11, 4.77) during the fit. In each panel, the upper half one is obtained if a null signal has been observed (i.e., N exp tot = B) in a future 0νββ decay experiment, while the lower half stands for the scenario where a positive excess of events has been observed (i.e., N exp tot − B > 0). The observed value of the test statistic is marked by the dashed vertical line. In the case of a positive excess, we assume the observed signal event number to be N exp tot − B ≈ 27 (can be read from Eq. (7)), which is the expectation value by taking m 3,true = 0 eV, ρ true = 90 • , and |M 0ν | true = 3. The model parameters of pseudo-experiments for all the hypotheses are located by maximizing the likelihood with respect to a null or positive signal.
further establish the connection between a given experimental configuration and its sensitivity to T 0ν 1/2 . As is well known, the sensitivity of a generic 0νββ decay experiment to the half-life T 0ν 1/2 can be estimated via the relation [20] T 0ν 1/2 = ln 2 · where N A is the Avogadro constant, ξ represents the exposure, is the event detection efficiency, m iso is the molar mass of the specified isotope. Here, S(B) represents the required expectation of signal event number within the ROI, for which a fraction q of a set of identical experiments can report a discovery of 0νββ decays at a C. L. larger than p. The dependence of S(B) on B = ξ · b with b being the background index is explicitly shown.

III. INFERENCE FROM A NULL OR POSITIVE SIGNAL
Without otherwise specified, we will adopt an idealized experimental setup for the numerical assessment, assum-ing a background index b = 1 ton −1 · yr −1 and a full detection efficiency = 100%. The chosen background index corresponds to one count of background events for one ton · yr of exposure, which represents an ultralow background level achievable in the future. The only remaining experimental parameter, which is allowed to freely vary, is the target exposure ξ.

A. Dirac or Majorana nature
The first problem we want to address is to what extent the Majorana nature of neutrinos will be verified or falsified if the neutrino mass ordering turns out to be inverted. Following the statistical framework outlined in Sec. II, our task is then to compare the Dirac neutrino hypothesis H D to the Majorana one H M . There are two possible experimental outcomes for the next-generation 0νββ decay experiments: • A null signal. Namely, the observed event number N exp tot coincides with the expectation of background B. During the fit, we first need to obtain the best-fit parameter values (e.g., m 3 , ρ and θ 12 ) The positive signal is taken to be the expectation by choosing m 3,true = 0 eV, ρ true = 90 • and |M 0ν | true = 3 with other oscillation parameters in their NuFIT 5.1 best fits. As usual, the background index is fixed as b = 1 ton −1 · yr −1 with a full detection efficiency. The 90% C. L. sensitivity to |m ee | for the corresponding exposure is indicated on the top axis. The sensitivities of LEGEND [97] and nEXO [33] to |m ee | with averaged NMEs are shown as the vertical orange lines, which should be interpreted with the top axis.
in H D and H M , respectively, by maximizing the total likelihood in Eq. (4). The observed value of the test statistic −2∆ ln L obs = −2 ln(L D /L M ) is then calculated with the best-fit parameters within each model. To find out the p-value of the observed test statistic for a given hypothesis, we generate −2∆ ln L as in Eq. (6) with pseudo-experiments assuming either H D or H M is true. An example of the probability distribution of −2∆ ln L is illustrated in the upper two panels of Fig. 2. The observed −2∆ ln L obs is marked by the vertical dashed lines. The experimental exposure is chosen to be ξ = 5 ton · yr. For the left panel, the NME is assumed to be completely known and taken as |M 0ν | = 3, and the p-value to accept the Majorana hypothesis is only 2%. We notice that the distributions of Dirac and Majorana hypotheses are well separated (meaning that we can distinguish them well) if we have a plausible |M 0ν | calculation. Whereas, for the right panel, the NME is varying freely in the range |M 0ν | ∈ (1.11, 4.77) during the fit, which severely dilutes the power of discrimina-tion * , and a larger exposure is required to well separate them.
• A positive signal. One can do the exercise conversely, by assuming that a positive signal excess indicating Majorana neutrinos arises. To fix N exp tot in this nominal analysis, we first assume that neutrinos are Majorana fermions and choose a benchmark parameter choice for H M : m 3,true = 0 eV, ρ true = 90 • and |M 0ν | true = 3 with other parameters in their NuFIT best fits. The resultant expectation of N th 0ν is then taken to be the observed signal event, i.e., N exp tot − B = 27 with B = 5, a rather optimistic expectation with IO. Similar to the null signal case, the Monte-Carlo results of the test statistic −2∆ ln L are given in the lower two panels of Fig. 2. In this example, the distributions of Dirac and Majorana hypotheses are almost completely separated, and one is able to reject Dirac (or accept Majorana) with a very high significance. Interestingly, even in the case with an unknown NME, we can reject the Dirac hypothesis with a significance similar to the case with a fixed NME. This is understandable as a positive signal will force the likelihood maximizer to find the true |M 0ν | for the Majorana hypothesis, in comparison to the case of null signal.
The above examples of discriminating Dirac and Majorana neutrinos have assumed a fixed experimental configuration. It is straightforward to explore the general sensitivity by varying the experimental exposure. We show in Fig. 3 the p-value as a function of the experimental exposure ξ, in order to reject H M over H D with a null signal, and vise versa. The corresponding 1σ, 2σ and 3σ significance levels of p-value are indicated by the horizontal gray lines. The red rhombus stands for the sensitivity to reject H M in the presence of the null signal by taking |M 0ν | = 3, while the black rhombus is for the varying NME |M 0ν | ∈ (1.11, 4.77). The blue triangle stands for the sensitivity to reject H D in the presence of the positive signal. The event number of the positive signal is taken to be the expectation by choosing m 3,true = 0 eV, ρ true = 90 • and |M 0ν | true = 3, which gives N exp tot − B ≈ 27 · ξ/(5 ton · yr). It should be noted that this parameter choice is only used to generate the nominal events, since the actual 0νββ decay experiments of concern have not yet given any data. This is needed for the demonstrative purpose, but a different choice of the signal event number will more or less alter the results. After the nominal events are fixed, during the analysis we vary the relevant model parameters freely and confine them with the total likelihood in Eq. (4). A different choice of m 3,true in generating the nominal events does not make too much difference for the major results, because we can see from Fig. 1 that the Planck experiment * The fit of the null signal will always drive NME to the smallest possible value |M 0ν | = 1.11. already pushed m 3 to the regime where the dependence of |m ee | on m 3 becomes weak. In Fig. 3, the sensitivity to |m ee | for the corresponding exposure with b = 1 ton −1 · yr −1 is indicated on the top axis. For a given exposure, this is obtained by converting T 0ν 1/2 in Eq. (7) with |M 0ν | = 3. We take the median sensitivity (q = 50%) to have a 90% C. L. discovery potential. For comparison, the sensitivities of two projects, i.e., LEGEND (8.5 − 19.4 meV) [97] and nEXO (5.7 − 17.7 meV) [33], to |m ee | with averaged NMEs are also indicated by the vertical lines. Note that the LEG-END and nEXO lines should be interpreted with the help of the top axis. The exact capability of LEGEND and nEXO in distinguishing Dirac and Majorana hypotheses should have been performed by using their detailed experimental configurations. However, as we mentioned, to simplify the analysis we have adopted the framework in Sec. II with the background index b = 1 ton −1 · yr −1 and the full detection efficiency = 100%. In this case, the vertical lines should stand for the experiments with our idealized setup, which have equivalent |m ee | sensitivities as the true LEGEND and nEXO experiments.
Some further remarks are as follows.
• If the uncertainties in NMEs are negligible in the future and take the average of current theoretical evaluations (|M 0ν | ≈ 3), the future projects such as LEGEND and nEXO can distinguish the Dirac and Majorana hypotheses with more than 2σ C. L. This can be seen by reading off the p-value of the red rhombus and blue triangle intersecting with the LEGEND and nEXO lines. The 2σ C. L. corresponds to a rejection p-value around 5%. • Pessimistically, one may presume that no further progress is achieved in the NME evaluation in the future. In this case, even with the sensitivity of nEXO one cannot tell whether neutrinos are Dirac or Majorana with a considerable significance, i.e., at most 1σ, if the Dirac hypothesis is true. Whereas, if the Majorana hypothesis is true, the discrimination significance will rely on how much the actual NME is. We do not necessarily need a precise knowledge of NME to reject the Dirac neutrino hypothesis if the actual NME is large, in which case the actual event excess will be large.

B. Majorana CP phases
If a positive signal excess is observed in the future, the most important message will be that neutrinos are Majorana fermions. Other than that, one can also attempt to extract additional unique information. Because neutrino oscillation parameters have been well measured and will be further improved by JUNO in the very near future, the observable quantity T 0ν 1/2 is mainly a function of the lightest neutrino mass m 3 , the CP phases and the NME |M 0ν |. Furthermore, the sum of neutrino masses is severely constrained by the cosmological observations if no new physics is present, which limits the magnitude of m 3 such that the dependence of T 0ν 1/2 on m 3 variation becomes weak. Hence, from the measurement of T 0ν 1/2 one can simultaneously constrain the Majorana phase ρ in the neutrino sector and |M 0ν | in the nuclear sector.
We start by exploring the sensitivity to ρ by making assumptions about the theoretical input of |M 0ν |. In principle, we will need a positive 0νββ decay signal to set bounds on ρ, since a null signal can always be interpreted as a sign of Dirac neutrinos. Assuming a nominal event number N exp tot has been observed, we can set constraints on ρ as follows. First, we scan over all the model parameters and calculate the corresponding likelihood L(θ 13 , θ 12 , ∆m 2 sol , ∆m 2 atm , m 3 , ρ, σ, |M 0ν |) with Eq. (4). Second, for every ρ we marginalize over all the other parameters to obtain the maximum L (or the minimum −2 ln L). Finally, we obtain the likelihood profile −2∆ ln L as a function of ρ by subtracting its global minimum. The likelihood profiles of ρ with various nominal inputs are given in Fig. 4 for the experimental exposure ξ = 20 ton · yr. Four benchmark choices are illustrated: assuming N exp tot − B to be vanishing (gray curves) or to be the expectation values of m 3,true = 0 eV and |M 0ν | true = 3 with ρ true = 0 • (dashed red curves), 90 • (dotted pink curves) and 180 • (solid red curves), respectively. The 3σ C. L. range of ρ can be obtained by setting −2 ln L M = 9. In the absence of any signal events, one can still constrain ρ because smaller |m ee | will be preferred. However, as has been mentioned, we will be favoring Dirac over Majorana hypothesis in this case.
As usual, we investigate two extreme scenarios of 136 Xe NME, i.e., |M 0ν | = 3 with a negligible uncertainty (left panel) and |M 0ν | ∈ (1.11, 4.77) (right panel). We find that the inclusion of the NME uncertainty will significantly dilute the sensitivity to ρ due to the understandable degeneracy. This result further encourages the efforts from the nuclear community to push forward the plausible calculation of |M 0ν |.
The dependence of the sensitivity on the nominal exposure ξ is shown in Fig. 5, where we give the 3σ allowed range of ρ with a signal event number calculated from m 3,true = 0 eV and |M 0ν | true = 3 with ρ true = 0 • (top panel), 90 • (middle panel) or 180 • (bottom panel). In the left panel, the NME is assumed to be known as |M 0ν | = 3, and next-generation 0νββ decay experiments will be able to rule out a significant amount of parameter space of the Majorana CP phase ρ. In comparison, for the right panel we observe that the lack of knowledge of |M 0ν | significantly reduces the constraining power to ρ, even with an exposure of 100 ton · yr and an ultra-low background of b = 1 ton −1 · yr −1 .
The measurement of Majorana CP phases is of great theoretical importance. For instance, the µ−τ reflection symmetry [98] by far remains one of the best flavor symmetry which is consistent with the neutrino oscillation data. This symmetry leads to very sharp predictions for the Majorana CP phases (ρ, σ = 0 • or 180 • ) [21,99]. From the left panel of Fig. 5, we notice that if one of the solution is true (e.g., ρ true = 0 • ), the remaining one (ρ true = 180 • ) can be ruled out at 3σ C. L. with an exposure of only 1 ton · yr. However, this capability is washed out, if the NME uncertainty persists in the future. In practice, the actual constraining power might be between those two extreme cases, with a reduced but not completely negligible NME uncertainty.

C. Nuclear matrix element
The NME, as an imperative theoretical input, can be conversely "measured" by 0νββ decay experiments. If IO turns out to be true, the 0νββ decay will become an excellent ruler for the NME, though the resolution of this ruler is limited by the unknown CP phases in |m ee |. For NO, the uncertainty of |m ee |, mainly due to the "well structure", will be too large to extract any upper bounds on NME. Nevertheless, deriving a lower bound is always possible for NO with a positive observation of 0νββ decays, which is, however, not the main concern of the 3σ Null signal Positive, |M 0ν ture =1.11 Positive, |M 0ν ture =3 Positive, |M 0ν ture =4.77 The 3σ sensitivities to the NME in terms of the 136 Xe exposure ξ. We assume that a positive signal has been observed in the experiment. The signal event number takes the expectation value of the choice m 3,true = 0 eV and ρ true = 90 • with |M 0ν | true = 1.11 (top panel), 3 (middle panel) or 4.77 (bottom panel). present work.
Following a procedure similar to fitting the CP phase, in Fig. 6 we generate the likelihood profile as a function of |M 0ν |. For demonstration, several benchmark choices are shown, including a null signal (gray curve) and a positive signal assuming m 3,true = 0 eV and ρ true = 90 • with |M 0ν | true = 1.11 (dashed blue curve), 3 (dotted blue curve) or 4.77 (solid blue curve).
In Fig. 7, we present the 3σ sensitivities to the NME with respect to the nominal exposure ξ. The true value of NME is taken to be |M 0ν | true = 1.11 (top panel), 3 (middle panel) or 4.77 (lower panel). We have several remarks.
• In the given example, if the true value of |M 0ν | is 3, a 3σ C. L. constraint 0.7 < |M 0ν | < 3.5 can be obtained with an exposure of ξ = 100 ton · yr, which is slightly narrower than the current prediction |M 0ν | ∈ (1.11, 4.77). • The resolution for NME is limited by the unknown Majorana CP phases. Hence, after the exposure increases to a certain level, the sensitivity will not be improved any further. Since there are not other foreseeable places to obtain the information of Majorana CP phases, this situation will persist regardless of the experimental achievements of 0νββ decays.

IV. CONCLUSION
The next-generation 0νββ decay experiments are designed for a target sensitivity around |m ee | ≈ 10 meV, which sets the lower boundary of |m ee | if the neutrino mass ordering is inverted. However, if NO is true and the lightest neutrino is unfortunately tiny, those nextgeneration projects will expect no 0νββ decay signals in the absence of new physics. While the mass ordering is still an open issue [63] to be addressed by the upcoming neutrino oscillation experiments, we present here a frequentist analysis about what we can extract from 0νββ decays if IO eventually wins. From a positive or negative observation of 0νββ decays in the future, the physical information we can extract is at least three-fold. First and foremost, the nature of massive neutrinos can be well discriminated in the assumption of IO. Second, the Majorana CP phase ρ can be probed with a positive 0νββ decay signal. Last, one can conversely constrain the NME in the presence of a signal. The role of NME in interpreting the unknowns in the neutrino sector is critical. To efficiently probe the nature of neutrinos and the Majorana CP phases, a robust NME evaluation beyond the state of the art is necessary. Taking 136 Xe for instance, if the actual NME is found to be |M 0ν | = 3 with a negligible uncertainty, the experimental proposals like nEXO with a 10 meV sensitivity can distinguish Dirac and Majorana neutrinos with 3σ C. L. or so, in the presence of a positive or null signal (as shown by Fig. 3). If |M 0ν | remains vague, for a null signal one cannot meaningfully reject the Majorana hypothesis with the next-generation design. Similar conclusions hold for the Majorana CP phases. The improved calculations of NMEs will make it possible to constrain one of the Majorana CP phases ρ (see Fig. 5). In addition, in the worst case with a completely unknown NME, one can do the exercise conversely by constraining the NME with a positive 0νββ decay signal in the future (see Figs. 6 and 7).