Unraveling the $0\nu\beta\beta$ Decay Mechanisms

We discuss the possibilities of distinguishing among different mechanisms of neutrinoless double beta decay arising in the effective field theory framework. Following the review and detailed investigation of the particular ways of discrimination, we conclude that the 32 different low-energy effective operators can be split into multiple groups that are in principle distinguishable from each other by measurements of the phase-space observables and by comparison of the decay rates obtained using different isotopes. This would require not only a substantial experimental precision but necessarily also a considerable improvement of the current theoretical knowledge of the underlying nuclear physics. Specifically, the limiting aspect in our approach turns out to be the currently unknown or uncertain values of low-energy constants. Besides the study adopting the effective field theory language we also look into several typical UV models.


Introduction
The unknown origin of neutrino masses, being one of the major puzzles of contemporary particle physics, strongly motivates the quest for lepton number violation in Nature. The prominent way of probing this symmetry is the search for neutrinoless double beta (0νββ) decay [1], observation of which would imply non-zero Majorana neutrino masses in accordance with the black-box theorem [2]. Besides the tight connection to neutrino masses as realized in the standard mass mechanism, neutrinoless double beta decay can be triggered in a variety of different ways, and thus potentially involve also other new physics. Generally, one can study higher-dimensional lepton-number-violating operators that can trigger 0νββ decay [3][4][5][6][7][8][9]. In fact, while the sole observation of 0νββ decay would indeed indicate that neutrinos acquire Majorana mass, it remains unclear whether the standard mechanism that gives a contribution proportional to the neutrino mass would be the dominant one. Examples of models beyond the standard model that can induce non-standard contributions to the 0νββ decay rate include, for instance, the left-right symmetric models [10][11][12][13] triggering several distinct mechanisms [5,14]. Sterile neutrinos can also contribute to 0νββ decay [14][15][16][17][18][19][20].
There is variety of experiments searching for 0νββ decay in different double-betadecaying isotopes [21][22][23][24][25][26][27][28][29][30]. Currently, the best limit on the half-life reaching 2.3×10 26 years, is claimed by the KamLAND-Zen collaboration [31] studying the decay of 136 Xe. The most stringent limit on the half-life of 0νββ decay of 76 Ge attains 1.8 × 10 26 years, as obtained by the GERDA collaboration [21]. Proposed next generation experiments such as LEGEND [32,33] ( 76 Ge), CUPID [34] ( 100 Mo), SNO+ [35] ( 130 Te) and nEXO [36] ( 136 Xe) aim towards testing half-lives of order of 10 27 − 10 28 years. Some experiments like NEMO-3 are also equipped with the technology to track individual electrons and measure the individual electron energy spectra and the opening angle between the two electrons, which can help to uncover new physics not only in 0νββ decay, but even in standard double beta decay [37]. Recent reviews of the experimental and theoretical efforts in the field of 0νββ decay can be found, e.g., in Refs. [38,39].
In this work we focus on different possibilities of experimental discrimination among different mechanisms inducing 0νββ decay. To do so, we adopt the effective field theory (EFT) framework developed in Refs. [6,7], which is briefly introduced in Section 2. In subsequent Section 3 we study the possible ways of distinguishing among the relevant set of low-energy EFT operators from 0νββ decay observables. After having discussed the single operator settings, we turn towards more complete models in Section 4. Finally, we summarize our findings in Section 5. This work has been carried out utilizing the upcoming NuBB code package [40].
2 EFT Approach to 0νββ Decay: The Master Formula

The half-life master formula
As we apply in this work the effective field theory approach introduced by [6,7], let us start by briefly summarizing the most important parts. Below the scale of electroweak symmetry breaking (EWSB) the 0νββ decay amplitude can be described in terms of an SU (3) C × U (1) Q invariant low-energy effective field theory (LEFT). Including operators up to LEFT dimension 9 the most relevant Lagrangians for 0νββ are given by [6,7] (2.1) and L (7) for the long-range part, where as well as the dimension 9 short-range Lagrangian i,R e R e c R + C (9) i,L e L e c L O i + C (2.5) Here α, β are color-indices and the t A are the generators of SU (3) in the fundamental representation given by the 8 Gell-Mann matrices λ A as t A = 1 2 λ A , A = 1...8. The operators O and O in (2.5) are related via parity transformation. Together with the standard mechanism of light Majorana neutrino-exchange, this framework contains 32 different LEFT operators that can trigger 0νββ decay. The transition from the quark level to the nuclear level can be achieved employing the chiral effective field theory (χEFT) [42]. The expected half-life contributed by the 32 effective operators is then captured by a "0νββ master-formula" combining the 32  coefficients, 6 different phase-space factors (PSFs) given in Table 1 and nuclear matrix elements (NMEs) summarized in Table 2. At the same time, low-energy constants (LECs) that describe the nuclear interactions within χEFT enter the formula -we summarize these in Table 3. The 0νββ half-life is then given in terms of different sub-amplitudes A i as T 0ν (2.6)

Sub-Amplitudes
The sub-amplitudes A i are categorized and defined via their corresponding leptonic currents. They each depend on the Wilson coefficients of different LEFT operators and can be  Table 2. NMEs used in our calculations based on the IBM2 model [43] written as SR , C T , C VL , C 1L , C 1L , C 2L , C 2L , C 3L , C 3L , C 4L , C 1R , C 1R , C 2R , C 2R , C 3R , C 3R , C 4R , C 5R , A me =M (6) me,L C VL + M (6) me,R C VR , 6 , C 7 , C 7 , C 8 , C 8 , C 9 , C (9) 9 . (2.7) The matrix elements M i depend on the different LECs and Wilson coefficients. We explicitly state the dependency on the different Wilson coefficients within the brackets in (2.7). A ν depends on the matrix elements  [47] g N N 1,6,7 O(1) g ππ 2 2.0 [47] g N N 2,3,4,5 [47] g N N V L O(1) g ππ 4 −1.9 [47] g N N where M ν represents the contribution from the standard mass mechanism. In contrast to the traditional approach employing the non-relativistic approximation, the EFT treatment contains also the contribution proportional to g N N ν , which parametrizes the contact-term contribution originating from the exchange of hard neutrinos [44,45]. A R is given by for A E the different contributions are

10)
A me is determined by (2.13) The short-range dimension-9 LEFT operators contribute to the C V,ππL,πN L,N N L couplings that appear in the chiral Lagrangian. They are given by (2.14) The two LECs g πN V andg πN V are defined as For the sake of convenience, we have marked all currently unknown LECs including g N N ν in bold within the above formulas.
In this work we will study the different LEFT operators at the matching scale of Λ = m W at which one would usually match the BSM model of interest onto LEFT. The running of the operators down to the scale of χPT at Λ χ 2 GeV is described in [7].

Relation to literature
A different basis to describe 0νββ decay developed first in [3] and [4] that is often used in the literature is defined by a set of 29 dimension-6 and dimension-9 lepton number violating LEFT operators given by for the long-range part with i, k ∈ {V ± A, S ± P, T L , T R , } and and the lepton currents j are given by

(2.19)
This framework does not include the dimension 7 operators of the framework utilized in our approach. While the remaining long-range part of the two descriptions can be related easily, the short-range operators are related to each other via Fierz transformations. One finds that (2.24) 1 We keep the two different types of indices for the short-range currents to stick with the literature (2.28) The main difference between the two set of operators is that the -basis contains short-range tensor operators instead of color octets.

Distinguishing the Effective Operators
Neutrinoless double-β-decay, if observed, would be characterized by several experimental observables, precise determination of which can give us some insight into the underlying BSM physics. Generally, the 0νββ decay experiments can be able to determine decay rate, single electron energy spectrum and angular correlation between the two emitted electrons. Additional information might be obtained by studying different ββ modes or by employing complementary information from other experiments such as the LHC. In the following, we will discuss possible ways of experimentally distinguishing among the 32 different 0νββ decay inducing LEFT operators with a focus on the limiting factors of a potential confirmation/exclusion of the existence of any additional non-standard scenario contributing to 0νββ decay alongside the standard mass mechanism.

Phase-Space Observables
While most experimental collaborations only attempt to measure the half-life of 0νββ decay, some experiments like NEMO-3 [51], or its future successor SuperNEMO [52], are designed to also measure the single electron energy spectrum and the angular correlation of the two outgoing electrons. These are associated with different electron currents and within the simplest approximation they can be calculated analytically. More exact solutions require numeric calculations of the exact electron wave functions [53]. The different PSFs G 0k can be written in the form [54] where p 1,2 and 1,2 are the momentum and energy of the first and second released electron, R is the radius of the final-state nucleus and E i,f denotes the energy of the initial-or finalstate nucleus, respectively. Here, we denote the part of the differential phase-space factor independent of the angle between the two outgoing electrons as g 0k , while h 0k is the angular correlation part proportional to the cosine of the opening angle θ. Additionally, G 04,06,09 have to be rescaled to comply with the definitions in [6,7] as The relations between the electron wave functions and the functions h 0k and g 0k are given in [54] to which we will refer here. We apply their simplest approximation scheme 'A' assuming a uniform charge distribution in the nucleus. Using Eq. (3.1) one can write the angular correlation coefficient a 1 /a 0 which is defined via as .

(3.5)
Here, ∆M Nuclei is the mass difference between the mother and daughter nuclei and Q ββ denotes the Q-value of the decay. The potential of utilizing the angular correlation of the outgoing electrons for discrimination between different mechanisms of 0νββ has been discussed e.g. in [55]. Similarly, the single electron spectra are given by Consequently, approximating the electron wave functions, we can easily calculate the expected angular correlation factor and single electron spectra for each of the 32 LEFT operators. The normalized single electron spectra as well as the angular correlations corresponding to each of the 6 distinct PSFs are shown in Figure 1. As we can see, using these observables the operators associated with distinct PSFs are in principle distinguishable from each other, provided substantial experimental accuracy is reached. However, distinguishing among different 0νββ mechanisms purely based on the phasespace observables has its obvious limitations. In fact, while G 06 is only induced in presence of multiple operators the dimension-6 vector operators both trigger several of the remaining PSFs. Taking this into account, we can identify 4 different groups of operators that are in principle distinguishable using the leptonic PSF observables, namely: C V R , the operators corresponding to G 01 and the ones corresponding to G 09 . The PSF observables that result from each of these 4 groups are shown in Figure 2. Here, we can see that the left-handed vector current operator C (6) V L and the operators corresponding to G 09 , while corresponding to distinct PSFs, are practically indistinguishable since the C (6) V L phase-space turns out to be dominated by the contribution from G 09 . The remaining groups are distinguishable from each other using at least one of the considered observables.
Note that while the electron wave functions depend on the charge of the daughter nucleus as well as on the decay energy, the general shape of the induced observables is not very dependent on the choice of the decaying isotope. In Figure 3 we show the single electron spectra and in Figure 4 the angular correlation coefficients corresponding to the 6 different PSFs in 4 different naturally occurring 0νβ − β − isotopes.

Decay Rate Ratios
The remaining 0νββ observable is the decay rate Γ itself. While the phase-space can be used to distinguish operators with different leptonic currents, information about the decay rates in various isotopes can be also applied to operators with distinct hadronic structures, as these give rise to different NMEs. The isotope dependence of the existing calculations of NMEs can be inferred from Table 2. Therefore, one can study the half-life ratios where The sums j,k are taken over all different PSFs generated by the operator O i and become relevant only for C  the unknown particle physics couplings, as was first discussed in [56] and shortly after also in [57]. Here, we take 76 Ge for the reference isotope. To be able to quantify how well one can distinguish two different operators O i,j from each other we can take the ratio Specifically, the ratios R im ββ relating the non-standard mechanisms with the standard mass mechanism will be of interest to compare the effect of different higher-dimensional operators and possibly identify the existence of additional exotic contributions to the 0νββ rate in experiments. Obviously, two operators O i,j would be indistinguishable via this method if the resulting ratio would equal unity, i.e., if R ij = 1. Vice versa, they would be perfectly distinguishable for either R ij → ∞ or R ij = 0, that is, for | log 10 (R ij )| → ∞.
Studying the decay rate ratios has several benefits. First of all, in case only one Wilson coefficient contributes at a time, it drops out. Therefore, the ratio corresponding to a certain operator and its Wilson coefficient is a constant that depends only on the corresponding NMEs, LECs and PSFs. If more Wilson coefficients contribute at the same time, then only the overall magnitude can be factored out. In this case, the relations between different coefficients can, of course, affect the resulting ratios. However, one can still utilize this method to study specific models and see if they are distinguishable from the standard mass mechanism. We will discuss this possibility in section 4. Additionally, when taking ratios of the half-lives, one can expect that the impact of correlated systematic relative errors on the NMEs decreases as they should (at least partially) cancel. In [58] it was shown that for the NME calculations using QRPA uncertainties arising from unknown g A quenching and nucleon-nucleon potentials are correlated among different isotopes. Half-life measurements in different isotopes as a tool to discriminate among different mechanisms of 0νββ decay have also been employed previously in [59][60][61][62]. Table 4. Operator groups that can possibly be distinguished via taking decay rate ratios. The choice of the groups depends on the knowledge of the LECs. If we set the unknown LECs to zero, the short-range scalar operator groups C S2−S5 become indistinguishable as well as the short-range vector operator groupsC (9) V andC (9) V . Improved knowledge of the LECs, assuming no fine tuning, would allow to distinguish among these operator groups.
Applying this approach to the master-formula framework one can identify 12 different groups of operators that can in principle be distinguished from each other. These groups are summarized in Table 4. However, the distinguishability of the short-range operators strongly depends on the currently unknown LECs. Taking the most of the unknown LECs to be zero while keeping g N N 6,7 = g πN V =g πN V = 1 (so that the contribution from the shortrange vector operators is not omitted) makes it impossible to distinguish the short-range scalar operators C (9) S2−S5 as well as the short-range vector operator groups C

Sensitivity on the unknown LECs
In Figure 5 we present the expected ratios R O i as well as the normalized ratios R im ββ defined in (3.7) and (3.8) for the above choice of LECs which will be our benchmark scenario. The ratios for the -basis are shown in Figure 6. Additionally, to study the uncertainties arising from the unknown LECs, the plots include 1000 points per operator group that each represent variations of the unknown LECs g i within the ranges − √ 10, −1/ √ 10 × |g i | and 1/ √ 10, √ 10 × |g i | , i.e. we vary the LECs within the range of values given by their expected order of magnitude shown in Table 3. For g N N ν which generates a short-range component into the standard mass-mechanism we take a variation of ±50%. The central values of the variation i.e. the median values are marked by crosses.
From the upper panel of Figure 5 one can infer that the half-life ratios R m ββ corresponding to the standard mass mechanism are not very sensitive to g N N ν (they are actually too small to be visible). In Figure 7 we explicitly show the impact of varying the g N N ν LEC on the expected half-life in 76 Ge for the standard mechanism. Again, compared to the impact of the unknown Majorana phases the effect of g N N ν is minor. However, it is important to note that the impact of g N N ν on the overall magnitude of the half-life cannot be ignored as easily. For comparison, we also present the case where g N N ν = 0 in Figure 7. For the remaining non-standard operators, however, we can see from Figure 5 that the values of the currently unknown LECs can have quite a significant impact on the expected ratios. Oftentimes, especially for the short-range C (9) i groups, the central values are significantly offset from our benchmark scenario with most unknown LECs turned off. Hence, for these operators the appearance of the unknown LECs has a significant impact on the corresponding 0νββ-decay rate. Although for some operator groups, such as C (9) 2S−5S , the spread of the values of the ratios obtained by varying the unknown LECs is relatively small, for other groups like the short-range vector contributions C (9) V ,C (9) V the variation of the unknown LECs results in a significant stretch around the central values. For these ratios the precise numerical value of the unknown LECs is of particular importance. The different sensitivities of the short-range scalar and vector operators arise from the fact that for the scalar operators some of the relevant LECs, namely those encoding pion-pion interactions g ππ i , are known, while for the short-range vector operators all relevant LECs are unknown. Since we do not fix the sign of the unknown LECs (except g N N ν ) there can be a gap within the LEC-varied ratios resulting in two visible central values for the operator groups, for which the ratios are sensitive to the sign of the LECs. The lower part of Figure 5 which displays the normalized R im ββ shows that the central values of the LEC-varied ratios are closer to 0 than the benchmark scenario. Therefore, the inclusion of the unknown LECs tends to impair the distinguishability from the standard mechanism.
The above discussion clearly shows the importance of determining the yet unknown LECs involved in the calculation. This can be achieved, for example, by lattice QCD calculations [63][64][65].

Distinguishing different operators
When studying the 0νββ decay rate the basic question to ask is whether and how well could a non-standard contribution be distinguished from the standard light-neutrino-exchange. Employing the half-life measurements in different isotopes, one can try to identify those that are most suitable for discrimination between the mass mechanism and an exotic 0νββ decay contribution triggered by a particular higher-dimensional operator. In the first row of Figure 8 we show the maximal ratios R max im ββ and the corresponding pair of isotopes obtained for each operator group. Here, we consider a "representative" scenario by studying the central values defined as the median ratio R im ββ of the range of values obtained from the variation of the LECs. At the same time, we identify the "worst-case scenario" ratio defined as the value within the range that is closest to unity, see the first column of Figure 8. In this context we consider only isotopes with existing experimental limits on the half-life, namely, the following: 76 [30]. Figure 8 also presents all the other ratios R max ij quantifying the mutual distinguishability of all the operator groups with the values corresponding to the representative scenario above the diagonal and the worst-case scenario below the diagonal. In addition, the dashed lines in Figure 8 mark the pairs of operators that could be discriminated using the phase-space observables. give the most distinct half-life ratios compared to the standard mass mechanism while the remaining long-range operator groups C (6) V L and C (6) T also result in sizable R max im ββ > 2. Additionally, both C (6) V L and C (6) V R could be identified by measuring the angular correlation of the emitted electrons. The non-standard short-range operators generally tend to have lower values of the ratios R max im ββ < 2; thus, there is less potential for identifying their contribution by experiments. However, the short-range vector operators C (9) V andC (9) V are associated with a different angular correlation than m ββ . On the contrary, the contributions from scalar short-range operators C (9) S2−S5 would be hardest to discriminate, as they do not manifest any significant isotope dependence on R im ββ and do not differ in the phase-space observables, either.
In the worst-case scenario, the operators in the group C (6,7) S,V , i.e., lepton number violating long-range scalar and vector interactions, are the only operators that result in R max im ββ > 2. In fact, this is the only operator group that is not affected by any unknown LECs. Besides C (6) T all remaining operators in this setting have expected ratios R max im ββ ≤ 1.3, which would require very precise measurements and accurate knowledge of the theoretically calculated half-lives to be able to claim a detection of any of these non-standard contributions. The V andC (9) V could be identified only based on measurements of the angular correlation and the scalar short-range operators, C (9) S1 , would be completely indistinguishable from the standard mechanism. In Appendix C we show for completeness the same results employing the full set of isotopes, for which there exist numerical values of NMEs computed using IBM2.
To be able to pinpoint the specific non-standard operator group O j contributing to 0νββ decay one needs to consider half-life ratios R ij for all different isotopes. Considering the central values, the best candidate to be clearly identified turns out to be the righthanded vector current C

The impact of nuclear uncertainties
The uncertainty induced by the nuclear part of the decay rate calculation i.e. the NMEs and LECs highly impacts and limits the above approach of distinguishing among different 0νββ mechanisms. The approach of comparing theoretically predicted ratios with experimentally measured ratios raises the question how well these theoretical uncertainties must be under control.
To study the impact of nuclear uncertainties, we can use the general formula for the half-life parameterized in terms of a Wilson coefficient C, the phase-space factor G and an  [67] and subsequent variants like the triaxial projected shell model (tpSM) [68] or realistic shell model (rSM) [69], the proton-neutron quasiparticle random phase approximation (pnQRPA) [70], the deformed QRPA (dQRPA) [71,72], the relativistic energy density functional method (rEDF) or covariant density functional theory (CDFT) [73,74], the non-relativistic energy density functional method (nrEDF) [75], the interacting boson model (IBM2) [43] and recently introduced ab initio approaches calculating NMEs from basic principles of χPT [76,77]. The grey bands mark the range of values covered by the different methods.
effective NME which we label M eff , Here, M eff is, generally, a weighted sum of combinations of different LECs and NMEs (see App. B for the explicit half-life equations of each single operator). If we consider the theoretical uncertainty of the half-life to be dominated by the uncertainty of M eff , we can determine the necessary theoretical accuracy of the nuclear physics. To estimate this, we assume M eff to be independent of the choice of the isotope, i.e., Then the necessary theoretical accuracy can be determined from the simple condition that the expected ratios should be distinguishable from unity within the theoretical uncertainty, Hence, the necessary theoretical accuracy for M eff reads Again, taking the central values as a baseline, the theoretical uncertainty on the overall nuclear part (that is from both LECs and NMEs) would need to be brought down to ∆M eff M eff ∼ 7% to be able to identify all possible non-standard contributions assuming single operator dominance. The exotic contribution easiest to identify using the half-life ratios would be the right-handed vector current C V R requiring an accuracy of ∆M eff M eff ∼ 22%. We want to emphasize, again, that this way of estimating nuclear uncertainties assumes that our calculations of the central ratios are a reasonable reflection of reality.
In Figure 9 we show the NME for the standard light-neutrino-exchange mechanism computed employing a variety of different numerical approaches. One can clearly see a significant variation of about a factor of ∼ 3 with some additional outliers corresponding to the rEDF (CDFT) approach. Given the distinct nature of individual nuclear structure computations the spread of the presented values clearly cannot be interpreted as theoretical uncertainty of the NMEs.
Reaching the estimated required accuracy on both LECs and NMEs seems to be rather challenging considering the current status of the relevant nuclear physics calculations. However, the recent advances in ab initio approaches to the computation of 0νββ decay NMEs seem to pave the path towards more reliable numerical values and clearer understanding of the theoretical uncertainties involved.

Other 0νββ Modes
Besides the usual 0νβ − β − -decay mode one could also make use of searches for neutrinoless modes of other ββ processes. In general there are four of these, ∆m ! > 0 (3.14) 2. β + β + : 3. ECβ + : 4. ECEC: Here, EC denotes electron capture. In principle, one could study all of these processes, as any of them would help to distinguish different mechanisms using the decay rate ratios considering nuclear uncertainties to be under control. However, all the other ββ processes listed in Table 6 are expected to have half-lives significantly longer than the usual 0νβ − β − decay and are therefore unlikely to show up in experiments. Despite that let us discuss their potential role in a bit more detail.

0νβ + β +
This process can be treated in a similar way as the usual 0νβ − β − decay, one only needs to consider a negative nuclear charge Z → −Z to calculate the positron wave functions. As such, the expected half-life will also be mainly determined by the PSF which goes with Q 5 . Looking at the second column of Table 6 we can see that the Q-values for naturally occurring isotopes are up to one order of magnitude smaller than usual 0νβ − β − Q-values. Additionally, the electromagnetic repulsion of the outgoing positrons deforms the wave functions and decreases the decay rate. Thus, we see that 0νβ + β + will be highly suppressed compared to 0νβ − β − . Also, given the similarities of the two decays there does not seem to be a natural way of enhancing the 0νβ + β + -decay rate with respect to 0νβ − β − . The relevant PSFs for the 0νβ + β + -decay of naturally occurring isotopes have been calculated in [78] and are about 3-5 orders of magnitude smaller than for 0νβ − β − decay.

0νECβ +
The PSFs for the neutrinoless mode of ECβ + were also calculated with good precision in [78] and are found to be 3-4 orders of magnitude smaller than those corresponding to 0νβ − β − decay. The reason is the same as in case of 0νβ + β + decay.

0νECEC
As this process has no particles other than the daughter isotope in the final state, mass degeneracy between the mother and daughter isotopes is required in order to satisfy conservation of energy and momentum. However, as the daughter isotope is a non-stable state due to the holes left in the electron shell after the electron capture, the corresponding decay width results in a resonance mechanism [79,80]. Resonances are often found when considering nuclear excitations in the final state isotope. However, the resonant enhancement strongly depends on the degeneracy between the initial and final state and hence small uncertainties in the mass measurements of these nuclei result in considerable uncertainties of the corresponding half-lives [79]. Existing studies tend to show that the resulting half-lives are still considerably longer than for 0νβ − β − [80]. Therefore, we do not consider this process in this work. However, it is fair to note that a close resonance might lead to half-lives comparable or even shorter than for 0νβ − β − decay. Recently, it was shown that further significant enhancement of the 0νECEC decay rate can be generated by a non-resonance shake mechanism [81]. In this case, the double electron capture is accompanied by emission of an electron from the shell of the final state isotope, which can carry away energy, thus making the whole process less dependent on the resonant behaviour. A dedicated review of the 0νECEC decay can be found in [82].

Bound-State 0νββ
Bound state 0νββ decay refers to a decay in which one or both of the two outgoing electrons end up in a bound energy level of the daughter isotope. It is usually referred to as 0νβEP and 0νEPEP for the one and two bound final state electrons, respectively, with EP denoting the electron production or electron placement. Being a reverse process of 0νECEC, also 0νEPEP requires mass degeneracy of the initial and final nuclei and the decay rate is described by a resonance-like mechanism. The explicit calculations show that the corresponding half-lives are even longer than those of double electron capture [79]. The reason is that there are no electron holes in the shell and only the (small) decay width of the nuclear excitation enters into the resonance. The single bound state double-β-decay 0νβEP was investigated in [83] and found to have PSFs 6-7 orders of magnitude smaller than those of 0νβ − β − decay. The decay rates can be significantly enhanced when considering fully ionized nuclei. In that case, the 0νβEP decay rate can for certain isotopes even exceed the one of 0νβ − β − decay [83]. Although this is an interesting idea, a full ionization of large number of isotopes represents an experimental challenge. Therefore, despite the enhanced decay rate, the number of available ions would be too small to reach the relevant experimental sensitivity.

Decay to excited final State Nuclei
Instead of utilizing different isotopes to determine the decay rate ratios one could also compare the ground state decay with the decay into an excited state final nucleus (0 + , 2 + ) using the same initial state isotope [84,85]. The potential benefit would be the possibility of studying this interplay within a single experiment. However, the excited state decays can be again expected to be highly suppressed due to the smaller phase-space resulting from the smaller Q-value. Additionally, previous studies tend to show that the NMEs for the decay into the excited final state are either of a similar size or smaller than those for the ground state decay [85][86][87], thus the half-lives would be rather further suppressed than enhanced by the nuclear part of the amplitude either.

Artificial Isotopes
Although there are 69 naturally occurring double-β-decaying isotopes, we found about ∼ 2700 possible 0νββ candidate isotopes when considering the full NIST list of elements [88]. Some of them have considerably larger Q-values of up to 50 MeV. 2 While such a large Qvalue of ∼ 50 MeV would result in a significant enhancement of the decay rate by ∼ 8 orders of magnitude, there are several fairly obvious experimental problems. Primarily, it is the artificial production of these isotopes, which would strongly limit the scale of the experiment. Again ton, kilogram and even gram scales would usually not be possible. Additionally, many artificial isotopes, especially those with large Q-values, come with additional decay modes that strongly dominate and often lead to extremely short half-lives such that storing them to study 0νββ decay would be impossible.
To sum up the above paragraphs, despite the fact that a variety of ββ processes exist, 0νβ − β − decay is largely the most relevant candidate to study, indeed. Therefore, other possible 0νββ modes would only become relevant in exotic scenarios leading either to their significant enhancements, or to strong suppression of 0νβ − β − decay. Given the similarity of all the ββ processes, such models would, however, seem to be rather unnatural from a particle physics point of view.

Distinguishing Specific Models
Following the discussion of possible discrimination among different LEFT operators, let us now have a brief look at complete models. As one would expect, lepton number violating BSM models will typically excite several LEFT operators at a time. While it would be challenging to identify a specific BSM model, as no finite set of BSM models exists and many different scenarios would result in the same low-energy physics, we do expect that, given fixed model parameters, one can at least check whether a model is consistent with the observed data and reject it if it is not. In the following paragraphs, we adopt and briefly discuss three different BSM scenarios that would lead to 0νββ decay. Each of the models will be compared with the standard mass mechanism predictions.

Minimal Left-Right Symmetric Model
The Standard Model is a chiral theory. That is, parity is explicitly broken due to the gauged SU (2) L symmetry and the missing right-handed neutrino. This particular choice of symmetries and particle content, additionally, results in vanishing neutrino masses. A simple approach to resolve these phenomena is to extend the Standard Model's gauge group to a left-right symmetric model SU (3) C × SU (2) L × SU (2) R × U (1) B−L [89][90][91] which is spontaneously broken to the Standard Model group SU (3) C × SU (2) L × U (1) Y . A comprehensive review of the minimal left-right symmetric Standard Model (mLRSM) is given in e.g. [92].
Extending the Standard Model to the left-right symmetric theory requires the existence of additional scalars and fermions. The conventional minimal setting includes two scalar triplets ∆ L ∈ (1, 3, 1, 2) and ∆ R ∈ (1, 1, 3, 2) as well as a scalar bidoublet Φ ∈ (1, 2, 2 * , 0) incorporating the SM Higgs doublet and the right-handed neutrinos ν R . The fermions are grouped into left-and right-handed doublets which under U ∈ SU (2) L,R transform as Figure 10. Feynman diagrams arising in the mLRSM that contribute to 0νββ. Here, ν i and N i represent the light and heavy neutrino mass eigenstates. It should be noted that, due to mixing of both left-and right-handed neutrinos and gauge bosons, each diagram (except the triplet-exchange diagram) comes with all possible combinations of the outgoing particle's chiralities. However, some diagrams are highly suppressed compared to others.
while the scalar fields transform as There are two discrete symmetries that one can impose onto a LR symmetric theory which can relate left-and right-handed fermions. These are parity P and charge conjugation C [93]. Thus, one can define two different discrete symmetry transformations Requiring either P or C invariance results in different constraints on the scalar potential as well as the Yukawa coupling matrices [93]. 3 The lepton number violation at low energy stems from the leptonic Yukawa interactions given by (4.7) 3 Note that a combination of both does not fit observational constraints [93].
After the neutral components of the scalars have acquired their VEVs one can infer the neutrino mass matrices from (4.7) (4.9) There are several diagrams in the mLRSM setting which can contribute to the 0νββ decay at tree level, see Figure 10. Detailed discussions of 0νββ decay within the mLRSM scenario can be found e.g. in [5,14,19]. The matching of the C-symmetric mLRSM onto SMEFT and, subsequently, onto the relevant LEFT operators has been discussed in [7]. Here, we will summarize their findings and study the distinguishability from the usual mass mechanism.
Integrating out the heavy fields with masses proportional to v R and matching the theory onto SMEFT results in the lepton number violating operators, where Φ SM is the Standard Model Higgs doublet. The matching scale corresponds to ∼ m W R and the Wilson coefficients at SMEFT level are given by LeudΦ , where v is the Standard Model Higgs doublets VEV, Here, C (5) corresponds to the usual seesaw formula. From the matching scale ∼ m W R the above coefficients have to be evolved down to m W ∼ 80 GeV, at which one can match onto the relevant LEFT operators by integrating out the remaining heavy particles with masses above m W . By doing so one obtains (4.13) Evolving the above coefficients down to the χPT scale of ∼ 2 GeV also generates a non-zero C (9) 5R coefficient since the RGEs of C 4,5 mix. The relevant Wilson coefficients are fixed by several physical parameters: the values of the triplet VEVs v L,R , the mass of the heavy right-handed triplet m ∆R as well as the masses of the three heavy neutrinos (m ν R1 , m ν R2 , m ν R3 ) and the lightest neutrino mass m ν min , the complex phases of the VEVs α and θ L and finally the left-right mixing parameter ξ. Here, we fix ξ = m b mt . The lightest neutrino mass together with the squared mass differences ∆m 2 ij that are known from oscillation data fix the remaining light neutrino masses for a given mass hierarchy. Taking we obtain Additionally, the mixing matrix U for the heavy neutrinos must be fixed. Here, we take U R = U PMNS for simplicity. In the C-symmetric case, one has (4.16) and the Dirac mass matrix can be derived as [94] as in [7], we can derive the LEFT Wilson coefficients in dependence on the minimal light neutrino mass m min , the Majorana phases entering U PMNS and the VEV phases θ L and α.
The resulting phase-space observables in this parameter setting of the mLRSM are hardly any different from the standard mechanism. This is because of the specific choice of Figure 11. Half-life ratios resulting from different mLRSM settings (different neutrino mass hierarchy and different minimal neutrino mass) when taking 76 Ge as the reference isotope. The ratios are compared to the standard mass mechanism. We vary both the unknown LECs as well as the unknown phases of the mLRSM model. parameters studied here which results in the scalar short-range contributions dominating over the long-range contributions. Hence, the phase-space is almost indistinguishable from the standard mechanism.
The resulting half-life ratios normalized to the standard mass mechanism are depicted in Figure 11. Here, additionally to varying the unknown LECs we also marginalized over the unknown phases. We can see that assuming inverted mass ordering there are only minor variations from the standard mechanism. In the case of normal ordering, the non-standard contributions alter the ratios notably only for small m min ≤ 10 −3 eV. In this region, as shown before in Figure 5, the central values of the variation differ significantly from the benchmark scenario. A similar behavior is manifested in Figure 12 displaying the half-life in dependence on the minimal neutrino mass m min for both orderings and in comparison with the standard mechanism on its own. One can see that in the case of inverted ordering the half-life is almost unaltered from the standard mechanism while for normal ordering the non-standard contributions start to play a substantial role below ∼ 10 −2 eV decreasing the expected half-life by about one order of magnitude compared to the standard scenario. In the same range of m min the uncertainties induced by the unknown LECs start to significantly influence the predicted half-life. On the other hand, the central values of the decay rate ratios alter for m min 10 −3 eV at most by a factor of R max im ββ ∼ 2.2 with 76 Ge as the reference isotope. The reason for this behaviour can be traced back to the dominance of the short-range contributions which (see Section 3) result in relatively small R im ββ despite the appearance of C  which relate the Standard Model's gauge bosons A a µ to their superpartner fermions Ψ a . One should note that since gauge bosons have 2 degrees of freedom (d.o.f.) and since this kind of transformation obviously cannot change the number of d.o.f., their superpartners Ψ a also have 2 degrees of freedom. Therefore, they are Majorana fermions. As particles within a supermultiplet must share the same mass, quantum numbers (except spin), interactions and couplings, SUSY must be broken at low energies to reproduce the experimentally confirmed SM predictions. Typically, after SUSY breaking there remains a discrete symmetry called R-parity (R p ) which can be assigned to every field, such that we have R p = +1 for Standard Model fields and R p = −1 for the superpartner fields. One can define R-parity as [96] R p = (−1) 2s+3(B−L) (4.21) where s is the spin and B and L are the corresponding baryon and lepton numbers of the field, respectively. If R p is a conserved quantity, it follows that the lightest superpartner cannot decay such that it becomes a candidate for explaining the origin of dark matter. However, R p conservation also comes with the conservation of both baryon and lepton number [96]. Thus, supersymmetric models aiming to explain the baryon asymmetry of the Universe via explicit violation of either lepton or baryon number need to break R p . This induces new lepton number violating terms [95] which can contribute to 0νββ decay. Contributions to 0νββ decay from R p −SUSY have been studied first in Refs. [95,97], the corresponding Feynman diagrams are shown in Figure 13. The relevant gluino (g) and neutralino (χ) -fermion interactions are given by [95,98,99] and One can obtain the low-energy effective Lagrangian by integrating out the heavy superfields as well as the Standard Model particles with masses m W . In doing so, one finds the different low-energy effective dimension-9 ∆L = 2 operators that contribute to 0νββ decay [97] (4.25) These can be matched onto the LEFT basis as (4.26) The coupling constants are given in terms of gluino, neutralino and squark masses as [95]  with Λ = √ 2π 3 (4.28) Both gluino-and neutralino-exchange diagrams contribute to the same low-energy operators. As pointed out in the previous section, distinguishing between the different contributions triggered by C 2R and C 3R is practically impossible due to the unknown LECs. Given that both operators contribute only to G 01 , the phase-space observables do not provide any additional information either. For completeness, we present the ratios normalized to the mass mechanism in Figure 14. Clearly, the R p -SUSY model we consider here follows the same pattern as the scalar short-range operators C (9) 2R and C (9) 2R already discussed in Section 3. In Figure 15 we show the expected half-life assuming the simultaneous existence of the standard mass mechanism and the R p -SUSY induced mechanisms. Here, we assume the masses of the non-neutralino superpartners to be given by the current experimental limits i.e. mẽ L = 410 GeV, mq L,R = 1600 GeV with q ∈ [u, d] and mg = 2260 GeV [100]. We fix the neutralino masses by requiring that the applied EFT framework holds, which necessitates m χ i ≥ Λ χ 2 GeV. For simplicity we take m χ 1 = 2 GeV and m χ i → ∞ for i = 1. Lighter neutralino masses in connection to 0νββ have also recently been studied [101]. We set the coupling constant to λ 111 = 2×10 −4 . Similarly to the mLRSM discussed above, the additional non-standard contributions hardly affect the inverted ordering setting. However, in the normal ordering case the non-standard contributions from R p -SUSY model start to significantly influence the expected half-lives decreasing them, again, by about one order of magnitude. While one would naively assume that this should result in significantly enhanced ratios, it is important to bear in mind that any enhancement in the decay rates which is independent of the isotope of interest will drop out when considering the decay rate ratios.

Leptoquark Models
Leptoquarks (LQs) are hypothetical bosons (3, X, Y ) with non-zero color charge which couple to both quarks and leptons. They appear in numerous Standard Model extensions such as technicolor and composite models [102,103] or grand unifications [104,105] and can be used to generate neutrino masses at 1-loop level [106]. For a comprehensive review on leptoquarks see e.g. [107].
Ignoring leptoquarks which do not directly couple to the Standard Model's particle content (i.e. without right-handed neutrinos), one can add up to 10 different leptoquarks obeying the Standard Model symmetries [108]. These are summarized in Table 5. By looking at the relevant Feynman diagrams in Figure 16 we can see that the contributions to 0νββ decay arise from leptoquarks with Q (1) = ±1/3 (Figure 16 left) and Q (2) = ±2/3 (Figure 16 right). The full set of renormalizable LQ-fermion interactions is given by [108] L and for the scalar (S) and vector (V) leptoquarks, respectively. We follow the notation of [108] distinguishing leptoquarks coupling to left-handed and right-handed quarks. In addition to the LQ-fermion interactions, one can write down the gauge invariant and renormalizable LQ-Higgs interactions, where the leptoquark triplets are defined as These LQ-Higgs interactions are essential when considering contributions to 0νββ decay because they result in non-zero correlation functions for, e.g., 33) where N is the mixing matrix which diagonalizes the mass matrix N T M 2 N = M 2 diag and I = N T I are the mass eigenstate fields. This particular example results in contributions captured by the right diagram in Figure 16.
After integrating out the heavy LQ degrees of freedom and rearranging the resulting EFT operators via Fierz transformations one arrives at the effective low-energy four-fermion interactions. The parts of the low-energy Lagrangian relevant for 0νββ decay are then given by [108] with the low-energy Wilson coefficients (4.37) Here, "common mass scales" M S and M V have been inserted for convenience. It should be noted that the exact choice of M S,V does not matter as they drop out. However, the exact LQ masses do enter into the calculation such that for leptoquark masses which are about the same order of magnitude one can choose M S,V to represent the suppression factors. Looking (4.38) We study the following 7 different settings of LQ contributions to 0νββ decay: 1. Full LQ Model: The left-handed scalar (SL) and left-handed vector (VL) models result in the same low-energy physics because they match onto the same LEFT operator. The same is true for SR and VR. In Figure 17 we show the corresponding single electron energy spectra and angular correlations corresponding to each of the above models and compare them with the standard mechanism scenario. When setting the unknown LECs to their orderof-magnitude estimates we find that except for the vector (V) scenario all other models give shapes distinguishable from the standard mass mechanism for at least one phase-space observable. The resulting half-life ratios normalized to the neutrino mass mechanism for each of the above scenarios are shown in Figure 18. Except for the SR and VR cases, for which the central values suggest somewhat weaker distinguishability, we find that the central values match fairly well the chosen benchmark scenario. Nonetheless, the spread in R im ββ is still significant for the full model as well as the SL and VL models. Considering the central values, the highest ratio when taking 76 Ge as the reference isotope is realized in the vector model with R max im ββ ∼ 4.5. Again, assuming that the calculated central values of the half-life ratios represent a reasonable estimate, this would correspond to a necessary theoretical accuracy on the nuclear part of the amplitude to satisfy ∆M eff M eff

19%.
In Figure 19 we show the expected half-lives for the simultaneous realization of the full LQ model and the standard mass mechanism. We assumed the suppression factors to be M S = M V = 10 7 GeV. One can see that in this setting the inverted mass ordering case is not altered significantly while the half-life in the normal ordering case is decreased such that the gap between the two mass orderings is closed.

Summary and Conclusions
Neutrinoless double beta decay is the best laboratory probe of lepton number violation and as such can naturally shed light on the generation of neutrino masses as well as associated UV physics. The implications of observation of this hypothetical nuclear process would largely depend on the mechanism responsible for the dominant contribution. In this paper we have performed a detailed analysis discussing the possibilities of experimental discrimination among the 32 different LEFT LNV operators of dimension ≤ 9 triggering 0νββ decay at low energy.
The main aim of our study is to understand the differences in various 0νββ decay mechanisms and to investigate the possibilities of identifying the potential exotic contribution in experiments. Assuming only one operator at a time, we found that the 32 different LEFT operators can be split into 12 groups which are distinguishable from each other by comparison of ratios of half-lives in different double-beta-decaying isotopes. We calculated the half-life ratios normalized to the standard mass mechanism R im ββ for each of the operator groups discussing the potential for their identification by experimental observations. Varying the currently unknown low-energy constants (LECs) around their order-of-magnitude estimates obtained using NDA we observed that their impact on the expected half-life ratios can be significant for most operator groups. To quantify this impact and temporarily eliminate it in our conclusions we focused on two different scenarios; namely, we identified the central values of the ratio ranges as well as the worst-case scenario considering the value of the ratios closest to 1 within each ratio range. In Figure 8 we summarized the potential of distinguishing among different operator groups for both of these scenarios considering all isotopes for which the experimental limit on 0νββ decay half-life exists. Based on the central-value scenario we estimated the required theoretical accuracy on the nuclear physics calculations, parameterized by an effective nuclear matrix element M eff , that would allow for identifying non-standard mechanisms in the single operator dominance scenario. We found that identifying all the non-standard mechanisms via half-life ratios would require (at least) a few-percent accuracy. While the required accuracy of the theoretical description is beyond the current status of nuclear uncertainties, advances in ab initio calculations of nuclear matrix elements may be able to deliver such precision in the future provided that the currently unknown LECs are fixed with similar accuracy as those that are already under control.
The additional information that can be inferred from the phase-space observables does not allow for distinguishing operators within the 12 operator groups corresponding to distinct half-life ratios. However, the phase-space observables are much less affected by nuclear uncertainties such that they can potentially deliver important insight into the underlying mechanism of 0νββ decay even if nuclear uncertainties remain significant. For operator groups such as C V (9) orC V (9) for which the expected half-life ratios do not differ significantly from the standard mass mechanism measurements, tracking the outgoing electrons would be a more promising approach of identification even if nuclear uncertainties are substantially reduced. Therefore, future experiments that would have the required technology such as SuperNEMO seem to be very relevant. The operator groups that could be distinguished by means of phase-space observables are also marked in Figure 8.
Besides the effective approach detailing individual operators, we focussed also on 0νββ decay contributions triggered by three different high-energy models. In each case, we identified and discussed the signatures that could help to distinguish these models in observations of 0νββ decay. Our approach can be easily extended to any other UV model that can be matched onto the applied EFT framework.
Based on the obtained results and following their discussion it becomes clear that although it might be possible to unravel an exotic contribution to 0νββ decay, pinpointing the dominant mechanism underlying this hypothetical nuclear process most probably would not be possible without other, complementary experiments. As discussed, the possibilities of employing different double-beta processes seems to be rather unlikely because of their phase-space suppression. On the other hand, the underlying mechanism could be identified by combining the 0νββ decay data with different experiments searching for lepton number violation, such as meson decays, tau decays, or collider searches, which can, however, more naturally verify a complete UV scenario rather than a specific effective operator. Useful information will be provided also by measurements aiming to determine the absolute neutrino mass scale, including the CMB data providing a constraint on the sum of neutrino masses, i m ν i . A detailed discussion of the interplay with complementary probes of neutrino masses and lepton number (non-)conservation is however beyond the scope of this paper. Similarly, we have not covered the possibility of existence of light sterile neutrinos and its implications for 0νββ decay. Same methods as employed in this study may allow -38 -to unravel additional contributions to 0νββ decay induced by light sterile neutrinos.
In Table 6 we present a list of all naturally occuring isotopes that decay via any of the 0νββ-modes i.e. 0νβ − β − , 0νβ + β + , 0νECβ + and 0νECEC. The isotopes as well as the corresponding Q-values are taken and calculated from the NIST list of elements [88].

B Contributions from each Operator
Assuming only one non-vanishing operator at a time, the half-life can be written in terms of a single Wilson coefficient, different phase-space factors and the nuclear contributions determined by the different NMEs and LECs. For convenience, we list the explicit decay rate equations for each operator below.

C Considering all Isotopes
While we have focussed our discussion on isotopes for which experimental limits on the halflives exist, we want to present our main findings of Figure 8 here again but now considering all naturally occuring 0νββ isotopes for which we have nuclear matrix elements available in the IBM2 framework. The corresponding results are presented in Figure 20. In Figures 21  and 22 we show the resulting ratios including variations of the unknown LECs similar to Figures 5 and 6 when considering the whole set of isotopes available.