1 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\nu \beta \beta \)) 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) \rightarrow (Z+2, A) + 2 e^{-}\), where Z and A represent the atomic and mass numbers of a nucleus, respectively. For a given \(0\nu \beta \beta \) decay process, the half-life of the isotope of concern can be calculated by [8]

$$\begin{aligned} \dfrac{1}{T^{0\nu }_{1/2}} = G^{}_{0\nu } \cdot \left| { M}^{}_{0\nu }\right| ^{2} \cdot \frac{\left| m_{e e} \right| ^{2} }{M^2_e } , \end{aligned}$$
(1)

where \(G_{0\nu }\) denotes the kinematic phase-space factor, \( M_{0\nu } \) is the nuclear matrix element (NME), and \(M^{}_{e} \approx 0.51~\mathrm{MeV}\) is the electron mass. In the standard three-flavor paradigm, the half-life of \(0\nu \beta \beta \) decays is controlled by the effective Majorana neutrino mass,

$$\begin{aligned} {m}^{}_{ee}\equiv & {} m^{}_1 \cos ^2 \theta ^{}_{13} \cos ^2 \theta ^{}_{12} e^{\mathrm{i}\rho } + m^{}_2 \cos ^2 \theta ^{}_{13} \sin ^2 \theta ^{}_{12} \nonumber \\&+ \ m^{}_3 \sin ^2 \theta ^{}_{13} e^{\mathrm{i}\sigma } , \end{aligned}$$
(2)

where \(m^{}_i\) (for \(i = 1, 2, 3\)) stand for the absolute masses of three neutrinos, \(\{\theta ^{}_{12},\theta ^{}_{13} \}\) are the leptonic mixing angles with the standard parametrization, and \(\{\rho , \sigma \}\) 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.

The importance of \(0\nu \beta \beta \) decay searches for neutrino physics is well recognized [12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30], as it can answer: (i) the nature of neutrinos, (ii) the Majorana CP phases, (iii) the absolute scale of neutrino masses, and (iv) the neutrino mass ordering. While the neutrino masses can be better probed by other observations (e.g., cosmology and oscillation experiments), the \(0\nu \beta \beta \) decay search seems to be the only reliable way towards the neutrino nature and the Majorana CP phases.

Fig. 1
figure 1

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\sigma \) uncertainties (\(\chi ^2_{\mathrm{osc}} = 9\)) caused by neutrino oscillation parameters, i.e., \(\theta ^{}_{13}\), \(\theta ^{}_{12}\), \(\varDelta m^{2}_{\mathrm{sol}}\) and \(\varDelta m^2_{\mathrm{atm}}\), and the lighter blue and red regions stem from the unknown Majorana CP phases \(\rho \) and \(\sigma \). The left two panels are calculated by adopting the NuFIT5.1 global-fit results, while the right two panels give the prospects of \(3\sigma \) 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^{}_{\mathrm{L}}\) is taken to be vanishing

Full size image

In order to establish the relation between the observable \(T^{0\nu }_{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\nu \beta \beta \) 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\nu \beta \beta \) 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\nu \beta \beta \) 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 ordering turns out to be normal, we do not expect any \(0\nu \beta \beta \) decay signals with a high probability in the near future. A positive signal would mean either new physics contributions [13, 46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63] in \(|m^{}_{ee}|\) or an unusual cosmology with large neutrino masses [64,65,66,67]. If the mass ordering turns out to be inverted, in contrast to the current global-fit preference [68,69,70], we will instead have a considerable discovery probability. A positive signal would not only mean neutrinos are Majorana fermions but also allow for constraining the Majorana CP phases and NMEs. A negative result would then favor Dirac neutrinos, contrary to the likes of most neutrino model builders.

In the upper two panels of Fig. 1, we give \(|m^{}_{ee}|\) as a function of the lightest neutrino mass \(m^{}_{\mathrm{L}}\) for NO (in blue) and IO (in red). The optimistic limit of KamLAND-Zen, \( |m_{ee} | < 36~\mathrm{meV}\) at 90% C. L., has touched for the first time the IO prediction with a vanishing \(m^{}_{\mathrm{L}}\) [71]. The next-generation proposals with a sensitivity of \(|m^{}_{ee}| \sim {\mathcal {O}}(10~\mathrm{meV})\) at 90% C. L., such as CUPID (8.4–14 meV) [72], JUNO-\(\beta \beta \) (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 \(\mathrm{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^{}_{\mathrm{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 \rightarrow 0 \) (\(m_3 \rightarrow 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 \(\{\theta ^{}_{12}, \varDelta m^{2}_{\mathrm{sol}}, \varDelta m^2_{\mathrm{atm}}\}\) after 6 years of running (see Table 1). The induced \(3\sigma \) 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~\mathrm{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\nu \beta \beta \) 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 Sect. 2, we describe the statistical framework for our purpose. Our quantitative results are presented in Sect. 3, including the implications for the nature of neutrinos (Sect. 3.1), the Majorana CP phases (Sect. 3.2) in the case of Majorana neutrinos as well as the NME (Sect. 3.3). Finally, we conclude in Sect. 4.

2 Statistical framework

While it is straightforward to establish the one-to-one analytical correspondence between \(T^{0\nu }_{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\nu \beta \beta \) 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 \(\rho \) or \(\sin {\rho }\) 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.

Table 1 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\sigma \) 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
Full size table

A simplified statistical analysis just counts the number of \(\beta \beta \) 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\nu \beta \beta \) decay experiment with a total observed event number \(N^{\mathrm{exp}}_{\mathrm{tot}}\), the likelihood of a hypothesis H can be evaluated by assuming the Poisson fluctuation,

$$\begin{aligned} {{\mathcal {L}}^{}_{0\nu } (N^{\mathrm{th}}_{0\nu }) }= \frac{(N^{\mathrm{th}}_{0\nu }+B)^{N^{\mathrm{exp}}_{\mathrm{tot}}}}{N^{\mathrm{exp}}_{\mathrm{tot}}!} \cdot e^{-(N^{\mathrm{th}}_{0\nu }+B)} . \end{aligned}$$
(3)

Here, \(N^{\mathrm{th}}_{0\nu }\) is the theoretical expectation for the \(0\nu \beta \beta \) 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 {{\mathcal {L}}^{}_{0\nu }}\) which can be identified as the chi-square \(\chi ^2\) if the data sample is large and approximately Gaussian distributed. On top of the log-likelihood from \(0\nu \beta \beta \) decays, we further add the information about neutrino masses and mixing angles from neutrino oscillation experiments, beta decays and cosmological observations, namely

$$\begin{aligned} -2 \ln {{\mathcal {L}}}= & {} -2\ln {{\mathcal {L}}^{}_{0\nu }} -2 \ln {{\mathcal {L}}^{}_{\mathrm{osc}}} \nonumber \\&-2 \ln {{\mathcal {L}}^{}_{ \beta }} -2 \ln {{\mathcal {L}}^{}_{\mathrm{cosmo}}} . \end{aligned}$$
(4)

The cosmological likelihood \({\mathcal {L}}^{}_{\mathrm{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 \(\Sigma m^{}_{i} <0.12~\mathrm{eV}\) [77]. The beta-decay 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 \(\chi ^2_{\mathrm{osc}} = -2 \ln {{\mathcal {L}}^{}_{\mathrm{osc}}}\) assuming Gaussian distributions. Our neutrino oscillation chi-square \(\chi ^2_{\mathrm{osc}}\) is taken as

$$\begin{aligned} \chi ^2_{\mathrm{osc}} \equiv \sum ^{}_{i} \frac{(\varTheta ^{i}_{\mathrm{osc}} - \varTheta ^{\mathrm{bf},i}_{\mathrm{osc}})^2}{\sigma ^2_{i}} , \end{aligned}$$
(5)

where \(\varTheta ^{i}_{\mathrm{osc}}\) include \(\{\sin ^2{\theta ^{}_{13}}, \sin ^2\theta ^{}_{12}, \varDelta m^{2}_{\mathrm{sol}}, \varDelta m^{2}_{\mathrm{atm}}\}\), \(\varTheta ^{\mathrm{bf},i}_{\mathrm{osc}}\) are the best-fit values, and \(\sigma ^{}_{i}\) are the \(1\sigma \) symmetrized errors. Our best-fit values and the error of \(\sin ^2{\theta ^{}_{13}}\) are taken from NuFIT 5.1 [70] for IO, whereas the errors of \(\{ \sin ^2\theta ^{}_{12}, \varDelta m^{2}_{\mathrm{sol}}, \varDelta m^{2}_{\mathrm{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.

Table 2 The nuclear matrix elements \(|M_{0\nu }|\) for \(^{136}\)Xe evaluated with the IBM, QRPA, EDF and NSM methods, adapted from Ref. [2]
Full size table

For completion, we summarize in Table 1 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\sigma \) fractional accuracy, defined as 1/6 of the \(3\sigma \) 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 \({\mathcal {L}}^{\mathrm{max}}({H^{}_0})/{\mathcal {L}}^{\mathrm{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

$$\begin{aligned} -2\varDelta \ln {\mathcal {L}} \equiv -2\ln \frac{{\mathcal {L}}(H^{}_{0})}{{\mathcal {L}}(H^{}_{1})} \approx \chi ^2(H^{}_{0})-\chi ^2(H^{}_{1}) , \end{aligned}$$
(6)

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\nu \beta \beta \) decay searches, especially for those proposals with very low background indices, such as LEGEND, for which an \({\mathcal {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\varDelta \ln {\mathcal {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\varDelta \ln {\mathcal {L}}^{\mathrm{obs}}\).

To be definitive, we fix the \(0\nu \beta \beta \)-decaying isotope as \(^{136}\mathrm{Xe}\). The theoretical expectation of events \(N^{\mathrm{th}}_{0\nu }\) 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 2 following Ref. [2], which can differ by a factor more than four. For our numerical analysis, we treat the \(^{136}\mathrm{Xe}\) NME in two ways:

  • The NME uniformly (no preference in likelihood) distributes in the range \(|M_{0\nu }| \in (1.11, 4.77)\).

  • The NME takes a definitive value, e.g., \(|M_{0\nu }| = 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\nu \beta \beta \) 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 further establish the connection between a given experimental configuration and its sensitivity to \(T^{0\nu }_{1/2}\). As is well known, the sensitivity of a generic \(0\nu \beta \beta \) decay experiment to the half-life \(T^{0\nu }_{1/2}\) can be estimated via the relation [20]

$$\begin{aligned} T^{0\nu }_{1/2} =\ln {2} \cdot \frac{N^{}_{\mathrm{A}} \cdot \xi \cdot \epsilon }{m^{}_{\mathrm{iso}} \cdot S^{}_{}(B)}, \end{aligned}$$
(7)

where \(N^{}_{\mathrm{A}}\) is the Avogadro constant, \(\xi \) represents the exposure, \(\epsilon \) is the event detection efficiency, \(m^{}_{\mathrm{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\nu \beta \beta \) decays at a C. L. larger than p. The dependence of S(B) on \(B = \xi \cdot b\) with b being the background index is explicitly shown.

3 Inference from a null or positive signal

Without otherwise specified, we will adopt an idealized experimental setup for the numerical assessment, assuming a background index \(b=1~\mathrm{ton^{-1} \cdot year^{-1}}\) and a full detection efficiency \(\epsilon = 100\%\). The chosen background index corresponds to one count of background events for one \(\mathrm{ton \cdot year}\) of exposure, which represents an ultra-low background level achievable in the future. The only remaining experimental parameter, which is allowed to freely vary, is the target exposure \(\xi \).

3.1 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 Sect. 2, our task is then to compare the Dirac neutrino hypothesis \({H}^{}_{\mathrm{D}}\) to the Majorana one \({H}^{}_{\mathrm{M}}\). There are two possible experimental outcomes for the next-generation \(0\nu \beta \beta \) decay experiments:

Fig. 2
figure 2

Probability distributions of the test statistic \(-2\ln ({\mathcal {L}}^{}_{\mathrm{D}} / {\mathcal {L}}^{}_{\mathrm{M}})\) obtained from the Monte-Carlo simulations by assuming neutrinos are Dirac (blue histograms) or Majorana (orange histograms) fermions. The experimental configuration is fixed as \(\xi = 5~\mathrm{ton \cdot year}\) and \(b = 1~\mathrm{ton^{-1} \cdot year^{-1}}\). In the left panel, the NME is assumed to be completely known and fixed as \(\left| { M}^{}_{0\nu }\right| = 3\), while in the right panel, the NME takes a flat likelihood in the range of \(\left| { M}^{}_{0\nu }\right| \in (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^{\mathrm{exp}}_{\mathrm{tot}} = B\)) in a future \(0\nu \beta \beta \) decay experiment, while the lower half stands for the scenario where a positive excess of events has been observed (i.e., \(N^{\mathrm{exp}}_{\mathrm{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^{\mathrm{exp}}_{\mathrm{tot}} - B \approx 27\) (can be read from Eq. (7)), which is the expectation value by taking \(m^{}_{3,\mathrm{true}}=0~\mathrm{eV}\), \(\rho ^{}_{\mathrm{true}} = 90^{\circ }\), and \(\left| { M}^{}_{0\nu }\right| ^{}_{\mathrm{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

Full size image

  • A null signal. Namely, the observed event number \(N^{\mathrm{exp}}_{\mathrm{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}\), \(\rho \) and \(\theta ^{}_{12}\)) in \({H}^{}_{\mathrm{D}}\) and \({H}^{}_{\mathrm{M}}\), respectively, by maximizing the total likelihood in Eq. (4). The observed value of the test statistic \(-2\varDelta \ln {\mathcal {L}}^{\mathrm{obs}} = -2\ln ({\mathcal {L}}^{}_{\mathrm{D}}/{\mathcal {L}}^{}_{\mathrm{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\varDelta \ln {\mathcal {L}}\) as in Eq. (6) with pseudo-experiments assuming either \({H}^{}_{\mathrm{D}}\) or \({H}^{}_{\mathrm{M}}\) is true. An example of the probability distribution of \(-2\varDelta \ln {\mathcal {L}}\) is illustrated in the upper two panels of Fig. 2. The observed \(-2\varDelta \ln {\mathcal {L}}^{\mathrm{obs}}\) is marked by the vertical dashed lines. The experimental exposure is chosen to be \(\xi = 5~\mathrm{ton \cdot year}\). For the left panel, the NME is assumed to be completely known and taken as \(|{\mathcal {M}}^{}_{0\nu }|=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 \(|{\mathcal {M}}^{}_{0\nu }|\) calculation. Whereas, for the right panel, the NME is varying freely in the range \(\left| { M}^{}_{0\nu }\right| \in (1.11,4.77)\) during the fit, which severely dilutes the power of discrimination,Footnote 1 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^{\mathrm{exp}}_{\mathrm{tot}}\) in this nominal analysis, we first assume that neutrinos are Majorana fermions and choose a benchmark parameter choice for \({H}^{}_{\mathrm{M}}\): \(m^{}_{3,\mathrm{true}}=0~\mathrm{eV}\), \(\rho ^{}_{\mathrm{true}} = 90^{\circ }\) and \(\left| { M}^{}_{0\nu }\right| ^{}_{\mathrm{true}} = 3\) with other parameters in their NuFIT best fits. The resultant expectation of \(N^{\mathrm{th}}_{0\nu }\) is then taken to be the observed signal event, i.e., \(N^{\mathrm{exp}}_{\mathrm{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\varDelta \ln {\mathcal {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 \(\left| { M}^{}_{0\nu }\right| \) for the Majorana hypothesis, in comparison to the case of null signal.

Fig. 3
figure 3

The significance level to reject the Majorana or Dirac neutrino hypothesis, if a null or positive signal has been observed, as a function of the \(^{136}\)Xe exposure. The positive signal is taken to be the expectation by choosing \(m^{}_{3,\mathrm{true}}=0~\mathrm{eV}\), \(\rho ^{}_{\mathrm{true}} = 90^{\circ }\) and \(\left| { M}^{}_{0\nu }\right| ^{}_{\mathrm{true}} = 3\) with other oscillation parameters in their NuFIT 5.1 best fits. As usual, the background index is fixed as \(b = 1~\mathrm{ton^{-1} \cdot year^{-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

Full size image

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 \(\xi \), in order to reject \({H}^{}_{\mathrm{M}}\) over \({H}^{}_{\mathrm{D}}\) with a null signal, and vise versa. The corresponding \(1\sigma \), \(2\sigma \) and \(3\sigma \) significance levels of p-value are indicated by the horizontal gray lines. The red rhombus stands for the sensitivity to reject \({H}^{}_{\mathrm{M}}\) in the presence of the null signal by taking \(\left| { M}^{}_{0\nu }\right| =3\), while the black rhombus is for the varying NME \(\left| { M}^{}_{0\nu }\right| \in (1.11, 4.77)\). The blue triangle stands for the sensitivity to reject \({H}^{}_{\mathrm{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,\mathrm{true}}=0~\mathrm{eV}\), \(\rho ^{}_{\mathrm{true}} = 90^{\circ }\) and \(\left| { M}^{}_{0\nu }\right| ^{}_{\mathrm{true}} = 3\), which gives \(N^{\mathrm{exp}}_{\mathrm{tot}} - B \approx 27 \cdot \xi /(5~\mathrm{ton \cdot year})\). It should be noted that this parameter choice is only used to generate the nominal events, since the actual \(0\nu \beta \beta \) 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,\mathrm{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 already pushed \(m^{}_{3}\) to the regime where the dependence of \(|m^{}_{ee}|\) on \(m^{}_{3}\) becomes weak.

Fig. 4
figure 4

The likelihood profile with respect to the Majorana CP phase \(\rho \), by assuming \(\left| { M}^{}_{0\nu }\right| = 3\) (left panel) or \(\left| { M}^{}_{0\nu }\right| \in (1.11, 4.77)\) (right panel) during the fit. Here, the experimental configuration is fixed as \(\xi = 20~\mathrm{ton \cdot year}\) and \(b = 1~\mathrm{ton^{-1} \cdot year^{-1}}\). The observed number of signal events \(N^{\mathrm{exp}}_{\mathrm{tot}} - B\) is assumed to be vanishing (gray curves) or to be the expectation values of \(m^{}_{3,\mathrm{true}} = 0~\mathrm{eV}\) and \(\left| { M}^{}_{0\nu }\right| ^{}_{\mathrm{true}} = 3\) with \(\rho ^{}_{\mathrm{true}} = \) \(0^{\circ }\) (dashed red curves), \(90^{\circ }\) (dotted pink curves) and \(180^{\circ }\) (solid red curves), respectively

Full size image
Fig. 5
figure 5

The \(3\sigma \) sensitivities to the Majorana CP phase \(\rho \) with respect to the \(^{136}\)Xe exposure \(\xi \), by assuming the NME to be completely known as \(\left| { M}^{}_{0\nu }\right| = 3\) (left panel), or spreading over a range \(\left| { M}^{}_{0\nu }\right| \in (1.11, 4.77)\) (right panel), during the fit. All the cases assume a positive signal of \(0\nu \beta \beta \) decays. The signal event number takes the expectation value of the choice \(m^{}_{3,\mathrm{true}} = 0~\mathrm{eV}\) and \(\left| { M}^{}_{0\nu }\right| ^{}_{\mathrm{true}} = 3\) with \(\rho ^{}_{\mathrm{true}} = 0^{\circ }\) (top panel), \(90^{\circ }\) (middle panel) or \(180^{\circ }\) (bottom panel)

Full size image

In Fig. 3, the sensitivity to \(|m^{}_{ee}|\) for the corresponding exposure with \(b=1~\mathrm{ton^{-1} \cdot year^{-1}}\) is indicated on the top axis. For a given exposure, this is obtained by converting \(T^{0\nu }_{1/2}\) in Eq. (7) with \(\left| { M}^{}_{0\nu }\right| = 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 LEGEND 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 Sect. 2 with the background index \(b=1~\mathrm{ton^{-1} \cdot year^{-1}}\) and the full detection efficiency \(\epsilon = 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 (\(\overline{\left| { M}^{}_{0\nu }\right| } \approx 3\)), the future projects such as LEGEND and nEXO can distinguish the Dirac and Majorana hypotheses with more than \(2\sigma \) 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\sigma \) 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\sigma \), 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.

3.2 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\nu }_{1/2}\) is mainly a function of the lightest neutrino mass \(m^{}_{3}\), the CP phases and the NME \(\left| { M}^{}_{0\nu }\right| \). 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\nu }_{1/2}\) on \(m^{}_{3}\) variation becomes weak. Hence, from the measurement of \(T^{0\nu }_{1/2}\) one can simultaneously constrain the Majorana phase \(\rho \) in the neutrino sector and \(\left| { M}^{}_{0\nu }\right| \) in the nuclear sector.

We start by exploring the sensitivity to \(\rho \) by making assumptions about the theoretical input of \(\left| { M}^{}_{0\nu }\right| \). In principle, we will need a positive \(0\nu \beta \beta \) decay signal to set bounds on \(\rho \), since a null signal can always be interpreted as a sign of Dirac neutrinos. Assuming a nominal event number \(N^{\mathrm{exp}}_{\mathrm{tot}}\) has been observed, we can set constraints on \(\rho \) as follows. First, we scan over all the model parameters and calculate the corresponding likelihood \({{\mathcal {L}}} (\theta ^{}_{13},\theta ^{}_{12},\varDelta m^{2}_{\mathrm{sol}},\varDelta m^2_{\mathrm{atm}}, m^{}_{3}, \rho , \sigma , \left| { M}^{}_{0\nu }\right| )\) with Eq. (4). Second, for every \(\rho \) we marginalize over all the other parameters to obtain the maximum \({\mathcal {L}}\) (or the minimum \(-2 \ln {{\mathcal {L}}}\)). Finally, we obtain the likelihood profile \(-2 \varDelta \ln {{\mathcal {L}}}\) as a function of \(\rho \) by subtracting its global minimum. The likelihood profiles of \(\rho \) with various nominal inputs are given in Fig. 4 for the experimental exposure \(\xi = 20~\mathrm{ton \cdot year}\). Four benchmark choices are illustrated: assuming \(N^{\mathrm{exp}}_{\mathrm{tot}}-B\) to be vanishing (gray curves) or to be the expectation values of \(m^{}_{3,\mathrm{true}} = 0~\mathrm{eV}\) and \(\left| { M}^{}_{0\nu }\right| ^{}_{\mathrm{true}} = 3\) with \(\rho ^{}_{\mathrm{true}} = \) \(0^{\circ }\) (dashed red curves), \(90^{\circ }\) (dotted pink curves) and \(180^{\circ }\) (solid red curves), respectively. The \(3\sigma \) C. L. range of \(\rho \) can be obtained by setting \(-2 \ln {{\mathcal {L}}^{}_{\mathrm{M}}} = 9\). In the absence of any signal events, one can still constrain \(\rho \) 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}\mathrm{Xe}\) NME, i.e., \(\left| { M}^{}_{0\nu }\right| = 3\) with a negligible uncertainty (left panel) and \(\left| { M}^{}_{0\nu }\right| \in (1.11, 4.77)\) (right panel). We find that the inclusion of the NME uncertainty will significantly dilute the sensitivity to \(\rho \) due to the understandable degeneracy. This result further encourages the efforts from the nuclear community to push forward the plausible calculation of \(\left| { M}^{}_{0\nu }\right| \).

Fig. 6
figure 6

The likelihood profile for the NME \(\left| { M}^{}_{0\nu }\right| \), in the presence of a null signal (gray curve) and a positive signal assuming \(m^{}_{3,\mathrm{true}}=0~\mathrm{eV}\) and \(\rho ^{}_{\mathrm{true}} = 90^{\circ }\) with \(\left| { M}^{}_{0\nu }\right| ^{}_{\mathrm{true}} = 1.11\) (dashed blue curve), 3 (dotted blue curve) or 4.77 (solid blue curve). The experimental configuration is taken as \(\xi = 20~\mathrm{ton \cdot year}\) and \(b = 1~\mathrm{ton^{-1} \cdot year^{-1}}\)

Full size image
Fig. 7
figure 7

The \(3\sigma \) sensitivities to the NME in terms of the \(^{136}\)Xe exposure \(\xi \). 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,\mathrm{true}}=0~\mathrm{eV}\) and \(\rho ^{}_{\mathrm{true}} = 90^{\circ }\) with \(\left| { M}^{}_{0\nu }\right| ^{}_{\mathrm{true}} = 1.11\) (top panel), 3 (middle panel) or 4.77 (bottom panel)

Full size image

The dependence of the sensitivity on the nominal exposure \(\xi \) is shown in Fig. 5, where we give the \(3\sigma \) allowed range of \(\rho \) with a signal event number calculated from \(m^{}_{3,\mathrm{true}} = 0~\mathrm{eV}\) and \(\left| { M}^{}_{0\nu }\right| ^{}_{\mathrm{true}} = 3\) with \(\rho ^{}_{\mathrm{true}} = 0^{\circ }\) (top panel), \(90^{\circ }\) (middle panel) or \(180^{\circ }\) (bottom panel). In the left panel, the NME is assumed to be known as \(\left| { M}^{}_{0\nu }\right| = 3\), and next-generation \(0\nu \beta \beta \) decay experiments will be able to rule out a significant amount of parameter space of the Majorana CP phase \(\rho \). In comparison, for the right panel we observe that the lack of knowledge of \(\left| { M}^{}_{0\nu }\right| \) significantly reduces the constraining power to \(\rho \), even with an exposure of \(100~\mathrm{ton \cdot year}\) and an ultra-low background of \(b=1~\mathrm{ton^{-1} \cdot year^{-1}}\).

The measurement of Majorana CP phases is of great theoretical importance. For instance, the \(\mu {-}\tau \) 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 (\(\rho , \sigma = \) \( 0^\circ \) or \( 180^\circ \)) [21, 99]. From the left panel of Fig. 5, we notice that if one of the solution is true (e.g., \( \rho _{\mathrm{true}} = 0^\circ \)), the remaining one (\( \rho _{\mathrm{true}} = 180^\circ \)) can be ruled out at 3\( \sigma \) C. L. with an exposure of only \(1~\mathrm{ton \cdot year}\). 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.

3.3 Nuclear matrix element

The NME, as an imperative theoretical input, can be conversely “measured” by \(0\nu \beta \beta \) decay experiments. If IO turns out to be true, the \(0\nu \beta \beta \) 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\nu \beta \beta \) decays, which is, however, not the main concern of the present work.

Following a procedure similar to fitting the CP phase, in Fig. 6 we generate the likelihood profile as a function of \(\left| { M}^{}_{0\nu }\right| \). For demonstration, several benchmark choices are shown, including a null signal (gray curve) and a positive signal assuming \(m^{}_{3,\mathrm{true}}=0~\mathrm{eV}\) and \(\rho ^{}_{\mathrm{true}} = 90^{\circ }\) with \(\left| { M}^{}_{0\nu }\right| ^{}_{\mathrm{true}} = 1.11\) (dashed blue curve), 3 (dotted blue curve) or 4.77 (solid blue curve).

In Fig. 7, we present the \(3\sigma \) sensitivities to the NME with respect to the nominal exposure \(\xi \). The true value of NME is taken to be \(\left| { M}^{}_{0\nu }\right| ^{}_{\mathrm{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 \(\left| { M}^{}_{0\nu }\right| \) is 3, a \(3\sigma \) C. L. constraint \(0.7< \left| { M}^{}_{0\nu }\right| < 3.5\) can be obtained with an exposure of \(\xi = 100~\mathrm{ton\cdot year}\), which is slightly narrower than the current prediction \(\left| { M}^{}_{0\nu }\right| \in (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\nu \beta \beta \) decays.

4 Conclusion

The next-generation \(0\nu \beta \beta \) decay experiments are designed for a target sensitivity around \(|m^{}_{ee}| \approx 10~\mathrm{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 next-generation projects will expect no \(0\nu \beta \beta \) 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\nu \beta \beta \) decays if IO eventually wins. From a positive or negative observation of \(0\nu \beta \beta \) 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 \(\rho \) can be probed with a positive \(0\nu \beta \beta \) 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}\mathrm{Xe}\) for instance, if the actual NME is found to be \(\left| { M}^{}_{0\nu }\right| = 3\) with a negligible uncertainty, the experimental proposals like nEXO with a \(10~\mathrm{meV}\) sensitivity can distinguish Dirac and Majorana neutrinos with \(3\sigma \) C. L. or so, in the presence of a positive or null signal (as shown by Fig. 3). If \(\left| { M}^{}_{0\nu }\right| \) 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 \(\rho \) (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\nu \beta \beta \) decay signal in the future (see Figs. 6, 7).