Improved Statistical Analysis for the Neutrinoless Double-Beta Decay Matrix Element of ¹³⁶Xe

Abstract

Neutrinoless double beta decay nuclear matrix element ($M_{0\nu }$) for ¹³⁶Xe was recently analyzed using a statistical approach (Phys. Rev. C 107, 045501 (2023)). In the analysis, three initial shell model effective Hamiltonians were randomly altered, and their results for 23 measured observables were used to infer credibility for the $M_{0\nu }$ nuclear matrix element (NME) based on a Bayesian Model Averaging approach. In that analysis, a reasonable Gamow-Teller quenching factor of 0.7 was assumed for each starting effective Hamiltonian. Given that the result of the statistical analysis was sensible to this choice, we are here improving that analysis by assuming that the Gamow-Teller quenching factor is also randomly chosen within reasonabe limits for all three starting Hamiltonians. The outcomes are slightly higher expectation values and uncertainties for the $M_{0\nu }$ NME.

1. Introduction

Neutrinoless double beta decay ($0\nu \beta \beta$) is undoubtedly one of the most intense area of research in nuclear and particle physics [1]. The eventual observation of this process would indicate that the lepton number is not conserved, unlike assumed in the Standard Model of particle physics, but it is changing by two units in this decay. This conclusion by itself would indicate physics beyond the Standard Model (BSM), which could lead to many ramifications, such as establishing the Majorana nature of the neutrino via the black-box theorems [2,3], addition of new Majorana mass terms in the BSM Lagrangian [4,5], and a possible explanation of the observed baryonic asymmetry through leptogenesis [6].

An associated process, is the two neutrino double beta decay ($2\nu \beta \beta$), which is allowed by the Standard Model, and it was observed for 11 isotopes [7]. In addition, three double electron capture decays with emissions of two neutrinos were also observed [7]. For both types of decays the inverse half-life (decay constant) can be described as a product of a phase space factors (PSF), a nuclear matrix element (NME) squared, and a BSM lepton number violation (LNV) parameter in the case of the $0\nu \beta \beta$. A good theoretical description of the $0\nu \beta \beta$ half-lives is essential for providing guidance for effective choices of the expensive experimental setups.

While the PSF can be calculated with increased accuracy [8,9,10,11,12,13] there are still significant discrepancies between the NME calculated with different nuclear structure models. Among the nuclear structure methods used for calculating the NME are the interacting shell model methods [14,15,16,17,18,19,20,21,22,23,24], the proton-neutron Quasiparticle Random Phase Approximation (pn-QRPA) methods [10,25,26,27,28,29,30], the Interacting Boson Approximation (IBA) methods [31,32], the Energy Density Functional method [33], the Projected Hartree Fock Bogoliobov (PHFB) method [34], the Coupled-Cluster method (CC) [35], the in-medium generator coordinate method (IM-GCM) [36], and the valence-space in- medium similarity renormalization group method (VS-IMSRG) [37]. As a general rule, each of these methods use a preferred methodology and a preferred effective nuclear Hamiltonian constructed in a preferred model space to describe the same NME.

Recently, there is a larger and consistent interest in the community to study these observables from a statistical point of view, thus considering fluctuation in the underlying nuclear dynamics and providing not only values for these observables, but also uncertainties [38]. In two recent papers [39,40] we applied this strategy to the $0\nu \beta \beta$ decay cases of ⁴⁸Ca and ¹³⁶Xe. In particular, in Ref. [40] we analyzed the case $M_{0\nu }$ of ¹³⁶Xe using a a Bayesian Model Averaging approach [38,41]. The outcome was a range of values at 90% confidence level for the $M_{0\nu }$ NME, an expectation value and an uncertainty.

In the analysis of Ref. [40], three starting shell model effective Hamiltonians were used, whose two-body matrix elements were randomly altered by $\pm 10$%, and their results for 23 measured observables were used to infer credibility for the $M_{0\nu }$ based on a Bayesian Model Averaging approach. In that analysis, a reasonable Gamow-Teller quenching factor of 0.7 was assumed for each starting effective Hamiltonian. Given that the results of the statistical analysis were sensible to the choice of this quenching factor, we here present an improved analysis assuming that the Gamow-Teller quenching factor is also random variable chosen within reasonable limits between 0.35 and 0.75 for all three starting effective Hamiltonians. The outcome will be a new range of values at 90% confidence level for the $M_{0\nu }$ NME, a new expectation value and uncertainty. As in Ref. [40] the calculations are performed using three independent starting effective Hamiltonians, i.e., the SVD [42], the jj55t [20] and the GCN5082 [43], all valid in the same single particle (s.p.) model space, jj55, that is appropriate for ¹³⁶Xe. All three effective Hamiltonians are based on theoretical Brueckner G-Matrix effective Hamiltonians that were fine-tuned to describe the experimental excitation energies of a large number of nuclei that can be described in this s.p. model space.

The paper is organized as follows. In Section 2 we present a quick review of the methodology of calculating the 23 observables used in this study, and the statistical model. In Section 3 we present new results and analyze their significance, followed by a statistical analysis based on the Bayesian Model Averaging approach. Section 4 is devoted to conclusions and outlook.

2. The Statistical Model

Following Ref. [40], we extend our statistical study the neutrinoless NME, $M_{0\nu }$, for the ¹³⁶Xe isotope that is vigorously investigated experimentally by several current and future double beta decay (DBD) experiments [44,45,46,47,48]. As in Ref. [40] the calculations are done with the interacting shell model in the $jj55$ model space consisting of the $0g_{7/2},1d_{5/2},1d_{3/2},2s_{1/2}$ and $0h_{11/2}$ orbitals, which assumes ¹⁰⁰Sn as a core, and can be used to describe nuclei with number of neutrons and protons between 50 and 82. Given that the $0g_{7/2}$ and $0h_{9/2}$ spin-orbit partner orbitals are missing, one could be concerned that in the $jj55$ model space some missing Gamow-Teller strength from the Ikeda sum-rule occurs, but prior studies show that this deficit can be compensated by using effective Gamow-Teller operators (in our case a less smaller quenching factor).

The starting effective Hamiltonians used in this study are the same as in Ref. [40], namely SVD [42] (short name $svd$), jj55t [20] (short name $j5t$), and GCN5082 [43] (short name $gcn$). A full description of these starting effective Hamiltonian can be found in section II of the above Ref. [40]. The “starting” qualifier for the effective Hamiltonians indicates that they were fine-tuned for a range of nuclei appropriate for the $jj55$ s.p. model space. In our statistical study, these three starting effective Hamiltonians are modified by adding $\pm 10\%$ random perturbations to their two-body matrix elements (TBME). The rationale for limiting the changes of the TBME of the starting effective Hamiltonians to $\pm 10\%$ is given in Ref. [39]. Here, as in Ref. [40] we generate 1000 effective Hamiltonians by adding random perturbations to the TBME of the starting effective Hamiltonians. As in our previous statistical studies for ⁴⁸Ca [39] and ¹³⁶Xe [40], the single-particle energies were kept unmodified, to avoid changes in the magicity of the ¹⁰⁰Sn core.

The 24 observables used in this study are described in section II of Ref. [40], but we also list them here for an easier reading of Figure 1, Figure 2 and Figure 3, below: $0\nu \beta \beta$ NME, $2\nu \beta \beta$ NME, the energies of the first $2^{+}$, $4^{+}$, and $6^{+}$ states in the parent (¹³⁶Xe) and daughter (¹³⁶Ba) nuclei, B(E2)↑ transition probabilities for ¹³⁶Xe and ¹³⁶Ba to the first $2^{+}$ states, the Gamow-Teller transition probability for the transition from ¹³⁶Xe and from ¹³⁶Ba to the $1^{+}$ excited state in ¹³⁶Cs, and the neutron and proton occupancies for ¹³⁶Xe and ¹³⁶Ba above the ¹⁰⁰Sn core in the jj55 model space shells. Details on how these observables are calculated can be also found in section II of Ref. [40].

One important ingredient in these calculation is the $2\nu \beta \beta$ NME entering into $2\nu \beta \beta$ decay half-life [17,39]:

$$T_{2\nu }^{-1}=G_{2\nu }{g}_A^4{\mid M_{2\nu }\mid }^2$$

(1)

Here, $G_{2\nu }$ is the phase space factor [11,13], and $g_A=1.276$ [49] is the free nucleon weak axial coupling constant. The $2\nu \beta \beta$ NME, $M_{2\nu }$, is calculated as

$$M_{2\nu }=\sum _k\frac{〈0_f^{+}\mid {\left(\sigma {\tau }^{-}\right)}^{eff}\mid 1_k^{+}〉〈1_k^{+}\mid {\left(\sigma {\tau }^{-}\right)}^{eff}\mid 0_i^{+}〉}{E_k-E_0},$$

(2)

where the summation is on the $1_k^{+}$ states in ¹³⁶Cs, and $E_0=Q_{\beta \beta }{/2+\Delta M(}^{136}$Sc${-}^{136}$Xe). Details on how the sum one the intermediate $1_k^{+}$ state is performed in the $2\nu \beta \beta$ Equation (2) are given in section IV of Ref. [50]. In these calculations, the Gamow-Teller ${\tau }^{-}\sigma$ operator in Equation (2) was quenched by the quenching factor q, which is generally found to be around 0.7 for most s.p. model spaces and most effective Hamiltonians. The interacting shell model point of view is that this quenching effect is due to the renormalization of the ${\tau }^{-}\sigma$ operator in reduced model spaces [51,52,53], while the axial coupling constant remains that for free nucleons. Therefore, in this study, the “renormalized” ${\tau }^{-}\sigma$ operator in Equation (2) is ${\left(\sigma {\tau }^{-}\right)}^{eff}=q{\tau }^{-}\sigma$. However, in the jj55 model one needs a quenching factor of about 0.4 for the GCN5082 Hamiltonian to describe the experimental value of the $2\nu \beta \beta$ NME (and an intermediate value between 0.4 and 0.7 for the jj55t Hamiltonian). In Ref. [40] we used the reasonable q-value of 0.7 for all three effective Hamiltonians. In this paper, we improve the analysis by allowing the quenching factor q to take random values between 0.35 and 0.75. The precise methodology on how this change is implemented is described in the next section.

For the $0\nu \beta \beta$ NME we cannot not use a similar effective operator described by a simple quenching factor of the Gamow-Teller operator, but that role is played by the short-range two-body correlation function [39,40]. In the case of ⁴⁸Ca presented in Ref. [39] we used a short range correlation function that produces results close to those of the more sophisticated ab-initio method. Given that there are not similar ab-initio results for ¹³⁶Xe, here, as in Ref. [40] we choose a short-range correlation function based on the widely utilized CD-Bonn parametrization (see e.g., [22,23,24]). The $0\nu \beta \beta$ NME is calculated using the closure approximation [21,54,55,56].

3. Results of the Statistical Analysis

3.1. Description of Experimental Data and Statistical Results

The experimental data used in this study is listed in the first two columns of Table I of Ref. [40]. The sources for this data and its description can be found in Section III of the same Ref. [40]. Table I of Ref. [40] also contains the expectation values and the standard deviations for all 24 observables calculated with the 1000 sample for each starting effective Hamiltonian.

More statistical results are shown here for the GCN5082 Hamiltonian in Figure 1, Figure 2 and Figure 3, while the similar results for the SVD starting effective Hamiltonian can be found in Figures 1–3 of Ref. [40]. In these figures we used short notations that identify each observable, whose descriptions can also be found in Ref. [40], but they are also briefly repeated here for an easier reading of these figures. The “parent nucleus” is ¹³⁶Xe, the “daughter nucleus” is ¹³⁶Ba, and the “intermediate nucleus” is ¹³⁶Cs for the $\beta \beta$ transitions. We use the letter “P” at the beginning of a label for an observable indicates that it is related to the parent nucleus, and the letter “D” indicates observables corresponding to the daughter nucleus. $M_{0\nu }$ represents the $0\nu \beta \beta$ NMEs, and $M_{2\nu }$ represents the $2\nu \beta \beta$ NMEs. P$GT$ and D$GT$ label the Gamow-Teller strengths to the first excited $1^{+}$ state in the ¹³⁶Cs intermediate nucleus from ¹³⁶Xe and from ¹³⁶Ba, respectively. P$B\left(E2\right)\uparrow$ and D$B\left(E2\right)\uparrow$ are the electric quadrupole transition probabilities ($0^{+}\to 2^{+}$) for ¹³⁶Xe and ¹³⁶Ba, respectively. P$E_{2^{+}}$, P$E_{4^{+}}$, P$E_{6^{+}}$ and D$E_{2^{+}}$, D$E_{4^{+}}$, and D$E_{6^{+}}$ denote the energies of the first $2^{+}$, $4^{+}$, and $6^{+}$ excited states, for ¹³⁶Xe and ¹³⁶Ba and respectively. POP$g_7$, POP$s_1$, POP$h_{11}$, and ${\mathrm{POP}}_d$ stand for the proton occupation probabilities of the $0g_{7/2},2s_{1/2}$, $0h_{11/2}$, and d orbitals in ¹³⁶Xe, while DOP$g_7$, DOP$s_1$, DOP$h_{11}$, and ${\mathrm{DOP}}_d$ are the proton occupation probabilities of the $0g_{7/2},2s_{1/2}$, $0h_{11/2}$, and d orbitals in ¹³⁶Ba. DVN$g_7$, DVN$s_1$, DVN$h_{11}$, and ${\mathrm{DVN}}_d$ represent the neutron vacancy probabilities in ¹³⁶Ba. The experimental proton occupancies for $1d_{5/2}$ and $1d_{3/2}$ orbitals cannot be separated experimentally, and they are listed here as $POPd$ and $DOPd$.

Figure 1, Figure 2 and Figure 3 present detailed statistical results obtained with the GCN5082. Similar figures for the SVD starting Hamiltonian can be found in Ref. [40], where a full description of these figures can be also found. Figure 1 shows the heat-map, i.e., the Pearson correlation coefficients, for all 24 observables in the case of the GCN5082 starting effective Hamiltonian. Figure 2 shows the correlation matrix only for the observables that have a Pearson coefficient larger than 0.5 relative to the $M_{0\nu }$ NME. Figure 3 presents in blue the Kernel Distribution Estimates (KDE) for all 24 observables (the diagonals in Figure 2). The red curves are Gaussians plotted using the experimental values as expectation values, and the experimental errors used as standard deviations. These experimental values are listed in Table I of Ref. [40], and are not repeated here.

The main difference between the GCN5082 statistical results presented here, and SVD results shown in Ref. [40] is that the Pearson correlation coefficient between the $M_{0\nu }$ NME and the $M_{2\nu }$ NME is not as strong for GCN5082 starting effective Hamiltonian as it was for SVD starting effective Hamiltonian, likely because of the fixed quenching factor of 0.7 used in these calculation. However, in Section 3.2 below we improve our analysis by letting the quenching factor q to take random values in the interval 0.35 to 0.75 in our Bayesian Model Averaging analysis. Other observables that significantly correlate with the $M_{0\nu }$ NME are: the $E_{2^{+}}$, $E_{6^{+}}$ and $E_{4^{+}}$ excitation energies of the daughter nucleus (¹³⁶Ba), the $E_{2^{+}}$, $E_{4^{+}}$ and $E_{6^{+}}$ excitation energies of the parent nucleus (¹³⁶Xe), and the daughter neutron vacancy of the $0g_{7/2}$ orbital.

3.2. The Bayesian Model Averaging

A thorough statistical analysis of the $0\nu \beta \beta$ NME could provide valuable insights into its range, mean value, and its uncertainty. It appears that the values of all observables listed in Figure 3 remain robust in response to small variations in the GCN5082 starting effective Hamiltonian. There is though no indication of any significant deviation from its starting values (e.g., sudden deviations indicating phase transitions, double bumps, etc.), similar to the results for the SVD effective Hamiltonian shown in Ref. [40]. However, the experimental data with errors (red curves) are not always on top of the theoretical distributions, especially the $2\nu \beta \beta$ NME that is proportional to the quenching factor squared ($q^2$). In Ref. [40] we used a Bayesian Model Averaging analysis that includes the statistical results of all three starting effective Hamiltonians to provide a range of values, a mean value and its uncertainty for the $0\nu \beta \beta$ NME. For that we introduced a common model distribution function that mixes the distributions of the $0\nu \beta \beta$ NME for each starting effective Hamiltonian depicted in Figure 3 using weighting factors $W_H$,

$$\begin{array}{cc}\hfill P\left(M_{0\nu }\right)=& W_{svd}{P}_{svd}\left(M_{0\nu }\right)+W_{gcn}{P}_{gcn}\left(M_{0\nu }\right)+W_{j5t}{P}_{j5t}\left(M_{0\nu }\right).\hfill \end{array}$$

(3)

The normalized weights $W_k$ with $k=svd,gcn,j5t$ can be obtained using the the Bayesian Model Averaging method [38,41] described in Section IV of Ref. [40]. The main difference between the present analysis and that of Ref. [40] is that here we allow the quenching factor q to vary uniformly in a reasonable range of values between 0.35 and 0.75. This range is chosen while for q = 0.4 the GCN5082 Hamiltonian better describes the $2\nu \beta \beta$ NME that is strongly correlated with the $0\nu \beta \beta$ NME, for q = 0.7 the SVD Hamiltonian better describes the $2\nu \beta \beta$ NME, and the jj55t Hamiltonian is favored by intermediate values of the quenching factor q.

As in Ref. [40] we associate the weights $W_k$ with the Bayesian posterior probabilities, $p\left({\mathcal{M}}_k\right|(y_e,{\sigma }_e)$, given some prior probabilities $\pi \left({\mathcal{M}}_k\right)$ (see Equations (3)–(5) of Ref. [40]). To estimate these posterior probabilities one needs to calculate the so-called evidence integrals (Equation (4) of Ref. [40]) using Monte Carlo techniques. In Ref. [40] the evidence integrals were obtained integrating on all $N_{TBME}$ number of two-body matrix elements describing each effective Hamiltonians, denoted by ${\theta }_j$. Here we also include the integration over q by randomly choosing its value for each Monte-Carlo sample:

$$p(y_e,{\sigma }_e|{\mathcal{M}}_k)=\int \prod _i^{N_{obs}}d{y}_i{p}_{y_e,{\sigma }_e}\left(y_i\right)\left[\int dq\prod _j^{N_{TBME}}d{\theta }_{j}p\left(y_i\right|{\theta }_j,q,{\mathcal{M}}_k)\pi ({\theta }_j,q|{\mathcal{M}}_k)\right],$$

(4)

Here $N_{obs}$ refers to those obsevables that have a Pearson coefficient larger than 0.5 relative to the $M_{0\nu }$ NME (see Figure 2). Assuming as in Ref. [40] that the priors are chosen democratically equal, i.e., 1/3 for each model described by a starting effective Hamiltonian, one gets for the posterior probabilities

$$p\left({\mathcal{M}}_k\right|y_e,{\sigma }_e)=\frac{1}{1+{\sum }_{j\ne k}\frac{p(y_e,{\sigma }_e|{\mathcal{M}}_j)}{p(y_e,{\sigma }_e|{\mathcal{M}}_k)}}.$$

(5)

where the indexes k and j go over the collection of model: $svd$, $gcn$, and $j5t$. Our new calculations of the evidence integrals provides the following ratios:

$$\frac{p(y_e,{\sigma }_e|{\mathcal{M}}_{svd})}{p(y_e,{\sigma }_e|{\mathcal{M}}_{gcn})}=6\times {10}^{-10}$$

(6)

and

$$\frac{p(y_e,{\sigma }_e|{\mathcal{M}}_{j5t})}{p(y_e,{\sigma }_e|{\mathcal{M}}_{gcn})}=2\times {10}^{-5}.$$

(7)

Based on these results one can conclude that the $gcn$ model has a much better posterior credence than the other two models, likely because it describes better most of the other observables, except for the $2\nu \beta \beta$ NME when one takes q much larger than 0.4. Therefore, the normalized posterior probabilities become $p\left({\mathcal{M}}_{gcn}\right|y_e,{\sigma }_e)\approx 1$, $p\left({\mathcal{M}}_{svd}\right|y_e,{\sigma }_e)\approx 0$, and $p\left({\mathcal{M}}_{j5t}\right|y_e,{\sigma }_e)\approx 0$. Given the ascribed connection between the posterior probabilities of the Bayesian Model Averaging method and the weights in Equation (3), one could surmise that all $W_k=0$, except the $W_{gcn}=1$. However, in the spirit of the predictor-corrector approach of a step-by-step evolution, as described in Ref. [40], we consider for the weights $W_k$ an average between the prior probabilities $\pi \left({\mathcal{M}}_k\right)$ and the posterior probabilities $p\left({\mathcal{M}}_k\right|y_e,{\sigma }_e)$, where $k=svd,gcn,jj5$. Therefore, we take $W_{gcn}=4/6,W_{svd}=W_{jj5}=1/6$. Certainly, if new relevant observables were to be considered, such as the ordinary muon capture [57], and they correlate significantly with the $M_{0\nu }$ NME, one could further evolve these posterior probabilities to finding improved ranges, mean values and uncertainties.

An updated Figure 4 of Ref. [40] is presented in Figure 4, which shows the probability distribution functions (PDF) for the three starting effective Hamiltonians. The red curve shows the weighted sum of Equation (3). The PDF were calculated using kernel-density estimates [58,59] for the histograms describing the $M_{0\nu }$ NME for each starting effective Hamiltonians, an example of which is in the upper-left panel of Figure 3. Using the results of our statistical analysis summarized in Figure 4 (the red curve) one can conclude that with 90% confidence the $0\nu \beta \beta$ NME lies in the range between 1.65 and 2.85, with a mean value of about 2.20 and a standard deviation of 0.48. These values are slightly higher than those obtained in Ref. [40].

4. Conclusions and Outlook

In conclusion, we presented an improved statistical analysis of the $0\nu \beta \beta$ NME for ¹³⁶Xe that builds up on the first such analysis presented in Ref. [40], where we observed that the results are sensitive to the choice of the quenching factor for the effective Gamow-Teller operator. The main difference between the present work and that reported in Ref. [40] is that we allow the Gamow-Teller quenching factor to become a statistical parameter, for which we choose a reasonable range of 0.35 to 0.75. In addition, we present detailed statistical results for the GCN5082 starting effective Hamiltonian, which is favored by the Bayesian Model Averaging when the quenching factor is included in the evidence integrals.

As in Ref. [40] in the analysis we used three known starting effective Hamiltonians that were widely tested for tin isotopes and other nuclei near ¹³²Sn, namely SVD, GCN5082 and jj55t, and we randomly change their two-body matrix elements uniformly within ±10% of their original values. Using sample sizes of 1000 Hamiltonians we analyzed for each starting effective Hamiltonian the correlations between $0\nu \beta \beta$ NME and other 23 observables that are accessible experimentally. We show that for the GCN5082 starting effective Hamiltonian, as in the case of the SVD effective Hamiltonian [40], the $0\nu \beta \beta$ NME still correlates significantly with the $2\nu \beta \beta$ NME, and the $2^{+}$, $4^{+}$, $6^{+}$ states in ¹³⁶Xe and ¹³⁶Ba. Finally, we compare the posterior probabilities for each starting effective Hamiltonian within the framework of the Bayesian Model Averaging method. Based on this statistical analysis we propose a common probability distribution function for the $0\nu \beta \beta$ NME, which favors the distribution of the GCN5082 starting effective Hamiltonian, it has a range of 1.65–2.85 at 90% confidence level, with a mean value of 2.20 and a standard deviation of 0.48. These results are somewhat larger from those of Ref. [40]: the predicted $0\nu \beta \beta$ NME just increased form 1.99 to 2.2, the uncertainty got a little larger from 0.37 to 0.48, so in sum the range is slightly larger and shifted towards larger values

In principle, the Bayesian Model Averaging results can be further improved if new reliable experimental observables become available and they show significant correlation with the $0\nu \beta \beta$ NME. In that case the present posterior probabilities can be used as prior probabilities, and the new data can be used to update them. One example of such data could be the ordinary muon capture rate [57]. However, one needs to first show that this observable is as robust as the other to small changes in the effective Hamiltonians. Another observable that is often though to be used in calibrating the $0\nu \beta \beta$ NME is the B(M1) [38]. Unfortunately, there are not many cases where B(M1) strengths available for the states involved in the $0\nu \beta \beta$ decays, i.e., the $0^{+}$ states in the parent and daughter and the $1^{+}$ states in the intermediate nucleus. However, one can still consider other M1 transitions if they correlate significantly with the $0\nu \beta \beta$ NME. One could think of investigating changes in the TBME random contributions from $\pm 10\%$ to say $\pm 15\%$. We plan to pursue these venues of research in the near future. Similar analyses for other $0\nu \beta \beta$ cases, such as ⁸²Se, will be reported in the near future.

Finally, one should mention that here and in Ref. [40] we only analyzed the $0\nu \beta \beta$ NME for the mass mechanism. It is well known that other mechanisms can contribute to the $0\nu \beta \beta$ decay [24,60]. One then could ask if contributions from different mechanisms could interfere and cancel each others. In Ref. [61] we showed that the maximum interference between the mass mechanism and so called $\eta$ mechanism could be close to 50% for ¹³⁶Xe. One should investigate if this potential cancellation could be even much higher if small changes in the effective Hamiltonians are allowed.