aainstitutetext: Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, 69117
Heidelberg, Germanybbinstitutetext: Department of Physics, University of
California, Berkeley, CA 94720, USAccinstitutetext: Department of Physics,
University of California, San Diego, CA 92093, USA

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

Lukáš Gráf<EMAIL_ADDRESS>a Manfred Lindner<EMAIL_ADDRESS>a
Oliver Scholer<EMAIL_ADDRESS>
###### Abstract

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.

###### Contents

  1. 1 Introduction
  2. 2 EFT Approach to $0\nu\beta\beta$ Decay: The Master Formula
    1. 2.1 The half-life master formula
    2. 2.2 Sub-Amplitudes
    3. 2.3 Relation to literature
  3. 3 Distinguishing the Effective Operators
    1. 3.1 Phase-Space Observables
    2. 3.2 Decay Rate Ratios
      1. 3.2.1 Sensitivity on the unknown LECs
      2. 3.2.2 Distinguishing different operators
      3. 3.2.3 The impact of nuclear uncertainties
    3. 3.3 Other $\mathbf{0\nu\beta\beta}$ Modes
      1. 3.3.1 $\mathbf{0\nu\beta^{+}\beta^{+}}$
      2. 3.3.2 $\mathbf{0\nu EC\beta^{+}}$
      3. 3.3.3 $\mathbf{0\nu ECEC}$
      4. 3.3.4 Bound-State $\mathbf{0\nu\beta\beta}$
      5. 3.3.5 Decay to excited final State Nuclei
      6. 3.3.6 Artificial Isotopes
  4. 4 Distinguishing Specific Models
    1. 4.1 Minimal Left-Right Symmetric Model
    2. 4.2 Gluino- and Neutralino-exchange in $\cancel{R}_{p}$ \- SUSY
    3. 4.3 Leptoquark Models
  5. 5 Summary and Conclusions
  6. A Double $\beta$ Modes
  7. B Contributions from each Operator
  8. C Considering all Isotopes

##  1 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\nu\beta\beta$) 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\nu\beta\beta$ decay [3, 4, 5, 6, 7, 8, 9]. In fact, while
the sole observation of $0\nu\beta\beta$ 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\nu\beta\beta$ 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\nu\beta\beta$ decay [14, 15, 16, 17, 18, 19, 20].

There is variety of experiments searching for $0\nu\beta\beta$ decay in
different double-beta-decaying isotopes [21, 22, 23, 24, 25, 26, 27, 28, 29,
30]. Currently, the best limit on the half-life reaching $2.3\times
10^{26}\,$years, is claimed by the KamLAND-Zen collaboration [31] studying the
decay of 136Xe. The most stringent limit on the half-life of $0\nu\beta\beta$
decay of 76Ge attains $1.8\times 10^{26}\,$years, as obtained by the GERDA
collaboration [21]. Proposed next generation experiments such as LEGEND [32,
33] (76Ge), CUPID [34] (100Mo), SNO+ [35] (130Te) and nEXO [36] (136Xe) 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\nu\beta\beta$ decay, but even in standard double beta
decay [37]. Recent reviews of the experimental and theoretical efforts in the
field of $0\nu\beta\beta$ 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\nu\beta\beta$ 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\nu\beta\beta$ 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\nu\beta\beta$ Decay: The Master Formula

###  2.1 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\nu\beta\beta$ decay
amplitude can be described in terms of an $SU(3)_{C}\times U(1)_{Q}$ invariant
low-energy effective field theory (LEFT). Including operators up to LEFT
dimension 9 the most relevant Lagrangians for $0\nu\beta\beta$ are given by
[6, 7]

$\displaystyle\begin{split}\mathcal{L}_{\Delta
L=2}^{(6)}=\frac{2G_{F}}{\sqrt{2}}\bigg{[}&C_{\text{VL}}^{(6)}\;\big{(}\overline{u_{L}}\gamma^{\mu}d_{L}\big{)}\;\big{(}\overline{e_{R}}\gamma_{\mu}\nu^{c}_{L}\big{)}+C_{\text{VR}}^{(6)}\;\big{(}\overline{u_{R}}\gamma^{\mu}d_{R}\big{)}\;\big{(}\overline{e_{R}}\gamma_{\mu}\nu^{c}_{L}\big{)}\\\
+&C_{\text{SL}}^{(6)}\;\big{(}\overline{u_{R}}d_{L}\big{)}\;\big{(}\overline{e_{L}}\nu^{c}_{L}\big{)}+C_{\text{SR}}^{(6)}\;\big{(}\overline{u_{L}}d_{R}\big{)}\;\big{(}\overline{e_{L}}\nu^{c}_{L}\big{)}\\\
+&C_{\text{T}}^{(6)}\;\big{(}\overline{u_{L}}\sigma^{\mu\nu}d_{R}\big{)}\;\big{(}\overline{e_{L}}\sigma_{\mu\nu}\nu^{c}_{L}\big{)}\bigg{]}+\text{h.c.}\end{split}$
(1)

and

$\displaystyle\begin{split}\mathcal{L}_{\Delta
L=2}^{(7)}=\frac{2G_{F}}{\sqrt{2}v}\bigg{[}&C_{\text{VL}}^{(7)}\;\big{(}\overline{u_{L}}\gamma^{\mu}d_{L}\big{)}\;\big{(}\overline{e_{L}}\overset{\leftrightarrow}{\partial}_{\mu}\nu^{c}_{L}\big{)}+C_{\text{VR}}^{(7)}\;\big{(}\overline{u_{R}}\gamma^{\mu}d_{R}\big{)}\;\big{(}\overline{e_{L}}\overset{\leftrightarrow}{\partial}_{\mu}\nu^{c}_{L}\big{)}\bigg{]}+\text{h.c.}\end{split}$
(2)

for the long-range part, where

$\displaystyle\alpha\overset{\leftrightarrow}{\partial}\beta=\alpha(\partial\beta)-(\partial\alpha)\beta.$
(3)

as well as the dimension 9 short-range Lagrangian

$\displaystyle\begin{split}\mathcal{L}_{\Delta
L=2}^{(9)}=\frac{1}{v^{5}}\sum_{i}\bigg{[}\bigg{(}&C_{i,R}^{(9)}\;\big{(}\overline{e_{R}}e_{R}^{c}\big{)}+C_{i,L}^{(9)}\;\big{(}\overline{e_{L}}e_{L}^{c}\big{)}\bigg{)}\;\mathcal{O}_{i}+C_{i}^{(9)}\;\big{(}\overline{e}\gamma_{\mu}\gamma_{5}e^{c}\big{)}\;\mathcal{O}_{i}^{\mu}\bigg{]}\end{split}$
(4)

with the scalar $\mathcal{O}_{i}$ and vector $\mathcal{O}_{i}^{\mu}$ four-
quark operators [41, 7]

$\displaystyle\begin{aligned}
\mathcal{O}_{1}&=\big{(}\overline{u_{L}}^{\alpha}\gamma_{\mu}d_{L}^{\alpha}\big{)}\;\big{(}\overline{u_{L}}^{\beta}\gamma^{\mu}d_{L}^{\beta}\big{)}\,,\qquad&\mathcal{O}_{1}{{}^{\prime}}&=\big{(}\overline{u_{R}}^{\alpha}\gamma_{\mu}d_{R}^{\alpha}\big{)}\;\big{(}\overline{u_{R}}^{\beta}\gamma^{\mu}d_{R}^{\beta}\big{)}\,,\\\
\mathcal{O}_{2}&=\big{(}\overline{u_{R}}^{\alpha}d_{L}^{\alpha}\big{)}\;\big{(}\overline{u_{R}}^{\beta}d_{L}^{\beta}\big{)}\,,\qquad&\mathcal{O}_{2}{{}^{\prime}}&=\big{(}\overline{u_{L}}^{\alpha}d_{R}^{\alpha}\big{)}\;\big{(}\overline{u_{L}}^{\beta}d_{R}^{\beta}\big{)}\,,\\\
\mathcal{O}_{3}&=\big{(}\overline{u_{R}}^{\alpha}d_{L}^{\beta}\big{)}\;\big{(}\overline{u_{R}}^{\beta}d_{L}^{\alpha}\big{)}\,,\qquad&\mathcal{O}_{3}{{}^{\prime}}&=\big{(}\overline{u_{L}}^{\alpha}d_{R}^{\beta}\big{)}\;\big{(}\overline{u_{L}}^{\beta}d_{R}^{\alpha}\big{)}\,,\\\
\mathcal{O}_{4}&=\big{(}\overline{u_{L}}^{\alpha}\gamma_{\mu}d_{L}^{\alpha}\big{)}\;\big{(}\overline{u_{R}}^{\beta}\gamma^{\mu}d_{R}^{\beta})\,,\\\
\mathcal{O}_{5}&=\big{(}\overline{u_{L}}^{\alpha}\gamma_{\mu}d_{L}^{\beta}\big{)}\;\big{(}\overline{u_{R}}^{\beta}\gamma^{\mu}d_{R}^{\alpha}\big{)}\,,\\\
\\\
\mathcal{O}_{6}^{\mu}&=\big{(}\overline{u_{L}}\gamma^{\mu}d_{L}\big{)}\;\big{(}\overline{u_{L}}d_{R}\big{)}\,,\qquad&\mathcal{O}_{6}^{\mu}{{}^{\prime}}&=\big{(}\overline{u_{R}}\gamma^{\mu}d_{R}\big{)}\;\big{(}\overline{u_{R}}d_{L}\big{)}\,,\\\
\mathcal{O}_{7}^{\mu}&=\big{(}\overline{u_{L}}t^{A}\gamma^{\mu}d_{L}\big{)}\;\big{(}\overline{u_{L}}t^{A}d_{R}\big{)}\,,\qquad&\mathcal{O}_{7}^{\mu}{{}^{\prime}}&=\big{(}\overline{u_{R}}t^{A}\gamma^{\mu}d_{R}\big{)}\;\big{(}\overline{u_{R}}t^{A}d_{L}\big{)}\,,\\\
\mathcal{O}_{8}^{\mu}&=\big{(}\overline{u_{L}}\gamma^{\mu}d_{L}\big{)}\;\big{(}\overline{u_{R}}d_{L}\big{)}\,,\qquad&\mathcal{O}_{8}^{\mu}{{}^{\prime}}&=\big{(}\overline{u_{R}}\gamma^{\mu}d_{R}\big{)}\;\big{(}\overline{u_{L}}d_{R}\big{)}\,,\\\
\mathcal{O}_{9}^{\mu}&=\big{(}\overline{u_{L}}t^{A}\gamma^{\mu}d_{L}\big{)}\;\big{(}\overline{u_{R}}t^{A}d_{L}\big{)}\,,\qquad&\mathcal{O}_{9}^{\mu}{{}^{\prime}}&=\big{(}\overline{u_{R}}t^{A}\gamma^{\mu}d_{R}\big{)}\;\big{(}\overline{u_{L}}t^{A}d_{R}\big{)}\,.\end{aligned}$
(5)

Here $\alpha,\beta$ are color-indices and the $t^{A}$ are the generators of
SU(3) in the fundamental representation given by the 8 Gell-Mann matrices
$\lambda^{A}$ as $t^{A}=\frac{1}{2}\lambda^{A},A=1...8$. The operators
$\mathcal{O}$ and $\mathcal{O}{{}^{\prime}}$ in (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\nu\beta\beta$ decay.

The transition from the quark level to the nuclear level can be achieved
employing the chiral effective field theory ($\chi$EFT) [42]. The expected
half-life contributed by the 32 effective operators is then captured by a
“$0\nu\beta\beta$ master-formula” combining the 32 LEFT Wilson 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 $\chi$EFT enter the
formula - we summarize these in Table 3. The $0\nu\beta\beta$ half-life is
then given in terms of different sub-amplitudes $\mathcal{A}_{i}$ as

$\displaystyle\begin{split}\left(T_{1/2}^{0\nu}\right)^{-1}=&g_{A}^{4}\bigg{[}G_{01}\left(\left|\mathcal{A}_{\nu}\right|^{2}+\left|\mathcal{A}_{R}\right|^{2}\right)-2\left(G_{01}-G_{04}\right)\real\left[\mathcal{A}_{\nu}^{*}\mathcal{A}_{R}\right]\\\
&+4G_{02}\left|\mathcal{A}_{E}\right|^{2}+2G_{04}\left(\left|\mathcal{A}_{m_{e}}\right|^{2}+\real\left[\mathcal{A}_{m_{e}}^{*}\left(\mathcal{A}_{\nu}+\mathcal{A}_{R}\right)\right]\right)\\\
&-2G_{03}\real\left[\left(\mathcal{A}_{\nu}+\mathcal{A}_{R}\right)\mathcal{A}_{E}^{*}+2\mathcal{A}_{m_{e}}\mathcal{A}_{E}^{*}\right]\\\
&+G_{09}\left|\mathcal{A}_{M}\right|^{2}+G_{06}\real\left[\left(\mathcal{A}_{\nu}-\mathcal{A}_{R}\right)\mathcal{A}_{M}^{*}\right]\bigg{]}\end{split}$
(6)

| $G_{01}$ | $G_{02}$ | $G_{03}$ | $G_{04}$ | $G_{06}$ | $G_{09}$
---|---|---|---|---|---|---
238U | 6.96 | 3.79 | 2.75 | 5.26 | 14.43 | 17.32
232Th | 2.70 | 0.73 | 0.76 | 1.83 | 6.35 | 7.14
198Pt | 1.23 | 0.55 | 0.44 | 0.90 | 2.64 | 3.10
160Gd | 1.33 | 1.68 | 0.73 | 1.12 | 2.25 | 3.08
154Sm | 0.44 | 0.28 | 0.18 | 0.34 | 0.87 | 1.08
150Nd | 8.82 | 40.15 | 7.00 | 8.25 | 9.83 | 18.78
148Nd | 1.36 | 2.10 | 0.79 | 1.17 | 2.15 | 3.09
136Xe | 1.88 | 4.64 | 1.26 | 1.69 | 2.58 | 4.14
134Xe | 0.08 | 0.02 | 0.02 | 0.05 | 0.18 | 0.20
130Te | 1.81 | 4.68 | 1.22 | 1.63 | 2.43 | 3.96
128Te | 0.07 | 0.02 | 0.02 | 0.05 | 0.17 | 0.19
124Sn | 1.13 | 2.42 | 0.72 | 1.01 | 1.62 | 2.51
116Cd | 2.06 | 6.51 | 1.46 | 1.89 | 2.59 | 4.47
110Pd | 0.58 | 0.95 | 0.34 | 0.50 | 0.89 | 1.30
100Mo | 1.89 | 6.80 | 1.36 | 1.75 | 2.25 | 4.06
96Zr | 2.42 | 10.43 | 1.81 | 2.26 | 2.68 | 5.15
82Se | 1.15 | 3.96 | 0.80 | 1.06 | 1.37 | 2.47
76Ge | 0.26 | 0.43 | 0.15 | 0.23 | 0.40 | 0.59
Table 1: The different PSFs in terms of $10^{-14}y^{-1}$ used in our
calculations.

###  2.2 Sub-Amplitudes

The sub-amplitudes $\mathcal{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 written as

$\displaystyle\begin{split}\mathcal{A}_{\nu}=&\frac{m_{\beta\beta}}{m_{e}}\mathcal{M}^{(3)}_{\nu}+\frac{m_{N}}{m_{e}}\mathcal{M}^{(6)}_{\nu}\left(C_{\text{SL}}^{(6)},C_{\text{SR}}^{(6)},C_{\text{T}}^{(6)},C_{\text{VL}}^{(7)},C_{\text{VR}}^{(7)}\right)\\\
&+\frac{m^{2}_{N}}{m_{e}v}\mathcal{M}^{(9)}_{\nu}\left(C_{\text{1L}}^{(9)},C_{\text{1L}}^{(9)}{{}^{\prime}},C_{\text{2L}}^{(9)},C_{\text{2L}}^{(9)}{{}^{\prime}},C_{\text{3L}}^{(9)},C_{\text{3L}}^{(9)}{{}^{\prime}},C_{\text{4L}}^{(9)},C_{\text{5L}}^{(9)}\right)\,,\end{split}$
$\displaystyle\begin{split}\mathcal{A}_{R}=&\frac{m^{2}_{N}}{m_{e}v}\mathcal{M}^{(9)}_{R}\left(C_{\text{1R}}^{(9)},C_{\text{1R}}^{(9)}{{}^{\prime}},C_{\text{2R}}^{(9)},C_{\text{2R}}^{(9)}{{}^{\prime}},C_{\text{3R}}^{(9)},C_{\text{3R}}^{(9)}{{}^{\prime}},C_{\text{4R}}^{(9)},C_{\text{5R}}^{(9)}\right)\,,\end{split}$
$\displaystyle\begin{split}\mathcal{A}_{E}=&\mathcal{M}^{(6)}_{E,L}\left(C_{\text{VL}}^{(6)}\right)+\,\mathcal{M}^{(6)}_{E,R}\left(C_{\text{VR}}^{(6)}\right)\,,\end{split}$
$\displaystyle\begin{split}\mathcal{A}_{m_{e}}=&\mathcal{M}^{(6)}_{m_{e},L}\left(C_{\text{VL}}^{(6)}\right)+\,\mathcal{M}^{(6)}_{m_{e},R}\left(C_{\text{VR}}^{(6)}\right)\,,\end{split}$
$\displaystyle\begin{split}\mathcal{A}_{M}=&\frac{m_{N}}{m_{e}}\,\mathcal{M}^{(6)}_{M}\left(C_{\text{VL}}^{(6)}\right)+\frac{m_{N}^{2}}{m_{e}v}\mathcal{M}^{(9)}_{M}\left(C_{\text{6}}^{(9)},C_{\text{6}}^{(9)}{{}^{\prime}},C_{\text{7}}^{(9)},C_{\text{7}}^{(9)}{{}^{\prime}},C_{\text{8}}^{(9)},C_{\text{8}}^{(9)}{{}^{\prime}},C_{\text{9}}^{(9)},C_{\text{9}}^{(9)}{{}^{\prime}}\right)\,.\end{split}$
(7)

The matrix elements $\mathcal{M}_{i}$ depend on the different LECs and Wilson
coefficients. We explicitly state the dependency on the different Wilson
coefficients within the brackets in (7).

| $M_{F}$ | $M_{GT}^{AA}$ | $M_{GT}^{AP}$ | $M_{GT}^{PP}$ | $M_{GT}^{MM}$ | $M_{T}^{AA}$ | $M_{T}^{AP}$ | $M_{T}^{PP}$ | $M_{T}^{MM}$ | $M_{Fsd}$ | $M^{AA}_{GTsd}$ | $M^{AP}_{GTsd}$ | $M^{PP}_{GTsd}$ | $M^{AP}_{Tsd}$ | $M^{PP}_{Tsd}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
76Ge | -0.78 | 6.06 | -0.86 | 0.17 | 0.20 | 0.0 | 0.24 | -0.06 | 0.04 | -1.20 | 4.18 | -1.24 | 0.29 | -0.77 | 0.23
82Se | -0.67 | 4.93 | -0.71 | 0.14 | 0.17 | 0.0 | 0.24 | -0.06 | 0.04 | -1.01 | 3.46 | -1.03 | 0.25 | -0.73 | 0.22
96Zr | -0.36 | 4.32 | -0.64 | 0.13 | 0.15 | 0.0 | -0.21 | 0.05 | -0.04 | -0.87 | 3.06 | -0.89 | 0.21 | 0.64 | -0.20
100Mo | -0.51 | 5.55 | -0.90 | 0.20 | 0.22 | 0.0 | -0.29 | 0.07 | -0.05 | -1.28 | 4.48 | -1.33 | 0.30 | 0.93 | -0.28
110Pd | -0.42 | 4.43 | -0.76 | 0.17 | 0.18 | 0.0 | -0.21 | 0.06 | -0.04 | -1.07 | 3.72 | -1.11 | 0.25 | 0.79 | -0.24
116Cd | -0.34 | 3.17 | -0.55 | 0.12 | 0.13 | 0.0 | -0.12 | 0.04 | -0.03 | -0.80 | 2.72 | -0.81 | 0.18 | 0.49 | -0.16
124Sn | -0.57 | 3.37 | -0.50 | 0.11 | 0.12 | 0.0 | 0.12 | -0.03 | 0.02 | -0.82 | 2.56 | -0.77 | 0.19 | -0.42 | 0.13
128Te | -0.72 | 4.32 | -0.64 | 0.13 | 0.15 | 0.0 | 0.12 | -0.04 | 0.03 | -1.03 | 3.24 | -0.98 | 0.24 | -0.52 | 0.16
130Te | -0.65 | 3.89 | -0.57 | 0.12 | 0.14 | 0.0 | 0.14 | -0.04 | 0.02 | -0.94 | 2.95 | -0.89 | 0.22 | -0.47 | 0.15
134Xe | -0.69 | 4.21 | -0.62 | 0.13 | 0.15 | 0.0 | 0.12 | -0.04 | 0.03 | -0.97 | 3.07 | -0.92 | 0.22 | -0.48 | 0.15
136Xe | -0.52 | 3.20 | -0.45 | 0.09 | 0.11 | 0.0 | 0.12 | -0.03 | 0.02 | -0.73 | 2.32 | -0.69 | 0.17 | -0.36 | 0.12
148Nd | -0.36 | 2.52 | -0.48 | 0.11 | 0.12 | 0.0 | -0.12 | 0.02 | -0.02 | -0.78 | 2.54 | -0.79 | 0.19 | 0.30 | -0.09
150Nd | -0.51 | 3.75 | -0.76 | 0.17 | 0.19 | 0.0 | -0.12 | 0.04 | -0.03 | -0.74 | 2.46 | -0.76 | 0.18 | 0.34 | -0.10
154Sm | -0.34 | 2.98 | -0.52 | 0.11 | 0.13 | 0.0 | -0.12 | 0.03 | -0.02 | -0.78 | 2.64 | -0.79 | 0.19 | 0.39 | -0.13
160Gd | -0.42 | 4.22 | -0.71 | 0.15 | 0.17 | 0.0 | -0.21 | 0.05 | -0.03 | -1.02 | 3.52 | -1.04 | 0.24 | 0.60 | -0.19
198Pt | -0.33 | 2.27 | -0.50 | 0.11 | 0.12 | 0.0 | -0.12 | 0.03 | -0.02 | -0.78 | 2.57 | -0.78 | 0.18 | 0.37 | -0.12
232Th | -0.44 | 4.17 | -0.76 | 0.17 | 0.18 | 0.0 | -0.21 | 0.05 | -0.04 | -1.08 | 3.80 | -1.11 | 0.25 | 0.69 | -0.22
238U | -0.52 | 4.96 | -0.90 | 0.20 | 0.21 | 0.0 | -0.21 | 0.06 | -0.04 | -1.29 | 4.51 | -1.32 | 0.30 | 0.82 | -0.25

Table 2: NMEs used in our calculations based on the IBM2 model [43]

$\mathcal{A}_{\nu}$ depends on the matrix elements

$\displaystyle\begin{split}\mathcal{M}_{\nu}^{(3)}&=-V_{ud}^{2}\left(-\frac{1}{g_{A}^{2}}M_{F}+\mathcal{M}_{GT}+\mathcal{M}_{T}+2\frac{m_{\pi}^{2}\mathbf{g_{\nu}^{NN}}}{g_{A}^{2}}M_{F,sd}\right),\\\
\mathcal{M}_{\nu}^{(6)}&=V_{ud}\left(\frac{B}{m_{N}}\left(C_{SL}^{(6)}-C_{SR}^{(6)}\right)+\frac{m_{\pi}^{2}}{m_{N}v}\left(C_{VL}^{(7)}-C_{VR}^{(7)}\right)\right)\mathcal{M}_{PS}+V_{ud}C_{\text{T}}^{(6)}\mathcal{M}_{T6},\\\
\mathcal{M}_{\nu}^{(9)}&=-\frac{1}{2m_{N}^{2}}C_{\pi\pi
L}^{(9)}\left(M_{GT,sd}^{AP}+M_{T,sd}^{AP}\right)-\frac{2m_{\pi}^{2}}{g_{A}^{2}m_{N}^{2}}C_{NNL}^{(9)}M_{F,sd},\end{split}$
(8)

where $\mathcal{M}_{\nu}^{(3)}$ 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_{\nu}^{NN}$, which parametrizes the contact-term
contribution originating from the exchange of hard neutrinos [44, 45].
$\mathcal{A}_{R}$ is given by

Known LECs |  | Unknown LECs
---|---|---
$g_{A}$ | $1.271$ |  | $\left|g_{T}^{\prime}\right|$ | $\mathcal{O}(1)$
$g_{S}$ |  $0.97$ [46] |  | $\left|g_{T}{\pi\pi}\right|$ | $\mathcal{O}(1)$
$g_{M}$ | $4.7$ |  | $\left|g_{1,6,7,8,9}^{\pi N}\right|$ | $\mathcal{O}(1)$
$g_{T}$ |  $0.99$ [46] |  | $\left|g_{VL}^{\pi N}\right|$ | $\mathcal{O}(1)$
$B$ | $2.7\,\mathrm{GeV}$ |  | $\left|g_{T}^{\pi N}\right|$ | $\mathcal{O}(1)$
$g_{1}^{\pi\pi}$ |  $0.36$ [47] |  | $\left|g_{1,6,7}^{NN}\right|$ | $\mathcal{O}(1)$
$g_{2}^{\pi\pi}$ |  $2.0$ [47] |  | $\left|g_{2,3,4,5}^{NN}\right|$ | $\mathcal{O}(16\pi^{2})$
$g_{3}^{\pi\pi}$ |  $-0.62$ [47] |  | $\left|g_{VL}^{NN}\right|$ | $\mathcal{O}(1)$
$g_{4}^{\pi\pi}$ |  $-1.9$ [47] |  | $\left|g_{T}^{NN}\right|$ | $\mathcal{O}(1)$
$g_{5}^{\pi\pi}$ |  $-8.0$ [47] |  | $\left|g_{VL,VR}^{E,m_{e}}\right|$ | $\mathcal{O}(1)$
$g_{\nu}^{NN}$ |  $-92.9\,\mathrm{GeV^{-2}}\pm 50\%$ [48, 49, 50] |  |  |
Table 3: Summary of the low-energy constants necessary to calculate the
$0\nu\beta\beta$ half-life for all 32 different operators. The table is taken
from [7] and restructured.

$\displaystyle\begin{split}\mathcal{M}_{R}^{(9)}&=\mathcal{M}_{\nu}^{(9)}|_{L\rightarrow
R},\end{split}$ (9)

for $\mathcal{A}_{E}$ the different contributions are

$\displaystyle\begin{split}\mathcal{M}_{E,L}^{(6)}&=-\frac{V_{ud}C_{VL}^{(6)}}{3}\left(\frac{g_{V}^{2}}{g_{A}^{2}}M_{F}+\frac{1}{3}\left(2M_{GT}^{AA}+M_{T}^{AA}\right)+\frac{6\mathbf{g_{V\;L}^{E}}}{g_{A}^{2}}M_{F,sd}\right),\\\
\mathcal{M}_{E,R}^{(6)}&=-\frac{V_{ud}C_{VR}^{(6)}}{3}\left(\frac{g_{V}^{2}}{g_{A}^{2}}M_{F}-\frac{1}{3}\left(2M_{GT}^{AA}+M_{T}^{AA}\right)+\frac{6\mathbf{g_{V\;R}^{E}}}{g_{A}^{2}}M_{F,sd}\right),\end{split}$
(10)

$\mathcal{A}_{m_{e}}$ is determined by

$\displaystyle\begin{split}\mathcal{M}_{m_{e},L}^{(6)}&=\frac{V_{ud}C_{VL}^{(6)}}{6}\bigg{(}\frac{g_{V}^{2}}{g_{A}^{2}}M_{F}-\frac{1}{3}\left(M_{GT}^{AA}-4M_{T}^{AA}\right)-3\left(M_{GT}^{AP}+M_{GT}^{PP}+M_{T}^{AP}+M_{T}^{PP}\right)\\\
&\qquad\qquad\;\;\;\;\;-\frac{12\mathbf{g_{V\;L}^{m_{e}}}}{g_{A}^{2}}M_{F,sd}\bigg{)},\\\
\mathcal{M}_{m_{e},R}^{(6)}&=\frac{V_{ud}C_{VR}^{(6)}}{6}\bigg{(}\frac{g_{V}^{2}}{g_{A}^{2}}M_{F}+\frac{1}{3}\left(M_{GT}^{AA}-4M_{T}^{AA}\right)+3\left(M_{GT}^{AP}+M_{GT}^{PP}+M_{T}^{AP}+M_{T}^{PP}\right)\\\
&\qquad\qquad\;\;\;\;\;-\frac{12\mathbf{g_{V\;R}^{m_{e}}}}{g_{A}^{2}}M_{F,sd}\bigg{)},\end{split}$
(11)

and finally $\mathcal{A}_{M}$ is given by

$\displaystyle\begin{split}\mathcal{M}_{M}^{(6)}&=V_{ud}C_{VL}^{(6)}\bigg{[}2\frac{g_{A}}{g_{M}}\left(M_{GT}^{MM}+M_{T}^{MM}\right)\\\
&\qquad\qquad\;\;\;\;\;+\frac{m_{\pi}^{2}}{m_{N}^{2}}\left(-\frac{2}{g_{A}^{2}}\mathbf{g_{VL}^{NN}}M_{F,sd}+\frac{1}{2}\mathbf{g_{VL}^{\pi
N}}\left(M_{GT,sd}^{AP}+M_{T,sd}^{AP}\right)\right)\bigg{]},\\\
\mathcal{M}_{M}^{(9)}&=\frac{m_{\pi}^{2}}{m_{N}^{2}}\bigg{[}-\frac{2}{g_{A}^{2}}\big{(}\mathbf{g_{6}^{NN}}C_{V}^{(9)}+\mathbf{g_{7}^{NN}}\tilde{C}_{V}^{(9)}\big{)}M_{F,sd}\\\
&\qquad\qquad+\frac{1}{2}\big{(}\mathbf{g_{V}^{\pi
N)}}C_{V}^{(9)}+\mathbf{\tilde{g}_{V}^{\pi
N}}\tilde{C}_{V}^{(9)}\big{)}\big{(}M_{GT,sd}^{AP}+M_{T,sd}^{AP}\big{)}\bigg{]}.\end{split}$
(12)

In the above formulas we have defined the combined NMEs

$\displaystyle\begin{split}\mathcal{M}_{GT}&=M_{GT}^{AA}+M_{GT}^{AP}+M_{GT}^{PP}+M_{GT}^{MM},\\\
\mathcal{M}_{T}&=M_{T}^{AP}+M_{T}^{PP}+M_{T}^{MM},\\\
\mathcal{M}_{PS}&=\frac{1}{2}M_{GT}^{AP}+M_{GT}^{PP}+\frac{1}{2}M_{T}^{AP}+M_{T}^{PP},\\\
\mathcal{M}_{T6}&=2\frac{\mathbf{g_{T}^{\prime}}-\mathbf{g_{T}^{NN}}}{g_{A}^{2}}\frac{m_{\pi}^{2}}{m_{N}^{2}}M_{F,sd}-\frac{8g_{T}}{g_{M}}\left(M_{GT}^{MM}+M_{T}^{MM}\right)\\\
&\quad+\mathbf{g_{T}^{\pi
N}}\frac{m_{\pi}^{2}}{4m_{N}^{2}}\left(M_{GT,sd}^{AP}+M_{T,sd}^{AP}\right)+\mathbf{g_{T}^{\pi\pi}}\frac{m_{\pi}^{2}}{4m_{N}^{2}}\left(M_{GT,sd}^{PP}+M_{T,sd}^{PP}\right).\end{split}$
(13)

The short-range dimension-9 LEFT operators contribute to the $C_{V,\pi\pi
L,\pi NL,NNL}^{(9)}$ couplings that appear in the chiral Lagrangian. They are
given by

$\displaystyle\begin{split}&C_{V}^{(9)}=C_{6}^{(9)}+C_{6}^{(9)}{{}^{\prime}}+C_{8}^{(9)}+C_{8}^{(9)}{{}^{\prime}}\,,\quad\tilde{C}_{V}^{(9)}=C_{7}^{(9)}+C_{7}^{(9)}{{}^{\prime}}+C_{9}^{(9)}+C_{9}^{(9)}{{}^{\prime}}\\\
\hfil\displaystyle\begin{split}&C_{\pi\pi
L}^{(9)}=g_{2}^{\pi\pi}\left(C_{2L}^{(9)}+C_{2L}^{(9)}{{}^{\prime}}\right)+g_{3}^{\pi\pi}\left(C_{3L}^{(9)}+C_{3L}^{(9)}{{}^{\prime}}\right)-g_{4}^{\pi\pi}C_{4L}^{(9)}-g_{5}^{\pi\pi}C_{5L}^{(9)}-\frac{5}{3}g_{1}^{\pi\pi}m_{\pi}^{2}\left(C_{1L}^{(9)}+C_{1L}^{(9)}{{}^{\prime}}\right)\end{split}\\\
&C_{\pi NL}^{(9)}=\left(\mathbf{g_{1}^{\pi
N}}-\frac{5}{6}g_{1}^{\pi\pi}\right)\left(C_{1L}^{(9)}+C_{1L}^{(9)}{{}^{\prime}}\right)\\\
\hfil\displaystyle\begin{split}&C_{NNL}^{(9)}=\mathbf{g_{1}^{NN}}\left(C_{1L}^{(9)}+C_{1L}^{(9)}{{}^{\prime}}\right)+\mathbf{g_{2}^{NN}}\left(C_{2L}^{(9)}+C_{2L}^{(9)}{{}^{\prime}}\right)+\mathbf{g_{3}^{NN}}\left(C_{3L}^{(9)}+C_{3L}^{(9)}{{}^{\prime}}\right)+\mathbf{g_{4}^{NN}}C_{4L}^{(9)}+\mathbf{g_{5}^{NN}}C_{5L}^{(9)}\end{split}\\\
&C_{\left\\{\pi\pi,\pi N,NN\right\\}R}=C_{\left\\{\pi\pi,\pi
N,NN\right\\}L}|_{L\rightarrow R}\;.\end{split}$ (14)

The two LECs $g_{V}^{\pi N}$ and $\tilde{g}_{V}^{\pi N}$ are defined as

$\displaystyle\begin{split}\mathbf{g_{V}^{\pi N}}=\mathbf{g_{6}^{\pi
N}}+\mathbf{g_{8}^{\pi N}}\\\ \mathbf{\tilde{g}_{V}^{\pi
N}}=\mathbf{g_{7}^{\pi N}}+\mathbf{g_{9}^{\pi N}}\end{split}$ (15)

For the sake of convenience, we have marked all currently unknown LECs
including $g_{\nu}^{NN}$ in bold within the above formulas.

In this work we will study the different LEFT operators at the matching scale
of $\Lambda=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 $\chi$PT at
$\Lambda_{\chi}\simeq 2\,$GeV is described in [7].

###  2.3 Relation to literature

A different basis to describe $0\nu\beta\beta$ 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

$\displaystyle\mathcal{L}_{6}=\frac{G_{F}}{\sqrt{2}}\sum_{i,k}\epsilon_{k}^{i}j_{i}J_{k}$
(16)

for the long-range part with $i,k\in\left\\{V\pm A,S\pm
P,T_{L},T_{R},\right\\}$ and

$\displaystyle\begin{split}\mathcal{L}_{9}=\frac{G_{F}^{2}}{2m_{N}}\sum_{l,m,n}\Big{[}&\epsilon_{1}^{lmn}J_{l}J_{m}j_{n}+\epsilon_{2}^{lmn}J^{\mu\nu}_{l}{J_{\mu\nu}}_{m}j_{n}+\epsilon_{3}^{lmn}J^{\mu}_{l}{J_{\mu}}_{m}j_{n}\\\
&+\epsilon_{4}^{lmn}J^{\mu}_{l}{J_{\mu\nu}}_{m}j^{\nu}+\epsilon_{5}^{lmn}J^{\mu}_{l}J_{m}j_{\mu
n}\Big{]}\end{split}$ (17)

for the short-range part with $l,m,n\in\left\\{L,R\right\\}$. Here,
$\epsilon_{k}^{i}$ and $\epsilon^{lmn}$ are the Wilson coefficients of the
different long- and short-range operators. The quark currents $J$ are given
by111We keep the two different types of indices for the short-range currents
to stick with the literature

$\displaystyle\begin{split}J_{S\pm
P}=J_{R,L}&=\overline{u}\left(1\pm\gamma_{5}\right)d\;,\quad J_{V\pm
A}=J_{R,L}^{\mu}=\overline{u}\gamma^{\mu}\left(1\pm\gamma_{5}\right)d,\\\
J_{T_{R,L}}=J^{\mu\nu}_{R,L}&=\overline{u}\sigma^{\mu\nu}\left(1\pm\gamma_{5}\right)d,\end{split}$
(18)

and the lepton currents $j$ are given by

$\displaystyle\begin{split}j_{S\pm
P}&=\overline{e}\left(1\pm\gamma_{5}\right)\nu^{c}\;,\quad j_{V\pm
A}=\overline{e}\gamma^{\mu}\left(1\pm\gamma_{5}\right)\nu^{c}\\\
j_{T_{R,L}}&=\overline{e}\sigma^{\mu\nu}\left(1\pm\gamma_{5}\right)\nu^{c}\,,\quad
j_{R,L}=\overline{e}\left(1\pm\gamma_{5}\right)e^{c}\\\
j_{R,L}^{\mu}&=\overline{e}\gamma^{\mu}\left(1\pm\gamma_{5}\right)e^{c}.\end{split}$
(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

$\displaystyle C_{1L}^{(9)}=\frac{2v}{m_{N}}\epsilon^{LLL}_{3},\quad
C_{1L}^{(9)}{{}^{\prime}}=\frac{2v}{m_{N}}\epsilon^{RRL}_{3},\quad
C_{1R}^{(9)}=\frac{2v}{m_{N}}\epsilon^{LLR}_{3},\quad
C_{1R}^{(9)}{{}^{\prime}}=\frac{2v}{m_{N}}\epsilon^{RRR}_{3},$ (20)
$\displaystyle\begin{split}&C_{2L}^{(9)}=\frac{2v}{m_{N}}\left(\epsilon_{1}^{LLL}-4\epsilon^{LLL}_{2}\right),\quad
C_{2L}^{(9)}{{}^{\prime}}=\frac{2v}{m_{N}}\left(\epsilon_{1}^{RRL}-4\epsilon^{RRL}_{2}\right),\\\
&C_{2R}^{(9)}=\frac{2v}{m_{N}}\left(\epsilon_{1}^{LLR}-4\epsilon^{LLR}_{2}\right),\quad
C_{2R}^{(9)}{{}^{\prime}}=\frac{2v}{m_{N}}\left(\epsilon_{1}^{RRR}-4\epsilon^{RRR}_{2}\right),\end{split}$
(21)
$\displaystyle\begin{split}&C_{3L}^{(9)}=-\frac{16v}{m_{N}}\epsilon_{2}^{LLL},\quad
C_{3L}^{(9)}{{}^{\prime}}=-\frac{16v}{m_{N}}\epsilon_{2}^{RRL},\\\
&C_{3R}^{(9)}=-\frac{16v}{m_{N}}\epsilon_{2}^{LLR},\quad
C_{3R}^{(9)}{{}^{\prime}}=-\frac{16v}{m_{N}}\epsilon_{2}^{RRR},\end{split}$
(22) $\displaystyle C_{4L}^{(9)}=\frac{2v}{m_{N}}\epsilon_{3}^{RLL},\quad
C_{4R}^{(9)}=\frac{2v}{m_{N}}\epsilon_{3}^{RLR},$ (23) $\displaystyle
C_{5L}^{(9)}=-\frac{v}{m_{N}}\epsilon_{1}^{RLL},\quad
C_{5R}^{(9)}=-\frac{v}{m_{N}}\epsilon_{1}^{RLR},$ (24) $\displaystyle
C_{6}^{(9)}=\frac{v}{m_{N}}\left(\epsilon_{5}^{LRR}+i\frac{5}{3}\epsilon_{4}^{LRR}\right),\quad
C_{6}^{(9)}{{}^{\prime}}=\frac{v}{m_{N}}\left(\epsilon_{5}^{RLR}+i\frac{5}{3}\epsilon_{4}^{RLR}\right),$
(25) $\displaystyle C_{7}^{(9)}=4i\frac{v}{m_{N}}\epsilon_{4}^{LRR},\quad
C_{7}^{(9)}{{}^{\prime}}=4i\frac{v}{m_{N}}\epsilon_{4}^{RLR},$ (26)
$\displaystyle
C_{8}^{(9)}=\frac{v}{m_{N}}\left(\epsilon_{5}^{LLR}-i\frac{5}{3}\epsilon_{4}^{LLR}\right),\quad
C_{8}^{(9)}{{}^{\prime}}=\frac{v}{m_{N}}\left(\epsilon_{5}^{RRR}-i\frac{5}{3}\epsilon_{4}^{RRR}\right),$
(27) $\displaystyle C_{9}^{(9)}=-4i\frac{v}{m_{N}}\epsilon_{4}^{LLR},\quad
C_{9}^{(9)}{{}^{\prime}}=-4i\frac{v}{m_{N}}\epsilon_{4}^{RRR}.$ (28)

The main difference between the two set of operators is that the
$\epsilon$-basis contains short-range tensor operators instead of color
octets.

##  3 Distinguishing the Effective Operators

Neutrinoless double-$\beta$-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\nu\beta\beta$
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 $\beta\beta$ 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\nu\beta\beta$ 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\nu\beta\beta$ decay alongside the standard mass mechanism.

###  3.1 Phase-Space Observables

While most experimental collaborations only attempt to measure the half-life
of $0\nu\beta\beta$ 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]

Figure 1: Comparison of the normalized single electron spectra (lower left)
and angular correlation coefficients (upper right) in 136Xe that result from
the 6 PSFs which appear in the $0\nu\beta\beta$ half-life “master-formula”.
Red curves correspond to the red-labelled PSFs on horizontal axis, while blue
curves represent the PSFs denoted in blue on the vertical axis. The x-axis
covers the range $\tilde{\epsilon}\in\left[0-1\right]$.

$\displaystyle\begin{split}G_{0k}=\frac{(G_{F}V_{ud})^{4}m_{e}^{2}}{64\pi^{5}\ln{2}R^{2}}&\int\delta\bigg{(}\epsilon_{1}+\epsilon_{2}+E_{f}-E_{i}\bigg{)}\\\
&\times\bigg{(}h_{0k}(\epsilon_{1},\epsilon_{2},R)\cos{\theta}+g_{0k}(\epsilon_{1},\epsilon_{2},R)\bigg{)}\\\
&\times
p_{1}p_{2}\epsilon_{1}\epsilon_{2}\,\text{d}\epsilon_{1}\,\text{d}\epsilon_{2}\,\text{d}(\cos{\theta}),\end{split}$
(29)

where $p_{1,2}$ and $\epsilon_{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 final-state 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 $\theta$. Additionally, $G_{04,06,09}$ have to be rescaled to
comply with the definitions in [6, 7] as

$\displaystyle\begin{split}&G_{04}\rightarrow\frac{9}{2}G_{04},\\\
&G_{06}\rightarrow\frac{m_{e}R}{2}G_{06},\\\
&G_{09}\rightarrow\left(\frac{m_{e}R}{2}\right)^{2}G_{09}.\end{split}$ (30)

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. (29) one can write the angular correlation coefficient
$a_{1}/a_{0}$ which is defined via

$\displaystyle\frac{\text{d}\Gamma}{\mathrm{d}\cos{\theta}\mathrm{d}\tilde{\epsilon}_{1}}=a_{0}\left(1+\frac{a_{1}}{a_{0}}\cos{\theta}\right)$
(31)

with

$\displaystyle\tilde{\epsilon}_{i}=\frac{\epsilon_{i}-m_{e}}{Q_{\beta\beta}}\in\left[0,1\right]$
(32)

as

$\displaystyle\frac{a_{1}}{a_{0}}(\tilde{\epsilon})=\frac{\sum_{i}|M_{i}|^{2}h_{0i}\left(\epsilon,\Delta
M_{\text{Nuclei}}-\epsilon,R\right)}{\sum_{j}|M_{j}|^{2}g_{0j}\left(\epsilon,\Delta
M_{\text{Nuclei}}-\epsilon,R\right)}\,.$ (33)

Here, $\Delta M_{\text{Nuclei}}$ is the mass difference between the mother and
daughter nuclei and $Q_{\beta\beta}$ 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\nu\beta\beta$ has been
discussed e.g. in [55]. Similarly, the single electron spectra are given by

$\displaystyle\begin{split}\frac{\text{d}\Gamma}{\text{d}\epsilon_{1}}=&\frac{(G_{F}V_{ud})^{4}m_{e}^{2}}{64\pi^{5}\ln{2}R^{2}}\bigg{(}\sum_{i}|M_{i}|^{2}g_{0i}\left(\epsilon,\Delta
M_{\text{Nuclei}}-\epsilon,R\right)\bigg{)}p_{1}p_{2}\epsilon\big{(}\Delta
M_{\text{Nuclei}}-\epsilon\big{)}\,.\end{split}$ (34)

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.

Figure 2: Normalized single electron spectra (lower left) and angular
correlation coefficients (upper right) for each of the 4 distinguishable
groups of operators. The shapes are shown for 136Xe assuming the NDA values
for the currently unknown LECs. However, the particular choice does not result
in a significant difference in the general shape of the plots. Red curves
correspond to the red-labelled operator group on the horizontal axis, while
blue curves represent the operator group denoted in blue on the vertical axis.
The x-axis covers the range $\tilde{\epsilon}\in\left[0-1\right]$.

However, distinguishing among different $0\nu\beta\beta$ mechanisms purely
based on the phase-space 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_{VL}^{(6)},C_{VR}^{(6)}$, 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_{VL}^{(6)}$ and the operators corresponding to
$G_{09}$, while corresponding to distinct PSFs, are practically
indistinguishable since the $C_{VL}^{(6)}$ 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\nu\beta^{-}\beta^{-}$ isotopes.

Figure 3:  The single electron spectra for 4 different naturally occurring
$0\nu\beta^{-}\beta^{-}$ isotopes are shown. While the exact quantitative
curves depend on the choice of the isotope, their shape is mostly independent
of this choice. As before, the x-axis shows the normalized electron energy
$\tilde{\epsilon}$.

###  3.2 Decay Rate Ratios

The remaining $0\nu\beta\beta$ observable is the decay rate $\Gamma$ 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

$\displaystyle
R^{\mathcal{O}_{i}}(^{\text{A}}\text{X})\equiv\frac{T_{1/2}^{\mathcal{O}_{i}}(^{\text{A}}\text{X})}{T_{1/2}^{\mathcal{O}_{i}}(^{76}\text{Ge})}=\frac{\sum_{j}|\mathcal{M}_{j}^{\mathcal{O}_{i}}(^{76}\text{Ge})|^{2}G_{j}^{\mathcal{O}_{i}}(^{76}\text{Ge})}{\sum_{k}|\mathcal{M}_{k}^{\mathcal{O}_{i}}(^{\text{A}}\text{X})|^{2}G_{k}^{\mathcal{O}_{i}}(^{\text{A}}\text{X})}$
(35)

where $T_{1/2}^{\mathcal{O}_{i}}(^{A}X)$ is the half-life induced by the
operator $O_{i}$ in the isotope ${}^{A}X$. The sums $\sum_{j,k}$ are taken
over all different PSFs generated by the operator $\mathcal{O}_{i}$ and become
relevant only for $C_{VL,VR}^{(6)}$ (see Eqs. (89) and (90)). Studying the
half-life ratio allows for elimination of the unknown particle physics
couplings, as was first discussed in [56] and shortly after also in [57].
Here, we take 76Ge for the reference isotope. To be able to quantify how well
one can distinguish two different operators $\mathcal{O}_{i,j}$ from each
other we can take the ratio

$\displaystyle
R_{ij}(^{A}X)=\frac{R^{\mathcal{O}_{i}}(^{A}X)}{R^{\mathcal{O}_{j}}(^{A}X)}.$
(36)

Specifically, the ratios $R_{im_{\beta\beta}}$ 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\nu\beta\beta$ rate in
experiments. Obviously, two operators $\mathcal{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}\rightarrow\infty$ or $R_{ij}=0$, that is, for
$|\log_{10}(R_{ij})|\rightarrow\infty.$

Figure 4: Comparison of the angular correlation coefficients in 4 different
isotopes as done for the single electron energy spectra in Fig. 3.

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\nu\beta\beta$ decay have also
been employed previously in [59, 60, 61, 62].

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_{6,7}^{NN}=g_{V}^{\pi N}=\tilde{g}_{V}^{\pi N}=1$ (so that the contribution
from the short-range vector operators is not omitted) makes it impossible to
distinguish the short-range scalar operators $C_{S2-S5}^{(9)}$ as well as the
short-range vector operator groups $C^{(9)}_{V}$ and $\tilde{C}^{(9)}_{V}$.

$m_{\beta\beta}$ | $C_{VL}^{(6)}$ | $C_{VR}^{(6)}$ | $C_{T}^{(6)}$ | $C_{S,V}^{(6,7)}$ | $C_{S1}^{(9)}$ | $C_{S2}^{(9)}$ | $C_{S3}^{(9)}$ | $C_{S4}^{(9)}$ | $C_{S5}^{(9)}$ | $C_{V}^{(9)}$ | $\tilde{C}_{V}^{(9)}$
---|---|---|---|---|---|---|---|---|---|---|---
$m_{\beta\beta}$ | $C_{VL}^{(6)}$ | $C_{VR}^{(6)}$ | $C_{T}^{(6)}$ | $C_{SL}^{(6)}$ | $C_{1L}^{(9)}$ | $C_{2L}^{(9)}$ | $C_{3L}^{(9)}$ | $C_{4L}^{(9)}$ | $C_{5L}^{(9)}$ | $C_{6}^{(9)}$ | $C_{7}^{(9)}$
- | - | - | - | $C_{SR}^{(6)}$ | $C_{1R}^{(9)}$ | $C_{2R}^{(9)}$ | $C_{3R}^{(9)}$ | $C_{4R}^{(9)}$ | $C_{5R}^{(9)}$ | $C_{6}^{(9)}{{}^{\prime}}$ | $C_{7}^{(9)}{{}^{\prime}}$
- | - | - | - | $C_{VL}^{(7)}$ | $C_{1L}^{(9)}{{}^{\prime}}$ | $C_{2L}^{(9)}{{}^{\prime}}$ | $C_{3L}^{(9)}{{}^{\prime}}$ | - | - | $C_{8}^{(9)}$ | $C_{9}^{(9)}$
- | - | - | - | $C_{VR}^{(7)}$ | $C_{1R}^{(9)}{{}^{\prime}}$ | $C_{2R}^{(9)}{{}^{\prime}}$ | $C_{3R}^{(9)}{{}^{\prime}}$ | - | - | $C_{8}^{(9)}{{}^{\prime}}$ | $C_{9}^{(9)}{{}^{\prime}}$
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}^{(9)}$ become indistinguishable as well as the short-range vector
operator groups $\tilde{C}_{V}^{(9)}$ and $\tilde{C}_{V}^{(9)}$. Improved
knowledge of the LECs, assuming no fine tuning, would allow to distinguish
among these operator groups.

####  3.2.1 Sensitivity on the unknown LECs

In Figure 5 we present the expected ratios $R^{\mathcal{O}_{i}}$ as well as
the normalized ratios $R_{im_{\beta\beta}}$ defined in (35) and (36) for the
above choice of LECs which will be our benchmark scenario. The ratios for the
$\epsilon$-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 $\left[-\sqrt{10},-1/\sqrt{10}\right]\times|g_{i}|$ and
$\left[1/\sqrt{10},\sqrt{10}\right]\times|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_{\nu}^{NN}$ which generates a short-range component into the
standard mass-mechanism we take a variation of $\pm 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_{\beta\beta}}$ corresponding to the standard mass mechanism are not very
sensitive to $g_{\nu}^{NN}$ (they are actually too small to be visible). In
Figure 7 we explicitly show the impact of varying the $g_{\nu}^{NN}$ LEC on
the expected half-life in 76Ge for the standard mechanism. Again, compared to
the impact of the unknown Majorana phases the effect of $g_{\nu}^{NN}$ is
minor. However, it is important to note that the impact of $g_{\nu}^{NN}$ on
the overall magnitude of the half-life cannot be ignored as easily. For
comparison, we also present the case where $g_{\nu}^{NN}=0$ in Figure 7.

Figure 5: The decay rate ratios $R^{O_{i}}$ (upper plot) and
$R_{im_{\beta\beta}}$ (lower plot) for the different operator groups are
shown. The larger markers represent the choice of vanishing unknown LECs with
$g_{6,7}^{NN}=g_{V}^{\pi N}=\tilde{g}_{V}^{\pi N}=1$. Isotopes with a PSF
$G_{0}>10^{-14}\,\mathrm{y^{-1}}$ are represented by stars while isotopes with
smaller PSFs are represented by round markers. The additional points represent
variations of the different unknown LECs $g_{i}$ randomly chosen from
$\left[-\sqrt{10},-1/\sqrt{10}\right]\times|g_{i}|$ and
$\left[+1/\sqrt{10},+\sqrt{10}\right]\times|g_{i}|$ except for $g_{\nu}^{NN}$
which is varied in a range of $\pm 50\%$. The crosses represent the central
values of the variation i.e. the median values. The reference isotope is
chosen to be 76Ge. Note that the variation of $g_{\nu}^{NN}$ does induce a
small variation of $R^{m_{\beta\beta}}$ which is, however, not visible in the
above plot.

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_{i}^{(9)}$ 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\nu\beta\beta$-decay rate. Although for some operator groups,
such as $C_{2S-5S}^{(9)}$, 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_{V}^{(9)},\tilde{C}_{V}^{(9)}$ 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}^{\pi\pi}$, 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_{\nu}^{NN}$) 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_{\beta\beta}}$ 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].

####  3.2.2 Distinguishing different operators

Figure 6: The decay rate ratios $R^{O_{i}}$ (upper plot) and
$R_{im_{\beta\beta}}$ (lower plot) for the different operator groups in the
$\epsilon$-basis similar to Figure 5 are shown.

When studying the $0\nu\beta\beta$ 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\nu\beta\beta$ decay
contribution triggered by a particular higher-dimensional operator. In the
first row of Figure 8 we show the maximal ratios
$R_{im_{\beta\beta}}^{\mathrm{max}}$ 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_{\beta\beta}}$ 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: 76Ge [21], 82Se
[22], 96Zr [23], 100Mo [24], 116Cd [25], 128Te [26], 130Te [27], 134Xe [66],
136Xe [29] and 150Nd [30]. Figure 8 also presents all the other ratios
$R_{ij}^{\mathrm{max}}$ 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.

Figure 7: Here we show the half-life for 76Ge in dependence on the minimal
neutrino mass $m_{\mathrm{min}}$ for both normal (NO) and inverted (IO)
neutrino mass ordering. The scatter points were obtained by marginalizing over
$g_{\nu}^{NN}$ and the unknown Majorana phases. The blue and red contours show
the possible half-life ranges when only the phases are varied, while
$g_{\nu}^{NN}$ is fixed. Additionally, the black contours correspond to
$g_{\nu}^{NN}=0$.

Considering the central values, the non-standard long-range operators
$C_{VR}^{(6)}$ and $C_{S,V}^{(6,7)}$ give the most distinct half-life ratios
compared to the standard mass mechanism while the remaining long-range
operator groups $C_{VL}^{(6)}$ and $C_{T}^{(6)}$ also result in sizable
$R_{im_{\beta\beta}}^{\mathrm{max}}>2$. Additionally, both $C_{VL}^{(6)}$ and
$C_{VR}^{(6)}$ 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_{im_{\beta\beta}}^{\mathrm{max}}<2$; thus,
there is less potential for identifying their contribution by experiments.
However, the short-range vector operators $C_{V}^{(9)}$ and
$\tilde{C}_{V}^{(9)}$ are associated with a different angular correlation than
$m_{\beta\beta}$. On the contrary, the contributions from scalar short-range
operators $C_{S2-S5}^{(9)}$ would be hardest to discriminate, as they do not
manifest any significant isotope dependence on $R_{im_{\beta\beta}}$ and do
not differ in the phase-space observables, either.

In the worst-case scenario, the operators in the group $C_{S,V}^{(6,7)}$,
i.e., lepton number violating long-range scalar and vector interactions, are
the only operators that result in $R_{im_{\beta\beta}}^{\mathrm{max}}>2$. In
fact, this is the only operator group that is not affected by any unknown
LECs. Besides $C_{T}^{(6)}$ all remaining operators in this setting have
expected ratios $R_{im_{\beta\beta}}^{\mathrm{max}}\leq 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 contributions from $C_{VL}^{(6)},C_{V}^{(9)}$ and
$\tilde{C}_{V}^{(9)}$ could be identified only based on measurements of the
angular correlation and the scalar short-range operators, $C_{S1}^{(9)}$,
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.

Figure 8: The maximal ratios $R_{ij}^{\mathrm{max}}$ for all operator
combinations $i,j$. The exact values and the corresponding isotopes are
displayed in each tile. Additionally, operator combinations that result in
different phase-space observables are marked by dashed-line shading. In the
upper right half of the plot we show the ratios considering the central values
from the variation of the LECs. In the lower left half we show the worst-case
scenario considering the values of ratios $R_{ij}$ that are closest to 1
within the range obtained by the variation of the LECs.

To be able to pinpoint the specific non-standard operator group
$\mathcal{O}_{j}$ contributing to $0\nu\beta\beta$ 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 right-
handed vector current $C_{VR}^{(6)}$ for which all the ratios
$R_{iC_{VR}^{(6)}}^{\mathrm{max}}$ are large i.e. $\gtrsim 7.6$.

####  3.2.3 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\nu\beta\beta$ 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 _effective_ NME which we label
$M_{\mathrm{eff}}$,

$\displaystyle
T_{1/2}^{-1}=\left|C\right|^{2}G\left|M_{\mathrm{eff}}\right|^{2}.$ (37)

Here, $M_{\mathrm{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).

Figure 9: Comparison of the standard NME
$\mathcal{M}_{0\nu}=-\frac{1}{g_{A}^{2}}M_{F}+M_{GT}+M_{T}$ resulting from
different calculation methods. Explicitly, we show NMEs obtained from the
interacting shell model (SM) [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 $\chi$PT [76, 77]. The grey bands
mark the range of values covered by the different methods.

If we consider the theoretical uncertainty of the half-life to be dominated by
the uncertainty of $M_{\mathrm{eff}}$, we can determine the necessary
theoretical accuracy of the nuclear physics. To estimate this, we assume
$M_{\mathrm{eff}}$ to be independent of the choice of the isotope, i.e.,

$\displaystyle\frac{\Delta
M_{\mathrm{eff}}}{M_{\mathrm{eff}}}(^{A}Z)=\frac{\Delta
M_{\mathrm{eff}}}{M_{\mathrm{eff}}}=\mathrm{const.}.$ (38)

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,

$\displaystyle\Delta R_{ij}\overset{!}{<}|R_{ij}-1|\;.$ (39)

Hence, the necessary theoretical accuracy for $M_{\mathrm{eff}}$ reads

$\displaystyle\frac{\Delta
M_{\mathrm{eff}}}{M_{\mathrm{eff}}}\overset{!}{<}\frac{1}{4}\frac{|R_{ij}-1|}{R_{ij}}.$
(40)

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 $\frac{\Delta M_{\mathrm{eff}}}{M_{\mathrm{eff}}}\sim 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_{VR}^{(6)}$
requiring an accuracy of $\frac{\Delta M_{\mathrm{eff}}}{M_{\mathrm{eff}}}\sim
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 $\sim 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\nu\beta\beta$ decay NMEs seem to pave the path
towards more reliable numerical values and clearer understanding of the
theoretical uncertainties involved.

###  3.3 Other $\mathbf{0\nu\beta\beta}$ Modes

Besides the usual $0\nu\beta^{-}\beta^{-}$-decay mode one could also make use
of searches for neutrinoless modes of other $\beta\beta$ processes. In general
there are four of these,

  1. 1.

$\mathbf{\beta^{-}\beta^{-}}$:

$\displaystyle(A,Z)$
$\displaystyle\rightarrow(A,Z+2)+2e^{-}(+2\overline{\nu}_{e})$ (41)
$\displaystyle\Delta m$ $\displaystyle\overset{!}{>}0$ (42)

  2. 2.

$\mathbf{\beta^{+}\beta^{+}}$:

$\displaystyle(A,Z)$ $\displaystyle\rightarrow(A,Z-2)+2e^{+}(+2\nu_{e})$ (43)
$\displaystyle\Delta m$ $\displaystyle\overset{!}{>}4m_{e}$ (44)

  3. 3.

$\text{EC}\mathbf{\beta^{+}}$:

$\displaystyle(A,Z)+e^{-}$ $\displaystyle\rightarrow(A,Z-2)+e^{+}(+2\nu_{e})$
(45) $\displaystyle\Delta m$ $\displaystyle\overset{!}{>}2m_{e}$ (46)

  4. 4.

ECEC:

$\displaystyle(A,Z)+2e^{-}$ $\displaystyle\rightarrow(A,Z-2)(+2\nu_{e})$ (47)
$\displaystyle\Delta m$ $\displaystyle\overset{!}{>}0\qquad\text{2$\nu$ mode}$
(48) $\displaystyle\Delta m$ $\displaystyle\overset{!}{=}0\qquad\text{0$\nu$
mode}.$ (49)

Here, $\mathrm{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 $\beta\beta$ processes listed in Table 6
are expected to have half-lives significantly longer than the usual
$0\nu\beta^{-}\beta^{-}$ decay and are therefore unlikely to show up in
experiments. Despite that let us discuss their potential role in a bit more
detail.

####  3.3.1 $\mathbf{0\nu\beta^{+}\beta^{+}}$

This process can be treated in a similar way as the usual
$0\nu\beta^{-}\beta^{-}$ decay, one only needs to consider a negative nuclear
charge $Z\rightarrow-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\nu\beta^{-}\beta^{-}$ $Q$-values. Additionally, the
electromagnetic repulsion of the outgoing positrons deforms the wave functions
and decreases the decay rate. Thus, we see that $0\nu\beta^{+}\beta^{+}$ will
be highly suppressed compared to $0\nu\beta^{-}\beta^{-}$. Also, given the
similarities of the two decays there does not seem to be a natural way of
enhancing the $0\nu\beta^{+}\beta^{+}$-decay rate with respect to
$0\nu\beta^{-}\beta^{-}$. The relevant PSFs for the
$0\nu\beta^{+}\beta^{+}$-decay of naturally occurring isotopes have been
calculated in [78] and are about 3-5 orders of magnitude smaller than for
$0\nu\beta^{-}\beta^{-}$ decay.

####  3.3.2 $\mathbf{0\nu EC\beta^{+}}$

The PSFs for the neutrinoless mode of $\text{EC}\beta^{+}$ were also
calculated with good precision in [78] and are found to be 3-4 orders of
magnitude smaller than those corresponding to $0\nu\beta^{-}\beta^{-}$ decay.
The reason is the same as in case of $0\nu\beta^{+}\beta^{+}$ decay.

####  3.3.3 $\mathbf{0\nu 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\nu\beta^{-}\beta^{-}$ [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\nu\beta^{-}\beta^{-}$ decay. Recently, it was shown that further
significant enhancement of the ${0\nu\mathrm{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\nu\mathrm{ECEC}$ decay can be found in [82].

####  3.3.4 Bound-State $\mathbf{0\nu\beta\beta}$

Bound state $0\nu\beta\beta$ 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\nu\beta\text{EP}$ and
$0\nu\text{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\nu\text{ECEC}$, also $0\nu\text{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-$\beta$-decay $0\nu\beta\text{EP}$ was
investigated in [83] and found to have PSFs 6-7 orders of magnitude smaller
than those of $0\nu\beta^{-}\beta^{-}$ decay. The decay rates can be
significantly enhanced when considering fully ionized nuclei. In that case,
the $0\nu\beta\text{EP}$ decay rate can for certain isotopes even exceed the
one of $0\nu\beta^{-}\beta^{-}$ 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.

####  3.3.5 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.

####  3.3.6 Artificial Isotopes

Although there are 69 naturally occurring double-$\beta$-decaying isotopes, we
found about $\sim 2700$ possible $0\nu\beta\beta$ candidate isotopes when
considering the full NIST list of elements [88]. Some of them have
considerably larger Q-values of up to $50\,\rm{MeV}$.222Considering only
isotopes without a single-$\beta$-decay mode already significantly reduces
this number down to 86. None of these has, however, a significantly enhanced
$Q$-value. While such a large $Q$-value of $\sim 50\,\mathrm{MeV}$ would
result in a significant enhancement of the decay rate by $\sim 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\nu\beta\beta$ decay would be impossible.

To sum up the above paragraphs, despite the fact that a variety of
$\beta\beta$ processes exist, $0\nu\beta^{-}\beta^{-}$ decay is largely the
most relevant candidate to study, indeed. Therefore, other possible
$0\nu\beta\beta$ modes would only become relevant in exotic scenarios leading
either to their significant enhancements, or to strong suppression of
$0\nu\beta^{-}\beta^{-}$ decay. Given the similarity of all the $\beta\beta$
processes, such models would, however, seem to be rather unnatural from a
particle physics point of view.

##  4 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\nu\beta\beta$ decay. Each of the models will be compared with the standard
mass mechanism predictions.

###  4.1 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}\times SU(2)_{L}\times SU(2)_{R}\times U(1)_{B-L}$
[89, 90, 91] which is spontaneously broken to the Standard Model group
$SU(3)_{C}\times SU(2)_{L}\times 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 $\Delta_{L}\in\left(1,3,1,2\right)$ and
$\Delta_{R}\in\left(1,1,3,2\right)$ as well as a scalar bidoublet
$\Phi\in\left(1,2,2^{*},0\right)$ incorporating the SM Higgs doublet and the
right-handed neutrinos $\nu_{R}$. The fermions are grouped into left- and
right-handed doublets

$\displaystyle L_{L}$ $\displaystyle=\left(\begin{matrix}\nu_{L}\\\
e_{L}\end{matrix}\right)\in\left(1,2,1,-1\right),\qquad
Q_{L}=\left(\begin{matrix}u_{L}\\\
d_{L}\end{matrix}\right)\in\left(3,2,1,1/3\right),$ (50) $\displaystyle L_{R}$
$\displaystyle=\left(\begin{matrix}\nu_{R}\\\
e_{R}\end{matrix}\right)\in\left(1,1,2,-1\right),\qquad
Q_{R}=\left(\begin{matrix}u_{R}\\\
d_{R}\end{matrix}\right)\in\left(3,1,2,1/3\right),$ (51)

which under $U\in SU(2)_{L,R}$ transform as

$\displaystyle\Psi_{L,R}\longrightarrow U_{L,R}\Psi_{L,R}$ (52)

while the scalar fields transform as

$\displaystyle\Phi\longrightarrow U_{L}\Phi
U^{\dagger}_{R},\qquad\Delta_{L}\longrightarrow
U_{L}\Delta_{L}U_{L}^{\dagger},\qquad\Delta_{R}\longrightarrow
U_{R}\Delta_{R}U_{R}^{\dagger}.$ (53)

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

$d$$u$$e^{-}$$e^{-}$$d$$u$$W_{L}$$\nu_{i}\,N_{i}$$W_{L}$$d$$u$$e^{-}$$e^{-}$$d$$u$$W_{R}$$\nu_{i}\,N_{i}$$W_{R}$$d$$u$$e^{-}$$e^{-}$$d$$u$$W_{R}$$\nu_{i}\,N_{i}$$W_{L}$$d$$u$$d$$u$$e^{-}$$e^{-}$$W_{L(R)}$$W_{L(R)}$$\Delta_{L(R)}$
Figure 10: Feynman diagrams arising in the mLRSM that contribute to
$0\nu\beta\beta$. Here, $\nu_{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.

$\displaystyle
P:\Psi_{L}\Longleftrightarrow\Psi_{R},\qquad\Phi\Longleftrightarrow\Phi^{\dagger},\qquad\Delta_{L,R}\Longleftrightarrow\Delta_{R,L}$
(54) $\displaystyle
C:\Psi_{L}\Longleftrightarrow\left(\Psi_{R}\right)^{c},\qquad\Phi\Longleftrightarrow\Phi^{T},\qquad\Delta_{L,R}\Longleftrightarrow\Delta_{R,L}^{*}.$
(55)

Requiring either $P$ or $C$ invariance results in different constraints on the
scalar potential as well as the Yukawa coupling matrices [93].333Note that a
combination of both does not fit observational constraints [93].

The lepton number violation at low energy stems from the leptonic Yukawa
interactions given by

$\displaystyle\mathcal{L}_{y}=\sum_{ij}\bigg{\\{}Y^{l}_{ij}{\overline{L_{L}}}_{i}\Phi
L_{R,j}+\tilde{Y}^{l}_{ij}\overline{L_{L}}_{i}\tilde{\Phi}L_{R,j}+Y^{L}_{ij}L^{T}_{L,i}Ci\tau_{2}\Delta_{L}L_{L,j}+{Y^{R}_{ij}}^{\dagger}L^{T}_{R,i}Ci\tau_{2}\Delta_{R}L_{R,j}\bigg{\\}}+\text{h.c.}$
(56)

After the neutral components of the scalars have acquired their VEVs

$\displaystyle\left<\Phi\right>=\frac{1}{\sqrt{2}}\left(\begin{matrix}\kappa&0\\\
0&\kappa^{\prime}e^{i\alpha}\end{matrix}\right),\qquad\left<\Delta_{L}\right>=\frac{1}{\sqrt{2}}\left(\begin{matrix}0&0\\\
v_{L}e^{i\theta_{L}}&0\end{matrix}\right),\qquad\left<\Delta_{R}\right>=\frac{1}{\sqrt{2}}\left(\begin{matrix}0&0\\\
v_{R}&0\end{matrix}\right),$ (57)

one can infer the neutrino mass matrices from (56)

$\displaystyle\begin{split}M^{\nu}_{D,ij}&=\frac{1}{\sqrt{2}}\left[Y^{l}_{ij}\kappa+\tilde{Y}^{l}_{ij}\kappa^{\prime}\exp{-i\alpha}\right],\quad\;\\\
{M^{\nu}_{L,ij}}^{\dagger}&=\sqrt{2}Y^{L}_{ij}v_{L}\exp{i\theta_{L}},\qquad\qquad\quad\;\;\;M^{\nu}_{R,ij}=\sqrt{2}Y^{R}_{ij}v_{R}\;.\end{split}$
(58)

There are several diagrams in the mLRSM setting which can contribute to the
$0\nu\beta\beta$ decay at tree level, see Figure 10. Detailed discussions of
$0\nu\beta\beta$ 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,

$\displaystyle\begin{split}\mathcal{L}_{\Delta
L}=&C^{(5)}\left(\left(L^{T}Ci\tau_{2}\Phi_{SM}\right)\left(\tilde{\Phi}_{SM}^{\dagger}L\right)\right)\\\
&+\left(L^{T}\gamma^{\mu}e_{R}\right)i\tau_{2}\Phi_{SM}\left(C^{(7)}_{Leu\overline{d}\Phi}\overline{d_{R}}\gamma_{\mu}u_{R}+C^{(7)}_{L\Phi
De}\Phi_{SM}^{T}i\tau_{2}(D_{\mu}\Phi_{SM})\right)\\\
&+\overline{e_{R}}e_{R}^{c}\bigg{(}C^{(9)}_{eeud}\overline{u_{R}}\gamma^{\mu}d_{R}\overline{u_{R}}\gamma_{\mu}d_{R}+C^{(9)}_{ee\Phi
ud}\overline{u_{R}}\gamma^{\mu}d_{R}\left(\left[iD_{\mu}\Phi_{SM}\right]^{\dagger}\tilde{\Phi}_{SM}\right)\\\
&\qquad\qquad+.C^{(9)}_{ee\Phi
D}\left(\left[iD_{\mu}\Phi_{SM}\right]^{\dagger}\tilde{\Phi}_{SM}\right)^{2}\bigg{)},\end{split}$
(59)

where $\Phi_{SM}$ is the Standard Model Higgs doublet. The matching scale
corresponds to $\sim m_{W_{R}}$ and the Wilson coefficients at SMEFT level are
given by

$\displaystyle\begin{split}C^{(5)}\;\;\;\;\,&=\frac{1}{v^{2}}\left({M_{D}^{\nu}}^{T}{M_{R}^{\nu}}^{-1}M_{D}^{\nu}-M_{L}^{\nu}\right)\,,\\\
C^{(7)}_{Leu\overline{d}\Phi}&=\frac{\sqrt{2}}{v}\frac{1}{v_{R}^{2}}\left(V_{R}^{ud}\right)^{*}\left(M_{D}^{\nu
T}{M_{R}^{\nu}}^{-1}\right)_{ee}\,,\qquad\qquad\quad C^{(7)}_{L\Phi
De}\,=\frac{2i\xi\exp{i\alpha}}{\left(1+\xi^{2}\right){V_{R}^{ud}}^{*}}C^{(7)}_{Leu\overline{d}\Phi}\,,\\\
C^{(9)}_{eeud}\;\;\,&=-\frac{1}{2v_{R}^{4}}{V_{R}^{ud}}^{2}\left[\left({M_{R}^{\nu}}^{\dagger}\right)^{-1}+\frac{2}{m^{2}_{\Delta_{R}}}M_{R}^{\nu}\right]\,,\qquad
C^{(9)}_{ee\Phi
ud}\,=-4\frac{\xi\exp{-i\alpha}}{\left(1+\xi^{2}\right)V_{R}^{ud}}C^{(9)}_{eeud}\,,\\\
C^{(9)}_{ee\Phi
D}\;&=4\frac{\xi^{2}\exp{-2i\alpha}}{\left(1+\xi^{2}\right)^{2}{V_{R}^{ud}}^{2}}C^{(9)}_{eeud}\,,\end{split}$
(60)

where $v$ is the Standard Model Higgs doublets VEV,

$\displaystyle v^{2}=\kappa^{2}+\kappa^{\prime 2}.$ (61)

Here, $C^{(5)}$ corresponds to the usual seesaw formula. From the matching
scale $\sim m_{W_{R}}$ the above coefficients have to be evolved down to
$m_{W}\sim 80\,\rm{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

$\displaystyle\begin{split}&m_{\beta\beta}=-v^{2}C^{(5)}_{ee}\\\
&C_{VL}^{(6)}=-iV_{L}^{ud}\frac{v^{3}}{\sqrt{2}}{C_{L\Phi
De}^{(7)}}^{*}\,,\qquad\qquad\quad\,C_{VR}^{(6)}=\frac{v^{3}}{\sqrt{2}}{C^{(7)}_{Leu\overline{d}\Phi}}^{*}\,,\\\
&C_{1R}^{(9)}(m_{W})=v^{5}{V_{L}^{ud}}^{2}C_{ee\Phi D}^{(9)}(m_{W})\,,\qquad
C_{1R}^{(9)}{{}^{\prime}}(m_{W})=v^{5}C_{eeud}^{(9)}(m_{W})\,,\\\
&C_{4R}^{(9)}(m_{W})=-v^{5}V_{L}^{ud}C_{ee\Phi ud}^{(9)}(m_{W}).\end{split}$
(62)

Evolving the above coefficients down to the $\chi$PT scale of $\sim
2\,\rm{GeV}$ also generates a non-zero $C_{5R}^{(9)}$ coefficient since the
RGEs of $C_{4,5}^{(9)}$ 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_{\Delta R}$ as well as the masses of the three heavy neutrinos
$\left(m_{\nu_{R1}},m_{\nu_{R2}},m_{\nu_{R3}}\right)$ and the lightest
neutrino mass $m_{\nu_{\text{min}}}$, the complex phases of the VEVs $\alpha$
and $\theta_{L}$ and finally the left-right mixing parameter $\xi$. Here, we
fix $\xi=\frac{m_{b}}{m_{t}}$. The lightest neutrino mass together with the
squared mass differences $\Delta m_{ij}^{2}$ that are known from oscillation
data fix the remaining light neutrino masses for a given mass hierarchy.
Taking

$\displaystyle\ket{\nu_{\alpha\;L,R}}=\sum_{i}U_{\alpha
i}^{*}\ket{\nu_{i\;L,R}}$ (63)

we obtain

$\displaystyle\begin{split}M_{\nu}&=v^{2}C^{(5)}=U_{\text{PMNS}}\,m_{\nu}\,U_{\text{PMNS}}^{T}\;,\qquad\;\;\,m_{\nu}\;=\text{diag}\left(m_{\nu_{1}},m_{\nu_{2}},m_{\nu_{3}}\right)\\\
M^{\nu}_{R}&=U_{R}\,m_{\nu R}\,U_{R}^{T}\;,\qquad\qquad\qquad\qquad\qquad
m_{\nu
R}=\text{diag}\left(m_{{\nu}_{R1}},m_{{\nu}_{R2}},m_{{\nu}_{R3}}\right)\end{split}$
(64)

Additionally, the mixing matrix $U$ for the heavy neutrinos must be fixed.
Here, we take $U_{R}=U_{\text{PMNS}}$ for simplicity. In the $C$-symmetric
case, one has

$\displaystyle M_{L}^{\nu}=\frac{v_{L}\exp{i\theta_{L}}}{v_{R}}M^{R}_{\nu}$
(65)

and the Dirac mass matrix can be derived as [94]

$\displaystyle
M^{\nu}_{D}=U_{\text{PMNS}}\;m_{\nu_{R}}\sqrt{\frac{v_{L}\exp{i\theta_{L}}}{v_{R}}\mathbb{1}_{(3\times
3)}-m_{\nu R}^{-1}m_{\nu}}\;U_{\text{PMNS}}^{T}\;.$ (66)

Figure 11: Half-life ratios resulting from different mLRSM settings (different
neutrino mass hierarchy and different minimal neutrino mass) when taking 76Ge
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.

Assuming $V_{L}^{ud}=V_{R}^{ud}$ and

$\displaystyle\begin{split}&m_{{\nu}_{R1}}=10\,\textrm{TeV}\;,\quad
m_{{\nu}_{R2}}=12\,\textrm{TeV}\;,\quad m_{{\nu}_{R3}}=13\,\textrm{TeV}\\\
&m_{\Delta_{R}}=4\,\textrm{TeV}\,,\quad v_{L}=0.1\,\textrm{eV}\quad
v_{R}=10\,\textrm{TeV}\end{split}$ (67)

as in [7], we can derive the LEFT Wilson coefficients in dependence on the
minimal light neutrino mass $m_{\text{min}}$, the Majorana phases entering
$U_{\text{PMNS}}$ and the VEV phases $\theta_{L}$ and $\alpha$.

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 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_{\mathrm{min}}\leq 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_{\mathrm{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
$\sim 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_{\mathrm{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_{\mathrm{min}}\lesssim
10^{-3}\,$eV at most by a factor of $R_{im_{\beta\beta}}^{\mathrm{max}}\sim
2.2$ with 76Ge 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_{\beta\beta}}$ despite the
appearance of $C_{VR}^{(6)}$. This ratio would translate to a necessary
accuracy on the nuclear part of the amplitude of $\frac{\Delta
M_{\mathrm{eff}}}{M_{\mathrm{eff}}}\lesssim 14\%$.

Figure 12: Here we present the half-life of 76Ge in the mLRSM model in
dependency on the minimal neutrino mass $m_{\mathrm{min}}$ for both normal
(NO) and inverted ordering (IO). The blue and red areas represent the scenario
marginalized over the unknown phases with the LECs fixed to their order of
magniture estimate while the scattered points show the additional variation of
the relevant LECs. The area inside the black borders represents the usual
mass-mechanism without any additional contributions.

###  4.2 Gluino- and Neutralino-exchange in $\cancel{R}_{p}$ \- SUSY

Supersymmetric theories contain supermultiplets of fermions and bosons which,
under supersymmetry, transform into each other. The most simple constructions
are chiral supermultiplets

$\displaystyle\left(\Psi_{L,R},\Phi_{L,R}^{\Psi}\right)$ (68)

which relate two component chiral spinors $\left(\Psi_{L,R}\right)$ and a
corresponding complex scalar $\Phi_{L,R}$. To construct a supersymmetric
version of the Standard Model, one also needs to consider gauge
supermultiplets

$\displaystyle\left(A_{\mu}^{a},\Psi^{a}\right)$ (69)

which relate the Standard Model’s gauge bosons $A_{\mu}^{a}$ to their
superpartner fermions $\Psi^{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 $\Psi^{a}$ also have 2
degrees of freedom. Therefore, they are Majorana fermions.

$d$$e^{-}$$u$$u$$d$$e^{-}$$\tilde{g}/\chi$$\tilde{d}_{R}$$\tilde{u}_{L}$$d$$u$$e^{-}$$e^{-}$$d$$u$$\tilde{e}_{L}$$\chi$$\tilde{e}_{L}$$d$$e^{-}$$u$$e^{-}$$d$$u$$\,\chi$$\tilde{d}_{R}$$\tilde{e}_{L}$$d$$e^{-}$$u$$e^{-}$$d$$u$$\tilde{u}_{L}$$\chi$$\tilde{e}_{L}$$d$$e^{-}$$u$$d$$e^{-}$$u$$\tilde{g}/\chi$$\tilde{d}_{R}$$\tilde{d}_{R}$$d$$e^{-}$$u$$u$$d$$e^{-}$$\tilde{u}_{L}$$\tilde{g}/\chi$$\tilde{u}_{L}$
Figure 13: Feynman diagrams of gluino- and neutralino-exchange contributing to
$0\nu\beta\beta$ within the $\cancel{R}_{p}$-MSSM [95].

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]

$\displaystyle R_{p}=\left(-1\right)^{2s+3(B-L)}$ (70)

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]

$\displaystyle\begin{split}\mathcal{L}_{\cancel{R}_{p}}^{\Delta
L=1}=-\lambda{{}^{\prime}}_{111}\bigg{[}\left(\overline{u}_{L},\overline{d}_{R}\right)\left(\begin{matrix}e^{c}_{R}\\\
-\nu^{c}_{R}\end{matrix}\right)\tilde{d}_{R}+\left(\overline{e}_{L},\overline{\nu}_{L}\right)d_{R}\left(\begin{matrix}\tilde{u}^{*}_{L}\\\
-\tilde{d}^{*}_{L}\end{matrix}\right)+\left(\overline{u}_{L},\overline{d}_{L}\right)d_{R}\left(\begin{matrix}\tilde{e}^{*}_{L}\\\
-\tilde{\nu}^{*}_{L}\end{matrix}\right)\bigg{]}+\text{h.c.},\end{split}$ (71)

which can contribute to $0\nu\beta\beta$ decay. Contributions to
$0\nu\beta\beta$ decay from $\cancel{R}_{p}-$SUSY have been studied first in
Refs. [95, 97], the corresponding Feynman diagrams are shown in Figure 13. The
relevant gluino ($\tilde{g}$) and neutralino ($\chi$) - fermion interactions
are given by [95, 98, 99]

$\displaystyle\mathcal{L}_{\tilde{g}}$
$\displaystyle=-\sqrt{2}g_{3}\sum_{a}\frac{\lambda_{\alpha\beta}^{(a)}}{2}\left(\overline{q_{L}}^{\alpha}\tilde{g}\tilde{q}_{L}^{\beta}-\overline{q_{R}}^{\alpha}\tilde{g}\tilde{q}_{R}^{\beta}\right)+\text{h.c.}$
(72)

and

$\displaystyle\mathcal{L}_{\chi}$
$\displaystyle=\sqrt{2}g_{2}\sum_{i=1}^{4}\left(\epsilon_{Li}\left(\Psi\right)\overline{\Psi_{L}}\chi_{i}\tilde{\Psi}_{L}+\epsilon_{Ri}\left(\Psi\right)\overline{\Psi_{R}}\chi_{i}\tilde{\Psi}_{R}+\right)+\text{h.c.}.$
(73)

One can obtain the low-energy effective Lagrangian by integrating out the
heavy superfields as well as the Standard Model particles with masses $\gtrsim
m_{W}$. In doing so, one finds the different low-energy effective dimension-9
$\Delta L=2$ operators that contribute to $0\nu\beta\beta$ decay [97]

$\displaystyle\begin{split}\mathcal{L}_{\cancel{R}_{p}}=\frac{G_{F}^{2}}{2m_{N}}\bigg{[}&\left(\eta_{\tilde{g}}+\eta_{\chi}\right)\left(\left[\overline{u}(1+\gamma^{5})d\right]\left[\overline{u}(1+\gamma^{5})d\right]-\frac{1}{4}\left[\overline{u}\sigma^{\mu\nu}(1+\gamma^{5})d\right]\left[\overline{u}\sigma_{\mu\nu}(1+\gamma^{5})d\right]\right)\\\
&+\left(\eta_{\chi\tilde{e}}+\eta_{\tilde{g}}{{}^{\prime}}-\eta_{\chi\tilde{f}}\right)\left[\overline{u}(1+\gamma^{5})d\right]\left[\overline{u}(1+\gamma^{5})d\right]\bigg{]}\;\left[\overline{e}\left(1+\gamma_{5}\right)e^{c}\right].\end{split}$
(74)

These can be matched onto the LEFT basis as

$\displaystyle\begin{split}C_{2R}^{(9)}{{}^{\prime}}&=\frac{2v}{m_{N}}\left[2\eta_{\tilde{g}}+2\eta_{\chi}+\eta_{\chi\tilde{e}}+\eta_{g}{{}^{\prime}}-\eta_{\chi\tilde{f}}\right]\\\
C_{3R}^{(9)}{{}^{\prime}}&=\frac{4v}{m_{N}}\left[\eta_{\tilde{g}}+\eta_{\chi}\right].\end{split}$
(75)

The coupling constants are given in terms of gluino, neutralino and squark
masses as [95]

Figure 14: Half-life ratios resulting from the $\cancel{R}_{p}$-SUSY
contributions to $0\nu\beta\beta$ normalized to the standard mass mechanism.

$\displaystyle\begin{split}\eta_{\tilde{g}}&=\alpha_{s}\Lambda^{2}\frac{m_{N}}{m_{\tilde{g}}}\left[1+\left(\frac{m_{\tilde{d}_{R}}}{m_{\tilde{u}_{L}}}\right)^{4}\right]\\\
\eta_{\tilde{g}}{{}^{\prime}}&=2\alpha_{s}\Lambda^{2}\frac{m_{N}}{m_{\tilde{g}}}\left(\frac{m_{\tilde{d}_{R}}}{m_{\tilde{u}_{L}}}\right)^{2}\\\
\eta_{\chi}&=\frac{3\alpha_{2}}{4}\Lambda^{2}\sum_{i=1}^{4}\frac{m_{N}}{m_{\chi_{i}}}\left[\epsilon_{Ri}^{2}(d)+\epsilon_{Li}^{2}(u)\left(\frac{m_{\tilde{d}_{R}}}{m_{\tilde{u}_{L}}}\right)^{4}\right]\\\
\eta_{\chi\tilde{e}}&=9\alpha_{2}\Lambda^{2}\left(\frac{m_{\tilde{d}_{R}}}{m_{\tilde{e}_{L}}}\right)^{4}\sum_{i=1}^{4}\epsilon_{Li}^{2}(e)\frac{m_{N}}{m_{\chi_{i}}}\\\
\eta_{\chi\tilde{f}}&=\frac{3\alpha_{2}}{2}\Lambda^{2}\left(\frac{m_{\tilde{d}_{R}}}{m_{\tilde{e}_{L}}}\right)^{2}\sum_{i=1}^{4}\frac{m_{N}}{m_{\chi_{i}}}\bigg{[}\epsilon_{Ri}(d)\epsilon_{Li}(e)+\epsilon_{Ri}(u)\epsilon_{Li}(d)\left(\frac{m_{\tilde{e}_{L}}}{m_{\tilde{u}_{L}}}\right)^{2}\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\;\;+\epsilon_{Ri}(u)\epsilon_{Li}(e)\left(\frac{m_{\tilde{d}_{R}}}{m_{\tilde{u}_{L}}}\right)^{2}\bigg{]}\end{split}$
(76)

with

$\displaystyle\Lambda=\frac{\sqrt{2\pi}}{3}\frac{\lambda{{}^{\prime}}_{111}}{G_{F}m_{\tilde{d}_{R}}^{2}}.$
(77)

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}^{(9)}{{}^{\prime}}$
and $C_{3R}^{(9)}{{}^{\prime}}$ 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 $\cancel{R}_{p}$-SUSY model we consider here follows the same
pattern as the scalar short-range operators $C_{2R}^{(9)}$ and
${C_{2R}^{(9)}}^{\prime}$ 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 $\cancel{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_{\tilde{e}_{L}}=410\,$GeV,
$m_{\tilde{q}_{L,R}}=1600\,$GeV with $q\in[u,d]$ and $m_{\tilde{g}}=2260\,$GeV
[100]. We fix the neutralino masses by requiring that the applied EFT
framework holds, which necessitates $m_{\chi_{i}}\geq\Lambda_{\chi}\simeq
2\,\mathrm{GeV}$. For simplicity we take $m_{\chi_{1}}=2\,\mathrm{GeV}$ and
$m_{\chi_{i}}\rightarrow\infty$ for $i\neq 1$. Lighter neutralino masses in
connection to $0\nu\beta\beta$ have also recently been studied [101]. We set
the coupling constant to $\lambda{{}^{\prime}}_{111}=2\times 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 $\cancel{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.

Figure 15: Half-lives in 76Ge for the standard mass mechanism accompanied by
the exchange of heavy neutralinos and gluinos in the $\cancel{R}_{p}$-SUSY
with $\lambda{{}^{\prime}}_{111}=2\times 10^{-4}$,
$m_{\tilde{e}_{L}}=410\,$GeV,
$m_{\tilde{q}_{L,R}}=1600\,$GeV,$m_{\tilde{g}}=2260\,$GeV,
$m_{\chi_{1}}=\Lambda_{\chi}=2\,\mathrm{GeV}$ and
$m_{\chi_{i}}\rightarrow\infty$ for $i\neq 1$

###  4.3 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\nu\beta\beta$ decay arise from
leptoquarks with $Q^{(1)}=\pm 1/3$ (Figure 16 left) and $Q^{(2)}=\pm 2/3$
(Figure 16 right). The full set of renormalizable LQ-fermion interactions is
given by [108]

LQ $(\Omega)$ | $SU(3)_{C}$ | $SU(2)_{L}$ | $U(1)_{Y}$ | Q
---|---|---|---|---
$S_{0}$ | 3 | 1 | -2/3 | -1/3
$\tilde{S}_{0}$ | 3 | 1 | -8/3 | -4/3
$S_{1/2}$ | $\overline{3}$ | 2 | -7/3 | $\left(-2/3,-5/3\right)$
$\tilde{S}_{1/2}$ | $\overline{3}$ | 2 | -1/3 | $\left(1/3,-2/3\right)$
$S_{1}$ | 3 | 3 | -2/3 | $\left(2/3,-1/3,-4/3\right)$
$V_{0}$ | $\overline{3}$ | 1 | -4/3 | -2/3
$\tilde{V}_{0}$ | $\overline{3}$ | 1 | -10/3 | -5/3
$V_{1/2}$ | 3 | 2 | -5/3 | $\left(-1/3,-4/3\right)$
$\tilde{V}_{1/2}$ | 3 | 2 | 1/3 | $\left(2/3,-1/3\right)$
$V_{1}$ | $\overline{3}$ | 3 | -4/3 | $\left(1/3,-2/3,-5/3\right)$
Table 5: List of possible scalar and vector leptoquarks and their
transformation properties under the Standard Model symmetries [108].

$\displaystyle\begin{split}\mathcal{L}_{S,f}=&\left(\lambda_{S_{0}}^{R}\right)_{ij}S_{0}^{R\dagger}\left[\overline{u^{c}}_{i}P_{R}e_{j}\right]+\left(\lambda_{\tilde{S}_{0}}^{R}\right)_{ij}\tilde{S}_{0}^{\dagger}\left[\overline{d^{c}}_{i}P_{R}e_{j}\right]\\\
+&\left(\lambda_{S_{1/2}}^{R}\right)_{ij}S_{1/2}^{R\dagger}\left[\overline{u}_{i}P_{L}L_{j}\right]+\left(\lambda_{\tilde{S}_{1/2}}^{R}\right)_{ij}\tilde{S}_{1/2}^{\dagger}\left[\overline{d}_{i}P_{L}L_{j}\right]\\\
+&\left(\lambda_{S_{0}}^{L}\right)_{ij}S_{0}^{L\dagger}\left[\overline{Q^{c}}_{i}P_{L}i\tau_{2}L_{j}\right]+\left(\lambda_{S_{1/2}}^{L}\right)_{ij}S_{1/2}^{L\dagger}\left[\overline{Q^{c}}_{i}P_{R}i\tau_{2}e_{j}\right]\\\
+&\left(\lambda_{S_{1}}^{L}\right)_{ij}\left[\overline{Q^{c}}_{i}P_{L}i\tau_{2}S_{1}^{\dagger}L_{j}\right]+\text{h.c.}\end{split}$
(78)

and

$\displaystyle\begin{split}\mathcal{L}_{V,f}=&\left(\lambda_{V_{0}}^{R}\right)_{ij}V_{0\mu}^{R\dagger}\left[\overline{d}_{i}\gamma^{\mu}P_{R}e_{j}\right]+\left(\lambda_{\tilde{V}_{0}}^{R}\right)_{ij}\tilde{V}_{0\mu}^{R\dagger}\left[\overline{u}_{i}\gamma^{\mu}P_{R}e_{j}\right]\\\
+&\left(\lambda_{V_{1/2}}^{R}\right)_{ij}V_{1/2\mu}^{R\dagger}\left[\overline{d^{c}}_{i}\gamma^{\mu}P_{L}L_{j}\right]+\left(\lambda_{\tilde{V}_{1/2}}^{R}\right)_{ij}\tilde{V}_{1/2}^{\dagger}\left[\overline{u^{c}}_{i}\gamma^{\mu}P_{L}L_{j}\right]\\\
+&\left(\lambda_{V_{0}}^{L}\right)_{ij}V_{0\mu}^{L\dagger}\left[\overline{Q}_{i}\gamma^{\mu}P_{L}L_{j}\right]+\left(\lambda_{V_{1/2}}^{L}\right)_{ij}V_{1/2\mu}^{L\dagger}\left[\overline{Q^{c}}_{i}\gamma^{\mu}P_{R}e_{j}\right]\\\
+&\left(\lambda_{V_{1}}^{L}\right)_{ij}\left[\overline{Q}_{i}\gamma^{\mu}P_{L}V_{1\mu}^{\dagger}L_{j}\right]+\text{h.c.}\end{split}$
(79)

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,

$\displaystyle\begin{split}\mathcal{L}_{LQ,\Phi}=&h_{S_{0}}^{i}\tilde{\Phi}^{\dagger}\tilde{S}_{1/2}S_{0}^{i}+h_{V_{0}}^{i}\tilde{\Phi}^{\dagger}\tilde{V}_{1/2}^{\mu}V_{0\mu}^{i}\\\
+&h_{S_{1}}\tilde{\Phi}^{\dagger}S_{1}\tilde{S}_{1/2}+h_{V_{1}}\tilde{\Phi}^{\dagger}V_{1}^{\mu}\tilde{V}_{1/2\mu}\\\
+&Y_{S_{1/2}}^{i}\left(\tilde{\Phi}^{\dagger}S_{1/2}^{i}\right)\left(\tilde{S}_{1/2}^{\dagger}\Phi\right)+Y_{V_{1/2}}^{i}\left(\tilde{\Phi}^{\dagger}V_{1/2}^{\mu
i}\right)\left(\tilde{V}_{1/2\mu}^{\dagger}\Phi\right)\\\
+&Y_{S_{1}}\left(\tilde{\Phi}^{\dagger}S_{1}^{\dagger}\Phi\right)\tilde{S}_{0}+Y_{V_{1}}\left(\tilde{\Phi}^{\dagger}V_{1\mu}^{\dagger}\Phi\right)\tilde{V}_{0}^{\mu}\\\
+&\kappa_{S}^{i}\left(\Phi^{\dagger}S_{1}\Phi\right){S_{0}^{i}}^{\dagger}+\kappa_{V}^{i}\left(\Phi^{\dagger}V_{1}^{\mu}\Phi\right){V_{0\mu}^{i}}^{\dagger}+\text{h.c.}\\\
-&\sum_{\Omega}\left(\eta_{\Omega}M^{2}_{\Omega}-g_{\Omega}^{i_{1}i_{2}}\Phi^{\dagger}\Phi\right){\Omega^{i_{1}}}^{\dagger}\Omega^{i_{2}},\end{split}$
(80)

where the leptoquark triplets are defined as

$\displaystyle V_{1}=\sum_{i}\tau_{i}V_{1i}\qquad
S_{1}=\sum_{i}\tau_{i}S_{1i}.$ (81)

These LQ-Higgs interactions are essential when considering contributions to
$0\nu\beta\beta$ decay because they result in non-zero correlation functions
for, e.g.,

$\displaystyle\left<{S_{0}^{i}\tilde{S}_{1/2}}\right>\propto\sum_{\tilde{I}}\mathcal{N}_{S_{0}^{i}\tilde{I}}\mathcal{N}_{\tilde{S}_{1/2}\tilde{I}},$
(82)

where $\mathcal{N}$ is the mixing matrix which diagonalizes the mass matrix
$\mathcal{N}^{T}\mathcal{M}^{2}\mathcal{N}=\mathcal{M}^{2}_{diag}$ and
$\tilde{I}=\mathcal{N}^{T}I$ are the mass eigenstate fields. This particular
example results in contributions captured by the right diagram in Figure 16.

$d$$u$$e^{-}$$e^{-}$$d$$u$$\nu$$V(S)$$W$$d$$e^{-}$$u$$e^{-}$$d$$u$$V(S)$$W$
Figure 16: Feynman diagrams of the vector $(V)$ and scalar $(S)$ leptoquark
interactions contributing to $0\nu\beta\beta$.

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\nu\beta\beta$ decay are then given by [108]

$\displaystyle\begin{split}\mathcal{L}_{LQ}=&\left[\overline{e}P_{L}\nu^{c}\right]\left\\{\frac{\epsilon_{S}}{M^{2}_{S}}\left[\overline{u}P_{R}d\right]+\frac{\epsilon_{V}}{M^{2}_{V}}\left[\overline{u}P_{L}d\right]\right\\}\\\
-&\left[\overline{e}\gamma^{\mu}P_{L}\nu^{c}\right]\left\\{\left(\frac{\alpha_{S}^{R}}{M_{S}^{2}}+\frac{\alpha_{V}^{R}}{M_{V}^{2}}\right)\left[\overline{u}\gamma_{\mu}P_{R}d\right]-\sqrt{2}\left(\frac{\alpha_{S}^{L}}{M_{S}^{2}}+\frac{\alpha_{V}^{L}}{M_{V}^{2}}\right)\left[\overline{u}\gamma_{\mu}P_{L}d\right]\right\\}+\text{h.c.},\end{split}$
(83)

Figure 17: Angular correlations (upper right) and single electron energy
spectra (lower left) resulting from the different LQ contributions as well as
the standard mass mechanism in 136Xe. The unknown LECs are set to their order-
of-magnitude estimates. The specific choice of the isotope does slightly
influence the shape of the angular correlation.

with the low-energy Wilson coefficients

$\displaystyle\epsilon_{I}$
$\displaystyle=2^{-\eta_{I}}\left[\lambda_{I_{1}}^{L}\lambda_{\tilde{I}_{1/2}}^{R}\left(\tilde{\theta}_{43}^{I}\left(Q_{I}^{1}\right)+\eta_{I}\sqrt{2}\tilde{\theta}_{41}^{I}\left(Q_{I}^{2}\right)\right)-\lambda_{I_{0}}^{L}\lambda_{\tilde{I}_{1/2}}^{R}\tilde{\theta}_{23}^{I}\left(Q_{I}^{1}\right)\right]$
(84) $\displaystyle\alpha_{I}^{L}$
$\displaystyle=\frac{2}{3+\eta_{I}}\lambda_{I_{1/2}}^{L}\lambda_{I_{1}}^{L}\tilde{\theta}_{24}^{I}\left(Q_{I}^{2}\right),\qquad\qquad\alpha_{I}^{R}=\frac{2}{r+\eta_{I}}\lambda_{I_{0}}^{R}\lambda_{\tilde{I}_{1/2}}^{R}\tilde{\theta}_{23}^{I}\left(Q_{I}^{1}\right),$
(85)

where

$\displaystyle\tilde{\theta}^{I}_{ij}=\sum_{k}\mathcal{N}_{ik}\mathcal{N}_{jk}\frac{M_{I}^{2}}{M_{I_{k}}^{2}}.$
(86)

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 at Eq. (84) and Eq. (85), there is a priori no reason from, e.g.,
naturalness arguments why any of the low energy coefficients $\alpha_{I}$ and
$\epsilon_{I}$ should be suppressed or enhanced compared to the others.
However, if the LQ interactions arise from a more complete model or if simply
not all possible LQ interactions are realized in nature, hierarchical
structures might appear. We will therefore study different settings in which
some couplings dominate over the others. From Eq. (83) we can match the Wilson
coefficients in Eq. (84) and Eq. (85) onto the LEFT basis arriving at

$\displaystyle\begin{split}C^{(6)}_{SL}&=\frac{v^{2}}{M_{V}^{2}}\epsilon_{V}\,,\qquad\qquad\qquad\qquad
C^{(6)}_{SR}=\frac{v^{2}}{M_{S}^{2}}\epsilon_{S}\\\
C^{(6)}_{VL}&=\sqrt{2}v^{2}\left(\frac{\alpha_{S}^{L}}{M_{S}^{2}}+\frac{\alpha_{V}^{L}}{M_{V}^{2}}\right)\,,\qquad
C^{(6)}_{VR}=-v^{2}\left(\frac{\alpha_{S}^{R}}{M_{S}^{2}}+\frac{\alpha_{V}^{R}}{M_{V}^{2}}\right).\end{split}$
(87)

We study the following 7 different settings of LQ contributions to
$0\nu\beta\beta$ decay:

  1. 1.

Full LQ Model:
$\epsilon_{S}=\epsilon_{V}=\alpha_{S}^{L}=\alpha_{S}^{R}=\alpha_{V}^{L}=\alpha_{V}^{R}=1$

  2. 2.

Scalar LQs (S): $\epsilon_{S}=\alpha_{S}^{L}=\alpha_{S}^{R}=1$

  3. 3.

Scalar LQs coupling to LH fermions (SL): $\alpha_{S}^{L}=1$

  4. 4.

Scalar LQs coupling to RH fermions (SR): $\alpha_{S}^{R}=1$

  5. 5.

Vector LQs (V): $\epsilon_{V}=\alpha_{V}^{L}=\alpha_{V}^{R}=1$

  6. 6.

Vector LQs coupling to LH fermions (VL): $\alpha_{V}^{L}=1$

  7. 7.

Vector LQs coupling to RH fermions (VR): $\alpha_{V}^{R}=1$

Figure 18: Half-life ratios resulting from different leptoquark settings when
taking 76Ge as the reference isotope. The ratios are compared to the standard
mass mechanism.

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 order-of-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_{\beta\beta}}$ is
still significant for the full model as well as the SL and VL models.
Considering the central values, the highest ratio when taking 76Ge as the
reference isotope is realized in the vector model with
$R_{im_{\beta\beta}}^{\mathrm{max}}\sim 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 $\frac{\Delta
M_{\mathrm{eff}}}{M_{\mathrm{eff}}}\lesssim 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.

Figure 19: Here we show the expected half-lives for the full LQ model with the
parameters fixed to
$\epsilon_{S}=\epsilon_{V}=\alpha_{S}^{L}=\alpha_{S}^{R}=\alpha_{V}^{L}=\alpha_{V}^{R}=1$
and the suppression scales $M_{S}=M_{V}=10^{7}\,$GeV

##  5 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 $\leq 9$ triggering
$0\nu\beta\beta$ decay at low energy.

The main aim of our study is to understand the differences in various
$0\nu\beta\beta$ 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_{\beta\beta}}$ 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\nu\beta\beta$
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_{\mathrm{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\nu\beta\beta$ decay even if nuclear uncertainties remain significant. For
operator groups such as $C_{V}{(9)}$ or $\tilde{C}_{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\nu\beta\beta$ 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\nu\beta\beta$ 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\nu\beta\beta$ 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\nu\beta\beta$ 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, $\sum_{i}m_{\nu_{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\nu\beta\beta$ decay. Same methods as employed in this
study may allow to unravel additional contributions to $0\nu\beta\beta$ decay
induced by light sterile neutrinos.

## Acknowledgments

We would like to thank Jordy de Vries for many fruitful discussions and useful
comments. We are also grateful to Frank F. Deppisch for valuable comments on
the final version of the manuscript. LG acknowledges support from the National
Science Foundation, Grant PHY-1630782, and to the Heising-Simons Foundation,
Grant 2017-228.

##  Appendix A Double $\beta$ Modes

$2\nu\beta^{-}\beta^{-}$ | $2\nu\beta^{+}\beta^{+}$ | $2\nu\mathrm{EC}\beta^{+}$ | $2\nu\mathrm{ECEC}$
---|---|---|---
${}^{A}Z$ |  Q [MeV] |  ${}^{A}Z$ |  Q [MeV] |  ${}^{A}Z$ |  Q [MeV] |  ${}^{A}Z$ |  Q [MeV]
46Ca |  0.99 |  78Kr |  0.80 |  50Cr |  0.15 |  36Ar |  0.43
48Ca |  4.27 |  96Ru |  0.67 |  58Ni |  0.90 |  40Ca |  0.19
70Zn |  1.00 |  106Cd |  0.73 |  64Zn |  0.073 |  50Cr |  1.17
76Ge |  2.04 |  124Xe |  0.82 |  74Se |  0.19 |  54Fe |  0.68
80Se |  0.13 |  130Ba |  0.57 |  78Kr |  1.82 |  58Ni |  1.93
82Se |  3.00 |  136Ce |  0.33 |  84Sr |  0.77 |  64Zn |  1.09
86Kr |  1.26 |  |  |  92Mo |  0.63 |  74Se |  1.21
94Zr |  1.14 |  |  |  96Ru |  1.69 |  78Kr |  2.85
96Zr |  3.35 |  |  |  102Pd |  0.15 |  84Sr |  1.79
98Mo |  0.11 |  |  |  106Cd |  1.75 |  92Mo |  1.65
100Mo |  3.03 |  |  |  112Sn |  0.90 |  96Ru |  2.71
104Ru |  1.30 |  |  |  120Te |  0.71 |  102Pd |  1.17
110Pd |  2.02 |  |  |  124Xe |  1.84 |  106Cd |  2.78
114Cd |  0.54 |  |  |  130Ba |  1.60 |  108Cd |  0.27
116Cd |  2.81 |  |  |  136Ce |  1.36 |  112Sn |  1.92
122Sn |  0.37 |  |  |  144Sm |  0.76 |  120Te |  1.73
124Sn |  2.29 |  |  |  156Dy |  0.98 |  124Xe |  2.86
128Te |  0.87 |  |  |  162Er |  0.82 |  126Xe |  0.92
130Te |  2.53 |  |  |  168Yb |  0.39 |  130Ba |  2.62
134Xe |  0.83 |  |  |  174Hf |  0.077 |  132Ba |  0.84
136Xe |  2.46 |  |  |  184Os |  0.43 |  136Ce |  2.38
142Ce |  1.42 |  |  |  190Pt |  0.36 |  138Ce |  0.69
146Nd |  0.070 |  |  |  |  |  144Sm |  1.78
148Nd |  1.93 |  |  |  |  |  152Gd |  0.056
150Nd |  3.37 |  |  |  |  |  156Dy |  2.01
154Sm |  1.25 |  |  |  |  |  158Dy |  0.28
160Gd |  1.73 |  |  |  |  |  162Er |  1.85
170Er |  0.66 |  |  |  |  |  164Er |  0.025
176Yb |  1.09 |  |  |  |  |  168Yb |  1.41
186W |  0.49 |  |  |  |  |  174Hf |  1.10
192Os |  0.41 |  |  |  |  |  180W |  0.14
198Pt |  1.05 |  |  |  |  |  184Os |  1.45
204Hg |  0.42 |  |  |  |  |  190Pt |  1.38
232Th |  0.84 |  |  |  |  |  196Hg |  0.82
238U |  1.14 |  |  |  |  |  |
Table 6: Complete list of natural double-$\beta$ elements and the
corresponding $Q$-values calculated from the NIST list of elements [88] using
the conditions 42, 44, 46 and 48. Overall there are 69 different natural
elements that can decay via at least one double-$\beta$ mode.

In Table 6 we present a list of all naturally occuring isotopes that decay via
any of the $0\nu\beta\beta$-modes i.e.
$0\nu\beta^{-}\beta^{-},0\nu\beta^{+}\beta^{+},0\nu\mathrm{EC}\beta^{+}$ and
$0\nu\mathrm{ECEC}$. The isotopes as well as the corresponding Q-values are
taken and calculated from the NIST list of elements [88].

##  Appendix 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.

$\displaystyle\begin{split}m_{\beta\beta}:\qquad
T_{1/2}^{-1}=&g_{A}^{4}\left|m_{\beta\beta}\right|^{2}G_{01}\left|\frac{V_{ud}^{2}}{m_{e}}\left(\frac{1}{g_{A}^{2}}M_{F}-M_{GT}-2\frac{m_{\pi}^{2}\mathbf{g_{\nu}^{NN}}}{g_{A}^{2}}M_{F,sd}\right)\right|^{2}\end{split}$
(88) $\displaystyle\begin{split}C_{VL}^{(6)}:\qquad
T_{1/2}^{-1}=&g_{A}^{4}\left|C_{VL}^{(6)}\right|^{2}\Bigg{\\{}4G_{02}\left|\frac{V_{ud}}{3}\left(\frac{g_{V}^{2}}{g_{A}^{2}}M_{F}+\frac{1}{3}\left(2M_{GT}^{AA}+M_{T}^{AA}\right)+\frac{6\mathbf{g_{VL}^{E}}}{g_{A}^{2}}M_{F,sd}\right)\right|^{2}\\\
&\qquad\qquad\,+2G_{04}\bigg{|}\frac{V_{ud}}{6}\bigg{(}\frac{g_{V}^{2}}{g_{A}^{2}}M_{F}-\frac{1}{3}\left(M_{GT}^{AA}-4M_{T}^{AA}\right)\\\
&\qquad\qquad\qquad\qquad-3\left(M_{GT}^{AP}+M_{GT}^{PP}+M_{T}^{AP}+M_{T}^{PP}\right)-\frac{12\mathbf{g_{VL}^{m_{e}}}}{g_{A}^{2}}M_{F,sd}\bigg{)}\bigg{|}^{2}\\\
&\qquad\qquad\,-4G_{03}\text{Re}\bigg{[}\frac{\left|V_{ud}\right|^{2}}{18}\left(\frac{g_{V}^{2}}{g_{A}^{2}}M_{F}+\frac{1}{3}\left(2M_{GT}^{AA}+M_{T}^{AA}\right)+\frac{6\mathbf{g_{VL}^{E}}}{g_{A}^{2}}M_{F,sd}\right)^{*}\\\
&\qquad\qquad\qquad\qquad\bigg{(}\frac{g_{V}^{2}}{g_{A}^{2}}M_{F}-\frac{1}{3}\left(M_{GT}^{AA}-4M_{T}^{AA}\right)\\\
&\qquad\qquad\qquad\qquad-3\left(M_{GT}^{AP}+M_{GT}^{PP}+M_{T}^{AP}+M_{T}^{PP}\right)-\frac{12\mathbf{g_{VL}^{m_{e}}}}{g_{A}^{2}}M_{F,sd}\bigg{)}\bigg{]}\\\
&\qquad\qquad\,+G_{09}\bigg{|}\frac{m_{N}}{m_{e}}V_{ud}\bigg{[}2\frac{g_{A}}{g_{M}}\left(M_{GT}^{MM}+M_{T}^{MM}\right)\\\
&\qquad\qquad\qquad\qquad+\frac{m_{\pi}^{2}}{m_{N}^{2}}\left(-\frac{2}{g_{A}^{2}}\mathbf{g_{VL}^{NN}}M_{F,sd}+\frac{1}{2}\mathbf{g_{VL}^{\pi
N}}\left(M_{GT,sd}^{AP}+M_{T,sd}^{AP}\right)\right)\bigg{]}\bigg{|}^{2}\Bigg{\\}}\end{split}$
(89) $\displaystyle\begin{split}C_{VR}^{(6)}:\qquad
T_{1/2}^{-1}=&g_{A}^{4}\left|C_{VR}^{(6)}\right|^{2}\Bigg{\\{}4G_{02}\bigg{|}\frac{V_{ud}}{3}\left(\frac{g_{V}^{2}}{g_{A}^{2}}M_{F}-\frac{1}{3}\left(2M_{GT}^{AA}+M_{T}^{AA}\right)+6\frac{\mathbf{g_{VR}^{E}}}{g_{A}^{2}}M_{F,sd}\right)\bigg{|}^{2}\\\
&\qquad\qquad\,+2G_{04}\bigg{|}\frac{V_{ud}}{6}\bigg{(}\frac{g_{V}^{2}}{g_{A}^{2}}M_{F}+\frac{1}{3}\left(M_{GT}^{AA}-4M_{T}^{AA}\right)\\\
&\qquad\qquad\qquad\qquad+3\left(M_{GT}^{AP}+M_{GT}^{PP}+M_{T}^{AP}+M_{T}^{PP}\right)-12\frac{\mathbf{g_{VR}^{me}}}{g_{A}^{2}}M_{F,sd}\bigg{)}\bigg{|}^{2}\\\
&\qquad\qquad\,+2G_{03}\text{Re}\bigg{[}\frac{\left|V_{ud}\right|^{2}}{18}\left(\frac{g_{V}^{2}}{g_{A}^{2}}M_{F}-\frac{1}{3}\left(2M_{GT}^{AA}+M_{T}^{AA}\right)+6\frac{\mathbf{g_{VR}^{E}}}{g_{A}^{2}}M_{F,sd}\right)^{*}\\\
&\qquad\qquad\qquad\qquad\quad\bigg{(}\frac{g_{V}^{2}}{g_{A}^{2}}M_{F}+\frac{1}{3}\left(M_{GT}^{AA}-4M_{T}^{AA}\right)\\\
&\qquad\qquad\qquad\qquad\quad+3\left(M_{GT}^{AP}+M_{GT}^{PP}+M_{T}^{AP}+M_{T}^{PP}\right)-12\frac{\mathbf{g_{VR}^{m_{e}}}}{g_{A}^{2}}M_{F,sd}\bigg{)}\bigg{]}\Bigg{\\}}\end{split}$
(90) $\displaystyle\begin{split}C_{SL,SR}^{(6)}:\quad
T_{1/2}^{-1}=&g_{A}^{4}\left|C_{SL,SR}^{(6)}\right|^{2}G_{01}\left|\frac{B}{m_{e}}V_{ud}M_{PS}\right|^{2}\end{split}$
(91) $\displaystyle\begin{split}C_{T}^{(6)}:\qquad
T_{1/2}^{-1}=&g_{A}^{4}\left|C_{T}^{(6)}\right|^{2}G_{01}\left|V_{ud}\left[2\frac{\mathbf{g_{T}^{\prime}}-\mathbf{g_{T}^{NN}}}{g_{A}^{2}}\frac{m_{\pi}^{2}}{m_{N}^{2}}M_{F,sd}-\frac{8g_{T}}{g_{M}}\left(M_{GT}^{MM}+M_{T}^{MM}\right)\right]\right|^{2}\end{split}$
(92) $\displaystyle\begin{split}C_{VL,VR}^{(7)}:\quad
T_{1/2}^{-1}=&g_{A}^{4}\left|C_{VL,VR}^{(7)}\right|^{2}G_{01}\left|\frac{m_{\pi}^{2}}{m_{e}v}V_{ud}M_{PS}\right|^{2}\end{split}$
(93) $\displaystyle\begin{split}C_{1L,1R}^{(9)}{(^{\prime})}:\quad
T_{1/2}^{-1}=&g_{A}^{4}\left|C_{1L,1R}^{(9)}{(^{\prime})}\right|^{2}G_{01}\bigg{|}\frac{5g_{1}^{\pi\pi}}{6m_{N}^{2}}\left(\frac{1}{2}M_{GT,sd}^{AP}+M_{GT,sd}^{PP}+\frac{1}{2}M_{T,sd}^{AP}+M_{T,sd}^{PP}\right)\\\
&\qquad\qquad\qquad\qquad\quad+\left(\mathbf{g_{1}^{\pi
N}}-\frac{5}{6}g_{1}^{\pi\pi}\right)\frac{m_{\pi}^{2}}{2m_{N}^{2}}\left(M_{GT,sd}^{AP}+M_{T,sd}^{AP}\right)\\\
&\qquad\qquad\qquad\qquad\quad-2\frac{\mathbf{g_{1}^{NN}}}{g_{A}^{2}}\frac{m_{\pi}^{2}}{m_{N}^{2}}M_{F,sd}\bigg{|}^{2}\end{split}$
(94) $\displaystyle\begin{split}C_{2,3^{,}L,R}^{(9)}{(^{\prime})}:\quad
T_{1/2}^{-1}=&g_{A}^{4}\left|C_{2,3^{,}L,R}^{(9)}{(^{\prime})}\right|^{2}G_{01}\bigg{|}\frac{5g_{2,3}^{\pi\pi}}{6m_{N}^{2}}\left(\frac{1}{2}M_{GT,sd}^{AP}+M_{GT,sd}^{PP}+\frac{1}{2}M_{T,sd}^{AP}+M_{T,sd}^{PP}\right)\\\
&\qquad\qquad\qquad\qquad\quad+2\frac{\mathbf{g_{2,3}^{NN}}}{g_{A}^{2}}\frac{m_{\pi}^{2}}{m_{N}^{2}}M_{F,sd}\bigg{|}\end{split}$
(95) $\displaystyle\begin{split}C_{4,5^{,}L,R}^{(9)}:\quad
T_{1/2}^{-1}=&g_{A}^{4}\left|C_{4,5^{,}L,R}^{(9)}\right|^{2}G_{01}\bigg{|}\frac{5g_{4,5}^{\pi\pi}}{6m_{N}^{2}}\left(\frac{1}{2}M_{GT,sd}^{AP}+M_{GT,sd}^{PP}+\frac{1}{2}M_{T,sd}^{AP}+M_{T,sd}^{PP}\right)\\\
&\qquad\qquad\qquad\qquad\quad-2\frac{\mathbf{g_{4,5}^{NN}}}{g_{A}^{2}}\frac{m_{\pi}^{2}}{m_{N}^{2}}M_{F,sd}\bigg{|}\end{split}$
(96) $\displaystyle\begin{split}C_{V}^{(9)}:\qquad
T_{1/2}^{-1}=&g_{A}^{4}\left|C_{V}^{(9)}\right|^{2}G_{09}\left|\frac{m_{\pi}^{2}}{m_{e}v}\left(-\frac{2}{g_{A}}\mathbf{g_{6}^{NN}}M_{F,sd}+\frac{1}{2}\mathbf{g_{V}^{\pi
N}}\left(M_{GT,sd}^{AP}+M_{T,sd}^{AP}\right)\right)\right|^{2}\end{split}$
(97) $\displaystyle\begin{split}\tilde{C}_{V}^{(9)}:\qquad
T_{1/2}^{-1}=&g_{A}^{4}\left|\tilde{C}_{V}^{(9)}\right|^{2}G_{09}\left|\frac{m_{\pi}^{2}}{m_{e}v}\left(-\frac{2}{g_{A}}\mathbf{g_{7}^{NN}}M_{F,sd}+\frac{1}{2}\mathbf{\tilde{g}_{V}^{\pi
N}}\left(M_{GT,sd}^{AP}+M_{T,sd}^{AP}\right)\right)\right|^{2}\end{split}$
(98)

##  Appendix C Considering all Isotopes

While we have focussed our discussion on isotopes for which experimental
limits on the half-lives exist, we want to present our main findings of Figure
8 here again but now considering all naturally occuring $0\nu\beta\beta$
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.

Figure 20: Same as Figure 8 but now for all isotopes with available NMEs in
the IBM2 framework: The maximal ratios $R^{\mathrm{max}}_{ij}$ for all
operator combinations i,j are shown. The exact values and the corresponding
isotopes are displayed in each tile. Additionally, operator combinations that
result in different phase-space observables are marked by dashed-line shading.
In the upper right half of the plot we show the ratios considering the central
values from the variation of the LECs. In the lower left half we show the
worst-case scenario considering the values of ratios $R_{ij}$ that are closest
to 1 within the range obtained by the variation of the LECs.

Figure 21: The decay rate ratios $R^{O_{i}}$ (upper plot) and
$R_{im_{\beta\beta}}$ (lower plot) for the different operator groups are
shown. The larger markers represent the choice of vanishing unknown LECs with
$g_{6,7}^{NN}=g_{V}^{\pi N}=\tilde{g}_{V}^{\pi N}=1$. Isotopes with a PSF
$G_{0}>10^{-14}\,\mathrm{y^{-1}}$ are represented by stars while isotopes with
smaller PSFs are represented by round markers. The additional points represent
variations of the different unknown LECs $g_{i}$ randomly chosen from
$\left[-\sqrt{10},-1/\sqrt{10}\right]\times|g_{i}|$ and
$\left[+1/\sqrt{10},+\sqrt{10}\right]\times|g_{i}|$ except for $g_{\nu}^{NN}$
which is varied in a range of $\pm 50\%$. The crosses represent the central
values of the variation i.e. the median values. The reference isotope is
chosen to be 76Ge.

Figure 22: The decay rate ratios $R^{O_{i}}$ (upper plot) and
$R_{im_{\beta\beta}}$ (lower plot) for the different operator groups in the
$\epsilon$-basis.

## References

  * [1] W. Furry, “On transition probabilities in double beta-disintegration,” Phys. Rev. 56 (1939) 1184–1193.
  * [2] J. Schechter and J. W. F. Valle, “Neutrinoless double-$\beta$ decay in su(2)$\times{}$u(1) theories,” Phys. Rev. D 25 (Jun, 1982) 2951–2954. https://link.aps.org/doi/10.1103/PhysRevD.25.2951.
  * [3] H. Pas, M. Hirsch, H. Klapdor-Kleingrothaus, and S. Kovalenko, “Towards a superformula for neutrinoless double beta decay,” Phys. Lett. B 453 (1999) 194–198.
  * [4] H. Pas, M. Hirsch, H. Klapdor-Kleingrothaus, and S. Kovalenko, “A Superformula for neutrinoless double beta decay. 2. The Short range part,” Phys. Lett. B 498 (2001) 35–39, arXiv:hep-ph/0008182.
  * [5] F. F. Deppisch, M. Hirsch, and H. Päs, “Neutrinoless Double Beta Decay and Physics Beyond the Standard Model,” J. Phys. G 39 (2012) 124007, arXiv:1208.0727 [hep-ph].
  * [6] V. Cirigliano, W. Dekens, J. de Vries, M. L. Graesser, and E. Mereghetti, “Neutrinoless double beta decay in chiral effective field theory: lepton number violation at dimension seven,” Journal of High Energy Physics 2017 no. 12, (Dec, 2017) . http://dx.doi.org/10.1007/JHEP12(2017)082.
  * [7] V. Cirigliano, W. Dekens, J. de Vries, M. L. Graesser, and E. Mereghetti, “A neutrinoless double beta decay master formula from effective field theory,” Journal of High Energy Physics 2018 no. 12, (Dec, 2018) . http://dx.doi.org/10.1007/JHEP12(2018)097.
  * [8] L. Graf, F. F. Deppisch, F. Iachello, and J. Kotila, “Short-Range Neutrinoless Double Beta Decay Mechanisms,” Phys. Rev. D 98 no. 9, (2018) 095023, arXiv:1806.06058 [hep-ph].
  * [9] F. F. Deppisch, L. Graf, F. Iachello, and J. Kotila, “Analysis of light neutrino exchange and short-range mechanisms in $0\nu\beta\beta$ decay,” Phys. Rev. D 102 no. 9, (2020) 095016, arXiv:2009.10119 [hep-ph].
  * [10] J. C. Pati and A. Salam, “Lepton Number as the Fourth Color,” Phys. Rev. D 10 (1974) 275–289. [Erratum: Phys.Rev.D 11, 703–703 (1975)].
  * [11] R. N. Mohapatra and J. C. Pati, “Left-Right Gauge Symmetry and an Isoconjugate Model of CP Violation,” Phys. Rev. D 11 (1975) 566–571.
  * [12] R. N. Mohapatra and J. C. Pati, “A Natural Left-Right Symmetry,” Phys. Rev. D 11 (1975) 2558.
  * [13] G. Senjanovic and R. N. Mohapatra, “Exact Left-Right Symmetry and Spontaneous Violation of Parity,” Phys. Rev. D 12 (1975) 1502.
  * [14] W.-C. Huang and J. Lopez-Pavon, “On neutrinoless double beta decay in the minimal left-right symmetric model,” Eur. Phys. J. C 74 (2014) 2853, arXiv:1310.0265 [hep-ph].
  * [15] C. Giunti and M. Laveder, “Short-Baseline Electron Neutrino Disappearance, Tritium Beta Decay and Neutrinoless Double-Beta Decay,” Phys. Rev. D 82 (2010) 053005, arXiv:1005.4599 [hep-ph].
  * [16] J. Barry, W. Rodejohann, and H. Zhang, “Light Sterile Neutrinos: Models and Phenomenology,” JHEP 07 (2011) 091, arXiv:1105.3911 [hep-ph].
  * [17] J. Barea, J. Kotila, and F. Iachello, “Limits on sterile neutrino contributions to neutrinoless double beta decay,” Phys. Rev. D 92 (2015) 093001, arXiv:1509.01925 [hep-ph].
  * [18] W. Dekens, J. de Vries, K. Fuyuto, E. Mereghetti, and G. Zhou, “Sterile neutrinos and neutrinoless double beta decay in effective field theory,” JHEP 06 (2020) 097, arXiv:2002.07182 [hep-ph].
  * [19] G. Li, M. Ramsey-Musolf, and J. C. Vasquez, “Left-Right Symmetry and Leading Contributions to Neutrinoless Double Beta Decay,” Phys. Rev. Lett. 126 no. 15, (2021) 151801, arXiv:2009.01257 [hep-ph].
  * [20] P. D. Bolton, F. F. Deppisch, and P. S. Bhupal Dev, “Neutrinoless double beta decay versus other probes of heavy sterile neutrinos,” JHEP 03 (2020) 170, arXiv:1912.03058 [hep-ph].
  * [21] M. Agostini et al., “Final results of gerda on the search for neutrinoless double-$\beta$ decay,” Physical Review Letters 125 no. 25, (Dec, 2020) . http://dx.doi.org/10.1103/PhysRevLett.125.252502.
  * [22] O. Azzolini et al., “Final result of cupid-0 phase-i in the search for the se82 neutrinoless double-$\beta$ decay,” Physical Review Letters 123 no. 3, (Jul, 2019) . http://dx.doi.org/10.1103/PhysRevLett.123.032501.
  * [23] J. Argyriades et al., “Measurement of the two neutrino double beta decay half-life of zr-96 with the nemo-3 detector,” Nuclear Physics A 847 no. 3-4, (Dec, 2010) 168–179. http://dx.doi.org/10.1016/j.nuclphysa.2010.07.009.
  * [24] E. Armengaud et al., “New limit for neutrinoless double-beta decay of mo100 from the cupid-mo experiment,” Physical Review Letters 126 no. 18, (May, 2021) . http://dx.doi.org/10.1103/PhysRevLett.126.181802.
  * [25] F. A. Danevich et al., “Search for double beta decay of116cd with enriched116cdwo4crystal scintillators (aurora experiment),” Journal of Physics: Conference Series 718 (May, 2016) 062009. http://dx.doi.org/10.1088/1742-6596/718/6/062009.
  * [26] C. Arnaboldi et al., “A calorimetric search on double beta decay of 130te,” Physics Letters B 557 no. 3-4, (Apr, 2003) 167–175. http://dx.doi.org/10.1016/S0370-2693(03)00212-0.
  * [27] D. Adams et al., “Improved limit on neutrinoless double-beta decay in te130 with cuore,” Physical Review Letters 124 no. 12, (Mar, 2020) . http://dx.doi.org/10.1103/PhysRevLett.124.122501.
  * [28] G. Anton et al., “Search for neutrinoless double-$\beta$ decay with the complete exo-200 dataset,” Physical Review Letters 123 no. 16, (Oct, 2019) . http://dx.doi.org/10.1103/PhysRevLett.123.161802.
  * [29] A. Gando et al., “Search for majorana neutrinos near the inverted mass hierarchy region with kamland-zen,” Physical Review Letters 117 no. 8, (Aug, 2016) . http://dx.doi.org/10.1103/PhysRevLett.117.082503.
  * [30] J. Argyriades et al., “Measurement of the double-$\beta$ decay half-life of nd150 and search for neutrinoless decay modes with the nemo-3 detector,” Physical Review C 80 no. 3, (Sep, 2009) . http://dx.doi.org/10.1103/PhysRevC.80.032501.
  * [31] KamLAND-Zen Collaboration, S. Abe et al., “First Search for the Majorana Nature of Neutrinos in the Inverted Mass Ordering Region with KamLAND-Zen,” arXiv:2203.02139 [hep-ex].
  * [32] N. Abgrall et al., “The large enriched germanium experiment for neutrinoless double beta decay (legend),” AIP Conference Proceedings (2017) . http://dx.doi.org/10.1063/1.5007652.
  * [33] LEGEND Collaboration, N. Abgrall et al., “Legend-1000 preconceptual design report,” 2021.
  * [34] E. Armengaud et al., “The cupid-mo experiment for neutrinoless double-beta decay: performance and prospects,” The European Physical Journal C 80 no. 1, (Jan, 2020) . http://dx.doi.org/10.1140/epjc/s10052-019-7578-6.
  * [35] V. Albanese et al., “The sno+ experiment,” Journal of Instrumentation 16 no. 08, (Aug, 2021) P08059. http://dx.doi.org/10.1088/1748-0221/16/08/P08059.
  * [36] nEXO Collaboration Collaboration, G. Adhikari et al., “nexo: Neutrinoless double beta decay search beyond $10^{28}$ year half-life sensitivity,” 2021.
  * [37] F. F. Deppisch, L. Graf, and F. Šimkovic, “Searching for New Physics in Two-Neutrino Double Beta Decay,” Phys. Rev. Lett. 125 no. 17, (2020) 171801, arXiv:2003.11836 [hep-ph].
  * [38] M. J. Dolinski, A. W. Poon, and W. Rodejohann, “Neutrinoless double-beta decay: Status and prospects,” Annual Review of Nuclear and Particle Science 69 no. 1, (Oct, 2019) 219–251. http://dx.doi.org/10.1146/annurev-nucl-101918-023407.
  * [39] M. Agostini, G. Benato, J. A. Detwiler, J. Menéndez, and F. Vissani, “Toward the discovery of matter creation with neutrinoless double-beta decay,” arXiv:2202.01787 [hep-ex].
  * [40] J. de Vries, L. Gráf, and O. Scholer To be published.
  * [41] M. L. Graesser, “An electroweak basis for neutrinoless double $\beta$ decay,” Journal of High Energy Physics 2017 no. 8, (Aug, 2017) . http://dx.doi.org/10.1007/JHEP08(2017)099.
  * [42] G. Prezeau, M. Ramsey-Musolf, and P. Vogel, “Neutrinoless double beta decay and effective field theory,” Phys. Rev. D 68 (2003) 034016, arXiv:hep-ph/0303205.
  * [43] F. F. Deppisch, L. Graf, F. Iachello, and J. Kotila, “Analysis of light neutrino exchange and short-range mechanisms in $0\nu\beta\beta$ decay,” Physical Review D 102 no. 9, (Nov, 2020) . http://dx.doi.org/10.1103/PhysRevD.102.095016.
  * [44] V. Cirigliano, W. Dekens, J. De Vries, M. L. Graesser, E. Mereghetti, S. Pastore, and U. Van Kolck, “New Leading Contribution to Neutrinoless Double-$\beta$ Decay,” Phys. Rev. Lett. 120 no. 20, (2018) 202001, arXiv:1802.10097 [hep-ph].
  * [45] V. Cirigliano, W. Dekens, J. De Vries, M. L. Graesser, E. Mereghetti, S. Pastore, M. Piarulli, U. Van Kolck, and R. B. Wiringa, “Renormalized approach to neutrinoless double- $\beta$ decay,” Phys. Rev. C 100 no. 5, (2019) 055504, arXiv:1907.11254 [nucl-th].
  * [46] T. Bhattacharya, V. Cirigliano, S. D. Cohen, R. Gupta, H.-W. Lin, and B. Yoon, “Axial, scalar, and tensor charges of the nucleon from2+1+1-flavor lattice qcd,” Physical Review D 94 no. 5, (Sep, 2016) . http://dx.doi.org/10.1103/PhysRevD.94.054508.
  * [47] A. Nicholson et al., “Heavy physics contributions to neutrinoless double beta decay from qcd,” Physical Review Letters 121 no. 17, (Oct, 2018) . http://dx.doi.org/10.1103/PhysRevLett.121.172501.
  * [48] V. Cirigliano, W. Dekens, J. de Vries, M. Hoferichter, and E. Mereghetti, “Toward complete leading-order predictions for neutrinoless double$\beta$ decay,” Physical Review Letters 126 no. 17, (Apr, 2021) . http://dx.doi.org/10.1103/PhysRevLett.126.172002.
  * [49] V. Cirigliano, W. Dekens, J. de Vries, M. Hoferichter, and E. Mereghetti, “Determining the leading-order contact term in neutrinoless double $\beta$ decay,” Journal of High Energy Physics 2021 no. 5, (May, 2021) . http://dx.doi.org/10.1007/JHEP05(2021)289.
  * [50] R. Wirth, J. M. Yao, and H. Hergert, “Ab initio calculation of the contact operator contribution in the standard mechanism for neutrinoless double beta decay,” arXiv:2105.05415 [nucl-th].
  * [51] R. Arnold and others., “Measurement of the $2\nu\beta\beta$ decay half-life and search for the $0\nu\beta\beta$ decay of cd116 with the nemo-3 detector,” Physical Review D 95 no. 1, (Jan, 2017) . http://dx.doi.org/10.1103/PhysRevD.95.012007.
  * [52] R. Arnold et al., “Probing new physics models of neutrinoless double beta decay with supernemo,” The European Physical Journal C 70 no. 4, (Dec, 2010) 927–943. https://doi.org/10.1140/epjc/s10052-010-1481-5.
  * [53] J. Kotila and F. Iachello, “Phase-space factors for double-$\beta$ decay,” Physical Review C 85 no. 3, (Mar, 2012) . http://dx.doi.org/10.1103/PhysRevC.85.034316.