Constraining Majorana CP Phase in Precision Era of Cosmology and Double Beta Decay Experiment

We show that precision measurement of (1) sum of neutrino masses by cosmological observation and (2) lifetime of neutrinoless double beta decay in ton-scale experiments, with supplementary use of (3) effective mass measured in single beta decay experiment, would allow us to obtain information on the Majorana phase of neutrinos. To quantify the sensitivity to the phase we use the CP exclusion fraction, a fraction of the CP phase parameter space that can be excluded for a given set of assumed input parameters, a global measure for CP violation. We illustrate the sensitivity under varying assumptions, from modest to optimistic ones, on experimental errors and theoretical uncertainty of nuclear matrix elements. Assuming that the latter can be reduced to a factor of ~1.5 we find that one of the two Majorana phases (denoted as alpha_(21)) can be constrained by excluding ~10-40% of the phase space at 2sigma CL even with the modest choice of experimental error for the lowest neutrino mass of 0.1 eV. The characteristic features of the sensitivity to alpha_(21), such as dependences on the true values of alpha_(21), are addressed.


I. INTRODUCTION
After all the angles in the lepton mixing matrix [1] and the mass squared differences of neutrinos are determined by the atmospheric, the solar, the reactor, as well as the accelerator neutrino oscillation experiments, we have entered into a new phase of neutrino physics [2][3][4]. Of course, we are still with the important unknowns in the lepton sector, CP violating phase of Kobayashi-Maskawa type [5] and the neutrino mass hierarchy. Yet, we do have relatively clearer view of how these questions are to be settled, either through the ongoing searches, or by new measurement to be carried out by the powerful apparatus described in various experimental proposals to date.
It is quite possible that even after these unknown quantities are measured by experiments there will still remain the questions of somewhat different category: What are the absolute masses of neutrinos? Are neutrinos Dirac or Majorana particles? If Majorana, which values of the Majorana CP phases [6][7][8] are chosen by nature? The first one, knowing the absolute mass scale, is crucial to understand physics behind the origin of neutrino mass, and the latter two are likely the key to understand baryon number asymmetry in the universe [9].
In this paper, we focus on detectability of the Majorana CP phase in the light of informations to be obtained by cosmological observation and neutrino experiments in the coming precision era. In particular, in the near future the double beta decay experiments aim to cover entire region of parameters allowed by the inverted mass hierarchy. (See Fig. 1.) We show that given the future perspective, when combined with precision measurement in cosmology, one can indeed expect the sensitivity, though limited one, to one of the phases (denoted as α 21 , see section II A). We hope that this result can trigger renewed interests in the problem of detectability of the Majorana phase despite the dominant pessimism which prevailed in the past decades. See, however, e.g., Refs. [10][11][12][13][14][15][16][17][18][19][20][21] for the foregoing efforts which examined the sensitivity to the Majorana phase.
We feel that now is the particularly right time to revisit these issues because of the following reasons. Firstly, cosmological observations entered into the precision era so that their data became sensitive to the absolute masses of neutrinos better than the current bounds from laboratory experiments, as seen most notably by the results obtained by the Planck collaboration [22]. See also Ref. [23] for overview and many relevant references. It makes us possible to employ a new strategy that the neutrino mass scale can be constrained mainly by cosmology and we can use data of neutrinoless double beta decay experiments, with supplementary use of single beta decay data, to constrain the Majorana phases of neutrinos. Even more interestingly a discrepancy between the cosmic microwave background (CMB) and lensing observations about cluster correlations would prefer massive neutrinos with mass of a few tenth of eV [24]. If confirmed, it would make the basis of our analysis more robust.
Secondly, as mentioned above, several new neutrinoless double beta decay experiments are in line to measure or constrain, the ee element of neutrino mass matrix in the flavor basis, so called the effective Majorana mass, to a few tens of meV [25][26][27][28][29][30][31][32]. For reviews and description of the other projects, see e.g., Refs. [33][34][35]. Needless to say, they are the prime candidates among various experiments that can measure or constrain the Majorana phases. We should also note that intensive efforts have been devoted in order to improve the calculations of the nuclear matrix element (NME) which is crucial to determine the value of effective Majorana mass from the measured half life time in double beta decay experiments [36][37][38][39]. See, also Ref. [35] for review of the recent progress in the calculations of NME.
Thirdly, the uncertainties in the relevant lepton mixing angles get decreased dramatically recently. The error of sin 2 θ 12 is now only ∼ 4% [40], and moreover it will be decreased to a sub-% level by future medium-baseline reactor experiments [41,42]. θ 13 , the mixing angle whose value was unknown until intriguing indication from T2K [43] by seeing ν e appearance events, is now precisely measured [44,45], and it will soon reach to a even greater precision, the error of sin 2 θ 13 would become ∼ 5% level. See Fig. 1 to know how small is the effect of the current uncertainties of the mixing parameters onto the allowed regions by double and single beta decay experiments. These developments are important to tighten up the constraint on Majorana CP phases and to place it on a firmer ground.
In addition, the single beta decay experiment KATRIN [46] (see also [47]) will constrain neutrino masses to a level of 0.2 eV at 90 % CL, or observe the effects consistent with masses of this order, in a manner quite complementary to cosmological observations. While enjoying the above new inputs, our analysis is nothing but a continuation of the numerous similar analyses attempted by the many authors, from which we are certainly benefited.
To display the sensitivities to the Majorana CP phases in a more global way, we use in this paper the CP exclusion fraction f CPX , which is defined as a fraction of the CP phase space that can be excluded for a given set of input parameters [48] 1 . In addition to f CPX plots, we also use the more conventional one, the allowed regions of phases for given input parameters and experimental setup which demonstrate how well the input parameters can be reproduced. The CP exclusion fraction is a useful tool particularly in the initial era of search for effects of the Majorana CP phase, in which even the partial exclusion of allowed range of the phase would be very valuable. On the other hand, f CPX has intimate connection to the error of CP phase, an appropriate measure in precision measurement era, as discussed in [48]. It suggests that the CP exclusion fraction can be a universally usable measure for CP sensitivity.
In the following two sections, II and III, we briefly review the basic features of the observables used and describe simple analytic estimation of the effects of Majorana phases. The expert readers must go directly to the results in sections VI and VII, starting from describing the analysis method in section IV.

II. ASSUMPTIONS AND OBSERVABLES USED IN THE ANALYSIS
In this work we assume that the neutrinos are of Majorana type, and the Majorana mass is the unique source of neutrinoless double beta decay. We also assume the standard three-flavor mixing scheme of neutrinos.
In our analysis we consider three observables [51]: (1) the effective neutrino mass (denoted as m 0νββ ) measured by neutrinoless double beta decay experiments, (2) the sum of the neutrino masses (denoted as Σ) which can be determined by cosmological observations, and (3) the effective neutrino mass (denoted as m β ) measured by single beta decay experiments. They are all sensitive to the absolute neutrino mass scale. The sensitivity to the Majorana CP violating phase, which is of our main concern in this paper, primarily comes from m 0νββ in (1). But we will see that the other two, in particular Σ in (2), greatly help enhancing the sensitivity. Let us briefly describe their current status.

A. Effective mass for neutrinoless double beta decay
The rate of neutrinoless double beta decay (hereafter denoted as 0νββ decay), or the inverse of the half life time of the 0νββ decay, T 0ν 1/2 , can be written as where G 0ν is the kinematic phase space factor, M (0ν) is the nuclear matrix element (NME) corresponding to the 0νββ transition and m e is the electron mass. The effective neutrino mass m 0νββ measured by 0νββ decay experiments is given by where m i (i = 1, 2, 3) are neutrino masses and U ek (k = 1, 2, 3) are the element of electron-row of the Maki-Nakagawa-Sakata (MNS) neutrino mixing matrix [1], α 21 and α 31 are the Majorana CP phases and c ij ≡ cos θ ij , s ij ≡ sin θ ij . Throughout this work, as we did in the second equality in Eq.
(2), we use the standard parameterization of the MNS matrix including the Majorana phases [52] 2 . We note that while G 0ν can be calculated very accurately as performed, e.g., in Refs. [53,54], M (0ν) is not known very well, the fact which should be taken into account as we will see below.
So far there is no convincing evidence for 0νββ decay or non-zero value of m 0νββ . Though it was claimed by the authors of Ref. [55] that they saw a positive signal at more than 4σ CL which is consistent with m 0νββ = 0.2 − 0.6 eV, it does not appear to become a consensus in the community. The result is not confirmed by GERDA [27] which used the same nuclei 76 Ge. Furthermore, a combined analysis of EXO [25] and KamLAND-Zen [26] data, both using 136 Xe, results in the upper limit of the effective mass which appears to exclude the above claim at a confidence level higher than 97.5% [26]. Fortunately, the sensitivity to m 0νββ in the ongoing and the next generation 0νββ decay experiments will go down to the region of a few tens of meV [33]. These experiments would eventually cover an entire region expected for the inverted mass hierarchy (see Fig. 1), and the controversy will certainly be settled on the way toward the goal. The largest uncertainty in translating the bounds on T 0ν 1/2 to that of m 0νββ comes from the uncertainty of the NME, M (0ν) , which produces typically a factor of ∼ 2 or more differences in terms of m 0νββ . Despite that a lot of efforts have been devoted, there are still some considerable differences among NME values obtained by different authors and by different calculation methods. Generally speaking, the difference between different authors is not very large if the same calculation technique is used. The best hope would be to normalize the coupling strength by using the 2ν decay by which the difference may go down to 30% level, as proposed within the framework of quasi-particle random phase approximation (QRPA) method [36,37]. But, notable difference persists to among results calculated by different techniques 3 .

B. Sum of neutrino masses
For our purpose, in addition to m 0νββ , it is crucial to take into account the sum of the neutrino masses, Σ, where the first limit was obtained by combining the Planck data with that coming from WMAP, CMB polarization data, and the second limit was obtained by adding also the data of baryon acoustic oscillation (BAO) (see Ref. [22] for details). In deriving these bounds, deviation from the flat universe was allowed.
However, more recently, some analyses are in favour of non-zero Σ which could be attributed to the sub-eV masses of the active or sterile neutrinos [24,63,64]. In Ref. [24] it was discussed that a discrepancy between CMB and lensing observations about cluster correlations prefers massive neutrinos leading to the preferred non-zero value as Σ = 0.32 ± 0.081 eV, which implies that m 1 ∼ m 2 ∼ m 3 ∼ 0.1 eV in the standard three flavour scheme. We expect that some future cosmological data either from galaxy surveys or weak lensing data by the Euclid satellite, with Planck constraints, may lead to the sensitivity to Σ 0.01-0.05 eV [65][66][67][68].
Despite the success in tightening up the neutrino mass bound by cosmological data we should keep in mind that the analysis has done (or needs to be done) with certain prior or within a cosmological model such as the standard ΛCDM model.

C. Effective mass for beta decay
In addition to the observables mentioned above, the absolute neutrino mass scale can also be probed by the single beta decay experiment, which measures the effective electron neutrino mass, m β , defined as providing us information independent of cosmology and 0νββ decay experiment. So far there is no evidence of non-zero value of m β and only the upper bounds, have been obtained by Mainz [69] and Troitsk [70] collaborations. In the near future, we expect that this bound can be improved by one order of magnitude by KATRIN [46] which is expected to reach the sensitivity of m β ∼ 0.2 eV at 90% CL.
D. Allowed regions of m 0νββ and m β as a function of lightest neutrino mass In this paper we use the lightest neutrino mass (rather than Σ), denoted as m 0 , as the relevant parameter (observable) to be determined by the experiments. Namely, m 0 ≡ m 1 and m 0 ≡ m 3 for the normal and the inverted mass hierarchies, respectively. All the other masses and their differences can be calculated by using the relations m 1 ≡ m 0 , m 2 = m 2 0 + ∆m 2 21 , m 3 = m 2 0 + ∆m 2 21 + ∆m 2 32 (normal mass hierarchy), where ∆ 2 m ij ≡ m 2 i − m 2 j (i, j = 1, 2, 3). Notice that ∆m 2 32 in Eq. (9) for the case of inverted mass hierarchy is a negative quantity.
While our results are insensitive to small variation of the mixing parameters we use, for definiteness, the values obtained in [71] will be used as our reference parameters throughout this work: ∆m 2 21 = 7.54 × 10 −5 eV 2 , sin 2 θ 12 = 0.308 or sin 2 2θ 12 = 0.853, for both mass hierarchies whereas ∆m 2 32 = 2.40 (−2.44) × 10 −3 eV 2 , sin 2 θ 13 = 0.0234 (0.0239) or sin 2 2θ 13 = 0.0914 (0.0933), (11) for the normal (inverted) mass hierarchy 4 . For later convenience we define ∆m 2 atm ≡ |∆m 2 32 |. Using these values of the parameters m 0νββ at its current status can be summarized as for the normal (inverted) mass hierarchy. In the case of normal (inverted) mass hierarchy the ranges are shown by green (blue) colour. In the both panels, the light (dark) coloured regions are computed by taking into account (without taking account) the current 1 σ uncertainties of the relevant mixing parameters. Also shown are the limits on m 0νββ coming from KamLAND-Zen and EXO [26] (by the light brown band and arrow), the bounds on m0 obtained by Planck collaboration [22] (by the magenta and the brown dotted lines), and also the upper limit on m β coming from Mainz [69] and Troitsk [70] experiments (by the black dotted line).
In Fig. 1, we show the allowed ranges of m 0νββ and m β as a function of m 0 for both the normal and the inverted mass hierarchies. The allowed ranges are calculated by the two ways, without and with the current 1σ uncertainties of the mixing parameters, which are indicated, respectively, by the dark and light colours, blue and green for the inverted and normal hierarchy. Notice that the effect of 1σ uncertainties of the mixing parameters is quite small by now. In contrast, additional variation over the Majorana phases α 21 and α 31 gives much larger impact on m 0νββ , not only producing sizeable width but also creating a down-going branch in region 10 −3 eV < ∼ m 0 < ∼ 10 −2 eV for the case of the normal mass hierarchy due to the strong cancellation of the three mass terms in Eq. (2).
The difference between the normal and the inverted mass hierarchies becomes clearly visible only in small m 0 region, typically m 0 < ∼ 0.1eV in the both plots. Therefore, onset to the almost degenerate mass regime, which is formally defined as m 0 ∆m 2 atm 0.05 eV, already takes place for m 0νββ and m β at around m 0 > ∼ 0.1eV. In the same plots, we also show the limits on m 0νββ coming from KamLAND-Zen and EXO [26] (by the light brown band and arrow), the bounds on m 0 obtained by Planck collaboration [22] (shown by the magenta and the brown dotted lines), and also the upper limit on m β coming from Mainz [69] and Troitsk [70] experiments (shown by the black dotted line). We note that the upper limit on m 0νββ provided by KamLAND and EXO is expressed by the band (not by line) because of the NME uncertainty.

III. ANALYTIC ESTIMATE OF THE EFFECT OF MAJORANA PHASE
In this section, we discuss qualitative features of effects of the Majorana phase onto the 0νββ decay observable m 0νββ . It is partly for pedagogical purpose but we hope that some new points presented serve for better understanding the results to be obtained in later sections. In particular, we aim at illuminating how and where the best sensitivity to the Majorana phase can be expected. After briefly mentioning a general properties of m 0νββ we will focus on this problem.

A. Reduction of Majorana phase space
We give here a general properties of m 0νββ in Eq. (2). Since the phases α 21 and α 31 enter into m 0νββ in the form of either one of cos α 21 , cos α 31 , and cos(α 31 − α 21 ) m 0νββ has the following reflection symmetry Notice that it is the simultaneous transformation of two phase variables and each one of α 21 or α 31 does not obey individually the symmetry under α ij → 2π − α ij .

B. The case of almost degenerate mass spectra
The best place to discuss clearly how and where the best sensitivity to the Majorana phase can be expected is in the almost degenerate regime, It means that the ratio of the maximum to the minimum value of effective mass changes by a factor of m max 0νββ /m min 0νββ 2.6 as the Majorana phase is varied. Therefore, if one can determine m 0νββ from experiment with errors smaller than the factor 2.6, after taking into account the uncertainty coming from the NME, it is in principle possible to constrain the Majorana phase α 21 . Then, the question is: In what regions of α 31 the best sensitivity to α 21 can be expected? We are primarily interested in the sensitivity to α 21 because the sensitivity to α 31 itself is quite low in our setting, due to the small value of θ 13 , as we will see in later sections. As the sensitivity to the phase arises in places where the NME uncertainty cannot obscure the effect of phases the best place is at the maximal and the minimal value of m 0νββ . One can show that m 0νββ in Eq. (14) is nearly periodic function of α 21 and α 31 , and has the following properties in regions 0 ≤ α ≤ π: (1) The function inside parenthesis in Eq. (14) is an extremely slowly decreasing function of α 31 for in the entire region of α 21 . (2) It is a decreasing function of α 21 in the entire region of α 31 , and varies from 1 to 0.15 as α 31 is varied from 0 to π. Therefore, it has the maximum (minimum) at α 21 = α 31 = 0 (α 21 = α 31 = π). Thus, the best sensitivity to α 21 is obtained at either α 31 = 0 or π. On the contrary, the worst sensitivity is expected at around the mean value of m 0νββ between their maximum and minimum values. Typically, they are at around α 21 2π/3 for the whole range of α 31 . These features will be confirmed explicitly by quantitative analysis in later sections. C. Hierarchical case; inverted mass spectrum Now let us consider the case of m 0 < ∆m 2 21 < ∼ 0.01 eV for the inverted mass ordering. In this case, m 1 ≈ m 2 ∼ ∆m 2 atm m 3 ≡ m 0 . This is the regime corresponding to the horizontal band for the inverted mass spectra in Fig. 1. It is interesting to observe that essentially the same treatment can go through for m 0νββ in these two entirely different mass spectra. In the case of strongly hierarchical inverted mass spectrum the m 3 term in Eq. (12) is suppressed than other two terms by a factor of m 3 /m 2 [∆m 2 21 /∆m 2 atm ] 1/2 0.18 and s 2 13 = 0.0234 (0.0239) for the normal (inverted) hierarchy. In the case of almost degenerate spectrum the suppression of m 3 term is only due to the latter factor.
Under the approximation of ignoring the m 3 term, m 0νββ is expressed as in the degenerate case, (14) without s 2 13 term, but replacing m 0 by ∆m 2 atm , Since the dependence of m 0νββ on the Majorana phase is essentially the same as in the case of the degenerate mass regime, the same results hold for the sensitivity to α 21 and the location of highest sensitivity. Under our assumption of the precision of δ(m 0νββ ) ≡ σ 0νββ = 0.01 eV (see section IV), however, the sensitivity to α 21 is quite limited in this regime when it is combined with the NME uncertainty. It will be demonstrated in the analysis in section VII.
Or numerically, a factor of ∼ 4 variation. Again, under our assumption of the precision of σ 0νββ = 0.01 eV it would be very difficult to probe the value of α 31 in this regime, as we will see in section VII.

IV. ANALYSIS METHOD
Throughout this paper, unless otherwise stated, we use the mixing parameters described in Eqs. (10) and (11). We do not consider the uncertainties of the mixing parameters apart from the currently unknown mass hierarchy, because the impact is rather small which is expected to become even smaller in the near future. This is because we expect significant improvement in the accuracies of these parameters by future oscillation experiments, in particular, by the medium baseline (∼ 50 − 60 km) reactor experiments such as JUNO [41] and RENO-50 [42] as demonstrated, e.g., in Refs. [72][73][74][75].

A. Assumption on experimental errors
As discussed in the previous sections, we use the three observables, sum of neutrino masses Σ, and the effective mass parameters m β and m 0νββ which are measured, respectively, in single and double beta decay experiments. Let us assume that they are measured with some uncertainties around the central values, where σ Σ , σ β and σ 0νββ denote the corresponding 1 σ uncertainties. We take the following values for the uncertainties, which will be used as the reference values in our analysis.
Here is some reasonings for our choices of the uncertainties in Eq. (23). It appears that future galaxy surveys [65] and weak lensing [66] both lead, under the Planck constraint, the sensitivity to Σ to the level of 0.03 − 0.05 eV. According to Ref. [68], in the most optimistic case, the sensitivity could even lead to ∼ 0.01 eV. Therefore, our choice would be in a conservative side. For KATRIN [47] they quote the error 0.025 eV 2 for m 2 β , which may be translated into the error of 0.063 eV for m β for the case where true value of m 0 = 0.2 eV, which justifies our choice in Eq. (23).
To have a reasonable estimate of the uncertainty in m 0νββ we need a detailed discussion because we must address the question of how the experimental uncertainty on T 0ν 1/2 is translated into that of m 0νββ , and the treatment heavily depend on to what extent the experimental backgrounds can be suppressed. We present such a discussion in Appendix A, based, e.g. on Refs. [35,76]. Under the assumption of dominance of statistical error, we will argue there that the uncertainty of 0.01 eV may be reachable in future experiments that are upgraded to a ton scale. Such experiments will be able to cover the entire allowed range corresponding to the inverted mass hierarchy.

B. Analysis procedure
As we emphasized in section II A the error of m 0νββ due to uncertainties of the NME would be a major obstacle to extract informations of the Majorana phases from 0νββ decay experiments. Therefore, its treatment is of key importance in our analysis to constrain the Majorana phases. We follow the procedure similar to the one used in Ref. [19].
Let us assume that the unknown true value of NME, denoted as M (0ν) , is restricted to some range as where M (0ν) min and M (0ν) max denote the minimum and maximum possible values of NME, respectively. Then, we introduce the NME uncertainty parameter r NME ≡ M min can be written as Notice that r NME = 1 implies no uncertainty on NME.
We take into account of the uncertainties of NME via the following way. Using Σ, m β , and m 0νββ in Eqs. (20) -(22), we define the χ 2 function as follows, where will be varied in the interval of [1/ √ r NME , √ r NME ], and the minimization will be done with respect to m 0 , α 21 , α 31 and ξ. In this way, the NME uncertainty is treated in a similar way as the systematic error. (It is very similar to the treatment of a normalization parameter to take into account of the flux uncertainty in reactor neutrino experiments.) In this work, we consider three values of the NME uncertainty parameter, r NME = 1.3, 1.5 and 2.0. (Exceptional use of an extreme value, r NME = 1.1, is made at the end of section VII.) We note that, strictly speaking, errors associated to the three observables we consider in this work are not Gaussian, in particular, the one comes from cosmology. It would be proper to do analysis based on Monte Carlo simulations, as performed, e.g., in Refs. [80,81]. However, we believe that the procedure we use produces approximately correct results, which are sufficient for our purpose. Currently, for a given nuclei, the typical ratio of the maximum to minimum NME calculated by different techniques are about the factor of 2. So we consider r NME = 2.0 as the current situation. We assume that in the future, theoretical uncertainty on the NME could be reduced to 50% or even to 30%, which may be possible if the main source/reason of the discrepancies of the NME values calculated by different techniques are identified. We expect that once the positive signal of 0νββ will be observed by using different isotopes, by comparing the real data from these different isotopes, it would be possible to check in a systematic way the validity of various theoretical calculations based on different models or methods, which may allow to reduce significantly the NME uncertainty.
Using the χ 2 function defined in Eq. (28), we will determine the allowed parameter space by imposing the condition, 3, 6.18 and 11.83 (1, 4 and 9) for 1, 2 and 3 σ CL for two (one) degree of freedom.

C. Plotting allowed regions and CP exclusion fraction
In the following two sections we will present our analysis results. In section VI we display the allowed regions for some combinations of the assumed true (input) values of m 0 , α 21 and α 31 . In section VII we present the contours of equal CP exclusion fraction f CPX . It is defined [48] as the fraction of values of α 21 ∈ [0, 2π] which can be excluded by the experiments at a given confidence level for each input point of the parameter space (m 0 , α 21 ). The f CPX contour plot can be defined for phase α 31 as well. It is thus a global measure which covers the entire input parameter space. For more precise definition of f CPX see Ref. [48].

V. HOW CAN THE MAJORANA PHASES BE CONSTRAINED?
One may ask why and how the Majorana phases can be constrained in the presence of NME uncertainty ξ in Eq. (27). Before presenting our main results let us make a short discussion and present some results of our analysis on this point because it illuminates how the NME uncertainty can be circumvented at least partly.
For clarity we focus on the case of almost degenerate mass spectrum discussed in section III B. In this case m 0νββ prepared by nature is related to the "observed" one, the inferred one with use of particular NME (using Eq. (14) without the terms which contains sin 2 2θ 13 ) as Since ξ enters into Eq. (30) in combination with m 0 , if we are completely ignorant about m 0 , there is no way to constrain the value of α 21 even in the case of no NME uncertainty (ξ = 1) no matter how accurately m 0νββ itself is determined. Therefore, in order to measure the Majorana phases, it is of crucial importance to determine m 0 from the other observation different from 0νββ decay 5 . We will show below that the major role is played by cosmological observation.
Assuming that m 0 is determined with good accuracy, change of m true 0νββ in Eq. (30) as α 21 is varied is of a factor of 2.6 (see section III B). The phase information can be gained only if the NME uncertainty ξ cannot compensate the variation. Thus, we expect that sensitivity to the Majorana phase is maximal at the true value of α 21 at 0 or π where m true 0νββ has extremal values in the case of degenerate mass spectrum. Let us examine relative contribution of each experiment in constraining the Majorana phase, α 21 . In Fig. 2  First we confirm that the contribution from the 0νββ decay alone can not determine or constrain at all the value of α 21 (see blue dotted lines in Fig. 2) which is true even if there is no NME uncertainty. Even if the result of beta decay is added to that of 0νββ, the restricting power of α 21 is very weak as we can see from the green dashed lines. However, once the cosmological observation of Σ is added ∆χ 2 jumps outside the region 0.6π − 1.4π of α 21 as shown by solid red curve in Fig. 2. Therefore, it is the precision measurement of Σ, of course together with 0νββ decay experiments, to produce nontrivial constraint on α 21 by circumventing the NME uncertainty. We note that at the bottom of ∆χ 2 in region α 21 = 0.75π − 1.25π, χ 2 min is very flat because of the NME uncertainty. This means that it is impossible to pin down, in the presence of the NME uncertainty, the value of the Majorana phase α 21 at the bottom no matter how accurately all the measurements are done.
We also observe that the sensitivity to the Majorana phase is higher for the inverted mass hierarchy. This is because, for a given value of m 0 , the first and second terms which contain m 1 and m 2 , respectively, in Eq. (2) are always larger in the case of inverted hierarchy than the ones in the normal one, as we can see from Eqs. (8) and (9). Assuming the difference between |∆m 2 32 | in both mass hierarchies is small, as is the case in Eq. (11), it makes ∆χ 2 for the inverted mass hierarchy larger, thereby making the inverted hierarchy case more sensitive to the change of α 21 . Let us also examine the contribution from each experiment to the determination of m 0 . In Fig. 3 we show ∆χ 2 obtained by marginalizing the other parameters as a function of the fitted value of m 0 for the case of inverted (left panel) and normal (right panel) mass hierarchies. Here, we present five (out of seven) possible combinations in the plot. Again the contribution from cosmology dominates the ∆χ 2 curve and the contribution from 0νββ decay is visible only in the region m 0 > ∼ 0.15 eV for the inverted mass hierarchy and m 0 < ∼ 0.03 eV and m 0 > ∼ 0.15 eV for the normal mass hierarchy.
We conclude that 0νββ experiments and cosmological observation work together in a very complementary fashion to produce nontrivial constrains on the Majorana phase. In the case of inverted mass hierarchy with hierarchical spectrum, m 0 in Eq. (30) is to be replaced by ∆m 2 atm , and the same treatment goes through.

VI. ANALYSIS RESULTS I: ALLOWED REGIONS OF m0, α21 AND α31
In this section, we study the sensitivity to the CP phase through the conventional plots of allowed regions, which indicate how accurately one can reproduce the assumed input parameters. Let us first consider, for the purpose of illustration, the case where the true values of relevant fundamental parameters are, m 0 = 0.1 (Σ 0.32 eV) and α 21 = α 31 = π. In Fig. 4 the allowed regions projected into the planes of (a) α 21 − m 0 (b) α 21 − α 31 and (c) m 0 − α 31 are drawn. The NME uncertainty is taken as r NME = 1.5. In the case of inverted mass hierarchy they are depicted by the filled colours, yellow, red, and light blue, corresponding, respectively, to 1, 2 and 3 σ CL for 2 degree of freedom (DOF). While for the normal mass hierarchy they are drawn by the black dotted (1σ), dashed (2σ) and the solid (3σ) curves.
We first notice that for both mass hierarchies, it is possible to restrict α 21 in some interval at 2σ CL (for 2 DOF) but no constraint on α 21 can be placed at 3σ CL as we can see from the panels (a) and (b) of Fig. 4. On the other hand, as expected, it is not possible to constrain at all the phase α 31 due to the small value of θ 13 , see Eq. (12). We have checked that this is true even if we consider different true values of α 31 , where the dependence of the allowed regions on α 31 is always very mild. We also note that the allowed regions in Fig. 4(b) are symmetric with respect to α 21 = π in a good approximation. It stems from the simultaneous reflection symmetry α 21 → 2π − α 21 and α 31 → 2π − α 31 in Eq. (13) and the weak dependence of ∆χ 2 on α 31 .
In Fig. 5 we show the same plots as in Fig. 4 but for different value of m 0 = 0.2 eV (Σ 0.61 eV). In this case, as we can see in Fig. 4(a) and (b) that we can constrain α 21 significantly better than the case of m 0 = 0.1 eV for the both mass hierarchies. We note that indeed the difference between the normal and inverted mass hierarchy is very small for m 0 = 0.2 eV, which is expected as it is in the almost degenerate regime.
Since the dependence on α 31 is quite mild, from now on we will show only the allowed regions projected into the plane of α 21 − m 0 , which is the most interesting combination for our purpose. In Fig. 6 we show the regions allowed  Fig. 6 are for the cases where the true values of α 21 = 0, α 21 = π/2, and α 21 = π, respectively. The input values are indicated by the asterisks. The allowed regions for the inverted mass hierarchy for r NME = 1.3, 1.5 and 2.0 are indicated by filled colors of yellow, red, and light blue, respectively, whereas that for normal mass hierarchy are indicated by the black dotted, dashed and solid curves, respectively. We first note that due to the marginalization with respect to α 31 , the allowed regions in Fig. 6 show the reflection symmetry around α 21 = π. Generally speaking the sensitivity to α 21 is higher for larger m 0 . The exception is the case of the inverted hierarchy with m 0 = 0 eV and the input α 21 = π, where the allowed contour is slightly smaller than that of the same hierarchy with m 0 = 0.05 eV. Among the input values of α 21 considered for Fig. 6, the sensitivity to α 21 is highest (though modest) for α 21 (true) = π and lowest for α 21 (true) = π/2, which is consistent with our expectation described in section III B. Even in the highest sensitivity case examined (top row) the allowed region cannot be smaller than a third of the range of α 21 . On the other hand, the sensitivity to α 21 still persists even for m 0 less than 0.05 eV. We note, however, if the allowed (or excluded) range of α 21 is evaluated for 1 DOF, which is reasonable if we are interested in restricting only the single parameter α 21 , as we will do in the next section, the sensitivity would become better.
As we can see in the first and the second rows of Fig. 6, the difference between normal and inverted hierarchies is quite small if the parameter m 0 is well within the degenerate regime, m 0 = 0.2 eV, but it starts to develop when m 0 is increased to 0.1 eV, as we can see in the first and second rows of Fig. 6. For m 0 smaller than 0.05 eV, the difference between the hierarchies become more visible, as can be seen in third and fourth rows of Fig. 6.
Let us finally comments also on the determination of m 0 . We first observe that the allowed ranges of m 0 for a given input value of m 0 do not depend essentially on the assumed input values of α 21 in agreement with our expectation as m 0 is mainly determined by cosmological observation coming from the constraint on Σ. We further note that if m 0 (true) = 0.05 eV and the true mass hierarchy is inverted, it is not possible to exclude the case of m 0 = 0 at 2σ CL even for r NME = 1.3. This is because if the hierarchy is inverted, the case of m 0 (true) = 0.05 eV, or equivalently Σ(true) = 0.19 eV, is compatible with m 0 (fit) = 0 since in this case the sum of the other 2 masses give m 1 + m 2 0.14 eV. On the other hand, for the normal hierarchy, m 3 0.05 eV alone is not sufficient to mimic the case of m 0 (true) = 0.05 eV, which is equivalent to Σ(true) 0.17 eV for the cases where α 21 = 0 and π/2 (see the panels (g) and (h) of Fig. 6).

VII. ANALYSIS RESULTS II: CP EXCLUSION FRACTION
In this section, we show the results on the CP exclusion fraction f CPX , which is defined as a fraction of the CP phase parameter space that can be excluded for a given set of input (true) parameter set of (m 0 , α 21 , α 31 ) and the NME uncertainty. As discussed in Ref. [48], this is a global measure of CP violation, which we believe particularly useful in the stage where the sensitivity to the CP phase is not so high.
We show in Fig. 7 the iso-contours of the CP exclusion fraction determined at 2σ (95.45%) CL in the plane spanned by the true values of α 21 and the lightest neutrino mass m 0 for r NME = 2.0, 1.5, and 1.3 from left to right, and the cases where the true values of α 31 is 0 and π shown, respectively, in the upper and lower panels 6 . Since the iso-contours of f CPX are symmetric under reflection around α 21 = π after marginalizing over α 31 due to the symmetry in Eq. (13), from now on, we show the results only for the range of 0 ≤ α 21 ≤ π. In Fig. 7 the iso-contours of f CPX for the inverted and the normal mass hierarchies are depicted, respectively, by the solid and the dashed lines, from 0.1 to 0.7 (with the step size of 0.1), using different colors as indicated in the legend. Notice that larger the values of f CPX , higher the sensitivity to α 21 because larger fraction of the phase space can be excluded. We note that for each point on the input parameter space, all the parameters, m 0 , α 21 , α 31 and ξ were varied in doing the fit.
In each panels in Fig. 7, the sensitivity is highest in region around α 21 π and next highest in region of α 21 near 0. The worst sensitivity is obtained at around α 21 ∼ 2π/3, in agreement with the qualitative discussion given in section III. As expected the difference between the normal and the inverted mass hierarchies are small in the degenerate regime, m 0 > ∼ 0.1 eV where the solid (inverted mass hierarchy) and dashed lines (normal mass hierarchy) overlap rather well.
It is also notable that the case of α 31 (true)= π (upper panels of Fig. 7) provides somewhat better sensitivity at around α 21 = π than in the case of α 31 (true)= 0 (lower panels) whereas this behaviour is contrary for the value around α 21 = 0 for which the results in the lower panels provides somewhat better sensitivity (see also the footnote 6). However, the differences are not so large, or the α 31 dependence is quite weak, which may be anticipated from Eq. (14) due to small θ 13 . Since we will not be able to know the true value of α 31 in our current set up (no restriction on α 31 ), in reality, we must assume that the real value of f CPX is somewhere between the results shown in the upper and lower panels of Fig. 7.
It is interesting to see that the sensitivity to α 21 strongly depends upon r NME . In particular the improvement of the sensitivity from r NME = 2.0 to r NME = 1.5 is remarkable. At m 0 = 0.15 eV, for example, while f CPX ≥ 0.3 region spans only 0.2 of α 21 space for r NME = 2.0, it jumps to 50% coverage for r NME = 1.5 for both cases of α 31 = 0 and π. It is also notable that in the inverted mass hierarchy case the sensitivity region reaches to the range 2, the tendency more and more visible for smaller values of r NME . Therefore, reducing the NME uncertainty is of crucial importance to have severer constrains (or to observe) the Majorana phase. Though it is encouraging to see some sensitivities to the Majorana phase thanks to precision measurement of Σ and 0νββ decay rate, we must remark that it is still at a relatively low level. For example, in the case of m 0 0.1 eV, which is preferred by some cosmological analyses, the sensitivity to α 21 is up to f CPX 0.1 − 0.3 level even in the most optimistic case of r NME = 1.3 except for the region close to α 21 = π (see panels (c) and (f) of Fig. 7). Therefore, we examine below two possibilities toward further improvement of the sensitivity to α 21 : (1) Higher precision measurement of neutrino masses in Cosmology as good as σ Σ = 0.02 eV, and (2) Revolutionary new technology of computing the NME of 0νββ decay which could lead to uncertainty of the level r NME = 1.1. As mentioned before, once the positive signal of 0νββ will be observed by using different isotopes, we hope that theoretical NME calculations can be "calibrated" to some extent using the real data. (We tentatively assume that 10% level uncertainty is unavoidable.) It is discussed that the former is within reach in the light of future cosmological observation as mentioned in section II B [65][66][67]. On the other hand, it remains to be seen how (2) can be realized.
We first consider the possibility (1). In the upper panels in Fig. 8 we show the same plots as shown in the upper panels of Fig. 7, the case where the true value of α 21 = 0, but generated only by changing the value of σ Σ from 0.05 eV to 0.02 eV to compute f CPX . By comparing results shown in the upper panels of Fig. 8 to the case of σ Σ = 0.05 eV shown in upper panels of Fig. 7, at 2 σ CL, in the parameter regions corresponding to α 21 < ∼ 0.6π, the improvement (increase of f CPX ) with σ Σ = 0.02 eV is not so large, typically less than 0.1. However, for α 21 > ∼ 0.6π, we observe that there are some regions, in m 0 < ∼ 0.2 eV, where the improvement is more significant, leading to ∆f CPX (Σ) ≡ f CPX (Σ = 0.02eV) − f CPX (Σ = 0.05eV) ∼ 0.1-0.3. For instance, at around m 0 0.1 eV and α 21 0.8π, for r NME =1.5, ∆f CPX (Σ) is ∼ 0.2 or so (see Fig. 8(e) and text below).
It is not so easy to see clearly, just by comparing by eyes the upper panels of Figs. 7 and 8, in which parameter regions in the α 21 − m 0 plane the impact of reducing the error for the determination of Σ is large, and how large is the effect. Therefore, we show in the lower panels of Fig. 8 the iso-contours of ∆f CPX (Σ), the difference between the results shown in the upper panels of Figs. 7 and 8, defined in the previous paragraph. The panels (d), (e) and (f) in Fig. 8 are for r NME =2.0, 1.5 and 1.3, respectively. They illuminate more clearly how much the sensitivity improves with reduced σ Σ (larger ∆f CPX (Σ) > 0 implies larger improvement). We confirm that for α 21 < ∼ (0.5-0.6)π, the impact of reducing σ Σ is rather small, ∆f CPX (Σ) < ∼ 0.1, whereas for α 21 > ∼ (0.5-0.6)π, there are some regions where the improvement is relatively large, or ∆f CPX (Σ) is larger than 0.2 (inside the regions delimited by the solid orange curves in the lower panels, (d), (e) and (f) of Fig. 8).
Finally, let us make a comparison, at 3 σ CL, by changing σ Σ from 0.05 eV to 0.02 eV and/or r NME from 1.3 to the extreme value of 1.1, the possibility (2) mentioned above. In Fig. 9 we show the cases of (σ Σ , r NME ) = (0.05 eV, 1.3), (0.05 eV, 1.1), (0.02 eV, 1.3) and (0.02 eV, 1.1), respectively, in the panels (a), (b), (d) and (e). As done for Fig. 8, we only consider the case of the true value of α 31 = 0 since the overall behaviours of theses plots does not depend so much on the assumed value of α 31 = 0. Since we are considering the extreme values of σ Σ and r NME we show our results at 3 σ CL.
The panels (g) and (h) are to show the difference of f CPX obtained when σ Σ , the uncertainty in Σ, get decreased from 0.05 eV to 0.02 eV, keeping r NME fixed; Exhibited in the panels (g) and (h) are the iso-contours of ∆f CPX (Σ) for the cases with r NME = 1.3 and 1.1, respectively. Whereas the panels (c) (the case of σ Σ = 0.05 eV) and (f) (σ Σ = 0.02 eV) of Fig. 9 are to show the difference of f CPX when r NME = 1.3 is varied to 1.1 by keeping σ Σ fixed, ∆f CPX (r NME ) ≡ f CPX (r NME = 1.1) − f CPX (r NME = 1.3), an ultimate improvement of the NME calculation. Finally, in the panel (i), the iso-contours of ∆f CPX (Σ, r NME ) ≡ f CPX (σ Σ = 0.02 eV, r NME = 1.1)−f CPX (σ Σ = 0.05 eV, r NME = 1.3) are displayed to illuminate the impact of reducing both σ Σ and r NME simultaneously.
(i) we can see that the simultaneous reduction of both σ Σ and r NME increase considerably the regions corresponding to ∆f CPX ≥ 0.3.

VIII. IS IT POSSIBLE TO CONSTRAIN THE NUCLEAR MATRIX ELEMENT?
Finally, let us consider the possibility to constrain the nuclear matrix element only from the experiments. So far we have been assuming that the NME is restricted to some range since, in principle, it can be computed theoretically. The question we ask in this section is: "Is it possible to obtain some non-trivial constraint on the NME just by the experimental consistency even though we let the NME vary completely freely in the fit?" This is of course under the assumption of neutrino mass mechanism of 0νββ decay.
In order to answer this question, we have repeated the same analyses but without imposing any condition on the NME uncertainty parameter ξ ≡ M is the reference NME value defined in Eq. (25) in section IV B, used to translate the measured half life to m 0νββ and M (0ν) is the true NME value which is now assumed to be completely unknown.
In Fig. 10 we show the similar results shown in Figs. 2 and 3 but for the NME uncertainty parameter ξ, namely, ∆χ 2 as a function of ξ for different true values of α 21 = 0, π/4, π/2 and π for the case of the inverted mass hierarchy with m 0 = 0.1 eV (left panel) and 0.2 eV (right panel). The ranges of ξ corresponding to the cases of r NME = 2.0, 1.5 and 1.3 are also shown, respectively by the black dashed, dotted and solid lines in the same plots. We only show the case of the inverted mass hierarchy since that results for the normal hierarchy case is similar. We note that the case of the inverted mass hierarchy can constrain somewhat better the NME uncertainty parameter ξ than the case of the normal hierarchy (not shown). From this plot, we can see that in some regions, α 21 < ∼ 0.5π or around π, it is possible to restrict (although not so severely) the possible range of the NME only from the consistency of the results of the experiments. In Fig. 11 we show the ranges of the fitted value of ξ allowed at 1, 2 and 3 σ CL, for 1 DOF, indicated, respectively, by filled colors of yellow, red and light blue, as a function of the true values of α 21 for the case of α 31 (true) = 0, m 0 (true) = 0.1 eV (left panel) and 0.2 eV (right panel) for the case of the inverted mass hierarchy. From this plot, if the true value of α 21 ∼ 0, we can constrain better the upper value of ξ whereas if α 21 ∼ π we can constrain better the lower limit of ξ though the restricting power is not so strong.
To summarize: Assuming the standard three flavor scheme and standard mechanism of light neutrino exchange for 0νββ decay, it is possible to perform a consistency check of the NME calculation and one may be able to exclude some version of the NME obtained. If good fit is obtained only for some value of ξ well beyond the range expected from the theoretical uncertainty, this could imply the presence of new physics such as, for example, the sterile neutrino contribution. See e.g, Refs. [82][83][84] where the impact of the sterile neutrinos on double beta decay was studied.

IX. CONCLUSIONS AND DISCUSSION
In this paper, we have studied how and to what extent the Majorana phases can be constrained under the assumption of highly accurate measurement of 0νββ decay rate and stringent constraint on sum of neutrino masses by cosmological observation. The remarkable feature of our results is that the synergy between cosmology and double beta decay experiments in constraining the Majorana phase is quite visible; We have demonstrated that it is possible to obtain highly non-trivial constraints on the Majorana phase, if accurate measurement of Σ is added, even if double beta decay experiment alone can not say anything about the Majorana phase. Notice that double beta decay experiment by itself cannot place any useful constraints on the Majorana phase in our current setting, not because they suffer from the NME (Nuclear Matrix Element) uncertainty (which of course makes it worse if it exists), but because they lack independent key information of the absolute neutrino mass scale in itself.
The general tendency of the sensitivity to the Majorana phases are as follows: • Dependence of m 0νββ on α 31 , one of the two Majorana phases, is very weak due to the suppression by s 2 13 . Hence, little sensitivity can be expected for α 31 , and the sensitivity to α 21 depend on α 31 but only very weakly.
• The regions of the best and the worst sensitivities to α 21 are located at around the true values α 21 0 or π, and at α 21 2π/3, respectively, in agreement with our analytic estimate in Sec. III B.
• For both mass hierarchies and within any uncertainties of the NME used, the better sensitivity to the CP phase α 21 is obtained for larger value of m 0 , apart from minor exception of region m 0 < ∼ 0.05 eV, mostly in the horizontal branch in Fig. 1 for the inverted mass hierarchy.
We took as the reference set up that the effective masses observed in 0νββ and β decay experiments can be measured with the uncertainty of 0.01 eV and 0.06 eV, respectively, whereas the sum of neutrino masses Σ can be determined with the uncertainty of 0.05 eV by cosmological observation. We have assumed the uncertainty r NME , defined as the ratio of the maximum and minimum values of the theoretically expected NME values, between factor of 1.3 − 2.0 and treated it in a different way from these experimental errors (which is similar to the treatment of flux normalization uncertainty). To show the dependence of these assumed errors on sensitivity to Majorana phase, we have examined in the last part of our analysis, the case with more aggressive errors, 0.02 eV for the error of Σ and r NME = 1.1 for the NME uncertainty. We note here that we took a fixed error (0.01 eV) for m 0νββ based on the discussion given in Appendix A. If a percentage error for m 0νββ is to be assumed, the sensitivity to α 21 in region m 0 < ∼ 0.05 eV would be greatly enhanced.
To display the sensitivities to the Majorana phase quantitatively we have used the two measures, (1) the conventional allowed region plots, and (2) CP exclusion fraction f CPX , a fraction of the CP phase space that can be excluded for a given set of input parameters. We believe that f CPX is a useful global measure, as the CP fraction widely used in the sensitivity studies for the long baseline neutrino oscillation experiments is, but in a different way. That is, f CPX is better suited to reveal (relatively poor) sensitivities in the early era of the search for the CP phase effect, as discussed in Ref. [48]. Since measuring the Majorana phase is very challenging in any means, we believe that even the partial exclusion of its possible range is quite useful.
We have shown with our reference setup and r NME = 1.5 that α 21 can be constrained by excluding 10 − 40% of the phase space of α 21 at 2σ CL for the lowest neutrino mass of 0.1 eV for both mass hierarchies. Even if we use 0.02 eV for the error of Σ the excluding fraction does not change much, and is 10 − 50%. Fortunately, the f CPX contours are rather stiff against change in CL from 2σ to 3σ so that they stay more or less the same, as shown in Fig. 8, indicating robust nature of sensitivity to the Majorana phase resistant to high CL demand. We also should note that the sensitivity to α 21 becomes significantly better when the uncertainty of the nuclear matrix element (r NME ) is reduced from the factor 2 to 1.5, independent of the assumed mass hierarchy and value of α 31 .
It is only recently that such a discussion started to obtain real perspective because cosmology entered into the precision era, as demonstrated most notably an epoch making precision measurement achieved by Planck. Yet, the accuracy of measurement of the sum of neutrino masses which warrants reasonably strong restriction to the Majorana phase is rather demanding, Σ 0.01 − 0.02 eV, which probably requires the Euclid satellite as well as the next generation galaxy surveys, in addition to the better understanding of the different systematics exist in the different set of the cosmological data. Of course, it also requires accurate measurement of lifetime of 0νββ decay in a ton scale experiment with very low background, which would allow uncertainty in measurement of m 0νββ as small as 0.01 eV.
One may argue that too much relying on cosmological observation in determining Σ is dangerous because, neutrinos being a minor player in the universe, its precise determination is possible only in a model dependent way. Though it is a valid point we must bear in mind that future precision observation itself offers an even more stringent test of the ΛCDM paradigm. Therefore, we have implicitly assumed as a prerequisite of the analysis that the paradigm (or even an alternative cosmological model) will be established in a sufficient robustness to allow reliable measurement of the neutrino component. The discussion in this paper would serve as a prototype of the similar analysis after we have acquired the established standard cosmological model.
Finally, we also studied the possibility to constrain the nuclear matrix element only from the consistency of the experimental results by varying its value completely freely. We found that in some limited situation where α 21 0 or π, it is possible to obtain non-trivial constraint on the nuclear matrix element, which can be used as a consistency check of its theoretical calculation. 100 kg · yr ε 2 det · m Ge · t exp 1 4 eV, whereas for 136 Xe (by using G 0ν = 14.58 × 10 −15 yr −1 [54]) we obtain, m min 0νββ ∼ 0.24 3.0 M (0ν) b 0.01 keV · kg · yr 1 4 ∆E 100 keV 1 4 100 kg · yr ε 2 det · m Xe · t exp 1 4 eV, where m Ge (m Xe ) is the total mass of the isotope of 76 Ge ( 136 Xe) to be used in the double beta decay experiments. If 0νββ decay will be actually observed, at first approximation, we assume that these values could be roughly corresponds to the uncertainty on the measurement of m 0νββ , or σ 0νββ .
As we can see, as the sensitivity improves only as the one of the fourth power of the size of the experiment, background and energy resolution, it is seems not so easy to improve the sensitivity and looks quite difficult to reach the level of σ 0νββ ∼ O(0.01) eV, even if we consider ∼ 1 ton scale experiment. Therefore, it would be necessary to realize the background free or very low background experiment to achieve such a level (see below).
Next let us consider the case that the background rate b is so low that N BG in Eq. (A3) is negligible, or namely, zero background experiment [35]. In this case, from the number of observed events, N 0νββ , the half life time can be estimated as and its uncertainty is estimated as We note that the result does not depend on the life time, and the uncertainty can be arbitrarily small as we increase the product of the source mass and exposure time of the experiment by considering only the statistical error, until the point that the background could no longer be neglected. By using the same numbers for G 0ν and M (0ν) we used to obtain m min 0νββ in Eqs. (A6) and (A7) and ), for 76 Ge, we obtain, σ 0νββ ∼ 0.06 100 kg · yr ε det · m Ge · t exp 1 2 whereas for 136 Xe, we obtain, σ 0νββ ∼ 0.04 100 kg · yr ε det · m Xe · t exp We expect that these values should coincide roughly with the minimum possible observable values (or sensitivities) of the experiments, or σ 0νββ ∼ m min 0νββ . Therefore, as long as the background is neglected, by considering the ∼ 1 ton size experiment, it seems possible to reach the level of σ 0νββ ∼ O(0.01) eV. In this work, for definiteness and simplicity, we assume that m 0νββ can be determined with an accuracy of 0.01 eV for a given reference value of the NME. See the section IV B how to take into account the NME uncertainty. We note, however, that in reality, the fate of the background would not be so simple to allow analytic treatment as ours. Therefore, most probably, one needs to perform detailed numerical simulations to reliably estimate the experimental uncertainties in measurement of the double beta decay lifetime.