#  USTC-ICTS/PCFT-21-40  Decomposition of $d=9$ short-range $0\nu\beta\beta$
decay operators at one-loop level

Ping-Tao Chen1,2111E-mail<EMAIL_ADDRESS>Gui-Jun Ding1,2222E-mail:
<EMAIL_ADDRESS>Chang-Yuan Yao3333E-mail<EMAIL_ADDRESS>
1Interdisciplinary Center for Theoretical Study and Department of Modern
Physics,
University of Science and Technology of China, Hefei, Anhui 230026, China
2Peng Huanwu Center for Fundamental Theory, Hefei, Anhui 230026, China 3School
of Physics, Nankai University, Tianjin 300071, China

###### Abstract

We perform a systematical study of the dimension-9 short-range
$0\nu\beta\beta$ decay operators at one-loop level. There are only six genuine
topologies which generate eight diagrams, and the recipe to identify the
possible one-loop realizations of the $0\nu\beta\beta$ decay operators is
sketched. Certain hypercharge assignments are excluded by the absence of tree-
level diagrams in a genuine one-loop model. The mediators of each
decomposition can generate Majorana neutrino masses which are discussed up to
two-loop level. We present an example of $0\nu\beta\beta$ decay model in which
the neutrino masses are generated at two-loop level, and the short-range
contribution can be comparable with the mass mechanism in some region of
parameter space.

##  1 Introduction

The discovery of neutrino oscillations shows that neutrinos have tiny masses.
In order to generate the tiny neutrino masses, one must extend the filed
content or gauge symmetry of the standard model, or depart from
renormalizability [1]. Depending on whether neutrinos are Dirac or Majorana
particles, the neutrino mass terms can be constructed in different ways. For
the case of Majorana neutrinos, neutrinos are identical with their
antiparticles and the lepton number would be violated by two units $\Delta
L=2$. As a consequence, the exchange of light Majorana neutrino can lead to
the well-known neutrinoless double beta ($0\nu\beta\beta$) decay

$(A,Z)\rightarrow(A,Z+2)+e^{-}+e^{-}\,,$ (1)

which violates the lepton number by two units as well. This is the standard
mass mechanism. Conversely, if the rare $0\nu\beta\beta$ decay is observed,
neutrinos must be Majorana particles [2]. The new physics giving rise to
$0\nu\beta\beta$ decay will always produce a non-zero mass of Majorana
neutrinos. Under the assumption that the Majorana neutrino mass is the unique
source of $0\nu\beta\beta$ decay, the decay rate is fixed by the phase space
factor, the nuclear matrix element and the effective Majorana mass
$m_{\beta\beta}$ with

$m_{\beta\beta}=\sum^{3}_{i}U_{ei}^{2}m_{i}\,,$ (2)

where $U$ is the lepton mixing matrix and $m_{i}(i=1,2,3)$ are the light
neutrino masses. The effective mass $m_{\beta\beta}$ is the (11) entry of the
light neutrino mass matrix in the basis where the charged lepton mass matrix
is diagonal. Hence $0\nu\beta\beta$ decay can provide valuable information on
both neutrino mass spectrum and Majorana CP phases. The experimental search
for the $0\nu\beta\beta$ decay is extremely challenging, and no positive
signal has been observed so far. The most stringent bound is given by the
KamLAND-Zen Collaboration with the half-life $T^{0\nu}_{1/2}>1.07\times
10^{26}$ years for 136Xe at $90\%$ confidence level, this bound implies a
limit on the effective Majorana neutrino mass $|m_{\beta\beta}|<(61-165)$ meV
[3]. The next generation tonne-scale $0\nu\beta\beta$ decay experiments such
as LEGEND [4] and nEXO [5] aim at the half-life sensitivity of approximately
$10^{28}$ years. If all the future results from $0\nu\beta\beta$ decays,
neutrino oscillation experiments and cosmology observation on the light
neutrino mass sum are compatible with each other, it would be a strong support
to the mass mechanism mentioned above.

However, it is notable that the $0\nu\beta\beta$ decay could possibly be
dominantly induced by other new physics effects beyond that of Majorana
neutrino masses. The contributions to the $0\nu\beta\beta$ decay rate can be
generally divided into the short-range part [6] and long-range part [7]. See
Refs. [8, 9, 10, 11, 12] for review. The short-range contributions arise from
the lepton number violating couplings or the Majorana mass terms of new
particles in possible extensions of the Standard Model (SM). The short-range
dynamics of the $0\nu\beta\beta$ decay is expected to be at a scale $\Lambda$
higher than the weak scale, and it can be described by a set of $\Delta L=2$
and baryon conserving operators which are invariant under the SM gauge
symmetry $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}$ . At leading order, there
are 11 dimension-9 operators which give rise to $0\nu\beta\beta$ decay [13,
14, 15]:

$\displaystyle\mathcal{O}^{0\nu}_{1}$ $\displaystyle=$
$\displaystyle\epsilon_{ij}(\overline{Q}_{i}\gamma^{\mu}Q_{m})(\overline{u}_{R}\gamma_{\mu}d_{R})(\overline{\ell}_{j}\ell^{c}_{m})\,,$
(3a) $\displaystyle\mathcal{O}^{0\nu}_{2}$ $\displaystyle=$
$\displaystyle\epsilon_{ij}(\overline{Q}_{i}\gamma^{\mu}\lambda^{A}Q_{m})(\overline{u}_{R}\gamma_{\mu}\lambda^{A}d_{R})(\overline{\ell}_{j}\ell^{c}_{m})\,,$
(3b) $\displaystyle\mathcal{O}^{0\nu}_{3}$ $\displaystyle=$
$\displaystyle(\overline{u}_{R}Q_{i})(\overline{u}_{R}Q_{j})(\overline{\ell}_{i}\ell^{c}_{j})\,,$
(3c) $\displaystyle\mathcal{O}^{0\nu}_{4}$ $\displaystyle=$
$\displaystyle(\overline{u}_{R}\lambda^{A}Q_{i})(\overline{u}_{R}\lambda^{A}Q_{j})(\overline{\ell}_{i}\ell^{c}_{j})\,,$
(3d) $\displaystyle\mathcal{O}^{0\nu}_{5}$ $\displaystyle=$
$\displaystyle\epsilon_{ij}\epsilon_{mn}(\overline{Q}_{i}d_{R})(\overline{Q}_{m}d_{R})(\overline{\ell}_{j}{\ell}^{c}_{n})\,,$
(3e) $\displaystyle\mathcal{O}^{0\nu}_{6}$ $\displaystyle=$
$\displaystyle\epsilon_{ij}\epsilon_{mn}(\overline{Q}_{i}\lambda^{A}d_{R})(\overline{Q}_{m}\lambda^{A}d_{R})(\overline{\ell}_{j}{\ell}^{c}_{n})\,,$
(3f) $\displaystyle\mathcal{O}^{0\nu}_{7}$ $\displaystyle=$
$\displaystyle(\overline{u}_{R}\gamma^{\mu}d_{R})(\overline{u}_{R}\gamma_{\mu}d_{R})(\overline{e}_{R}e^{c}_{R})\,,$
(3g) $\displaystyle\mathcal{O}^{0\nu}_{8}$ $\displaystyle=$
$\displaystyle\epsilon_{ij}(\overline{u}_{R}\gamma^{\mu}d_{R})(\overline{Q}_{i}d_{R})(\overline{\ell}_{j}\gamma_{\mu}e^{c}_{R})\,,$
(3h) $\displaystyle\mathcal{O}^{0\nu}_{9}$ $\displaystyle=$
$\displaystyle\epsilon_{ij}(\overline{u}_{R}\gamma^{\mu}\lambda^{A}d_{R})(\overline{Q}_{i}\lambda^{A}d_{R})(\overline{\ell}_{j}\gamma_{\mu}e^{c}_{R})\,,$
(3i) $\displaystyle\mathcal{O}^{0\nu}_{10}$ $\displaystyle=$
$\displaystyle(\overline{u}_{R}\gamma^{\mu}d_{R})(\overline{u}_{R}Q_{i})(\overline{\ell}_{i}\gamma_{\mu}e^{c}_{R})\,,$
(3j) $\displaystyle\mathcal{O}^{0\nu}_{11}$ $\displaystyle=$
$\displaystyle(\overline{u}_{R}\gamma^{\mu}\lambda^{A}d_{R})(\overline{u}_{R}\lambda^{A}Q_{i})(\overline{\ell}_{i}\gamma_{\mu}e^{c}_{R})\,,$
(3k)

where $Q=(u_{L},~{}d_{L})^{T}$ and $\ell=(\nu_{L},~{}l_{L})^{T}$ stand for the
quark and lepton doublets of $SU(2)_{L}$ respectively, and we denote the
charge conjugation fields $e^{c}_{R}=(e_{R})^{c}$ and
$\ell^{c}=((\nu_{L})^{c},~{}(l_{L})^{c})^{T}$. The Roman letters $i,j,m,n=1,2$
refer to $SU(2)_{L}$ indices, and $\epsilon$ is the 2-component antisymmetric
tensor with $\epsilon_{12}=-\epsilon_{21}=1$, and $\lambda^{A}$ are the Gell-
Mann matrices. Here we only consider operators involving first generation
fields, since we are concerned with the operators which directly contribute to
the $0\nu\beta\beta$ decay. In addition to $0\nu\beta\beta$ decays, these
operators in Eq. (3) can produce the signals of same sign dilepton, single
lepton $+$ MET or 2jets+MET at hadron colliders [13]. The different possible
ways to generate the above dimension nine ($d=9$) $0\nu\beta\beta$ decay
operators at tree level have been studied in Ref. [16]. All the tree-level
decomposition of the long-range $0\nu\beta\beta$ decay operators has been
discussed in Ref. [17]. The one-loop contribution to the $0\nu\beta\beta$
amplitude was studied in the radiative neutrino mass models [18, 19, 20, 21].
In this paper, we shall study the systematical decomposition of the eleven
operators $\mathcal{O}^{0\nu}_{1,2,\ldots,11}$ at one-loop level by using the
diagrammatic approach [22, 23].

In this work, we consider the cases that the mediators are either scalar or
fermion fields. The topologies of the diagrams for scalar and vector are the
same, nevertheless vector fields are expected to be the gauge bosons of some
new symmetry and their masses are generated through the spontaneous breaking
of that symmetry. That is to say the SM gauge group should be extended to
accommodate the new vector bosons, consequently we only discuss scalar
messengers here. Furthermore, we assume the new fermions are vector-like so
that the resulting models are anomaly free.

This paper is organized as follows. In section 2, the general procedures of
decomposing a non-renormalizable operator are sketched, and we discuss how to
generate the topologies, diagrams and models of the short-range
$0\nu\beta\beta$ operators. In section 3, the relation between
$0\nu\beta\beta$ models and neutrino masses is discussed. In section 4, we
give an example of one-loop model for short-range $0\nu\beta\beta$ decay, the
leading contribution to neutrino masses arises at two-loop level, and the
prediction for half-life is analyzed. Section 5 summarizes our main results.
The complete set of $0\nu\beta\beta$ decay operators below the electroweak
scale as well as the general formula for half-life are given in Appendix A. We
list the non-renormalizable topologies and finite non-genuine diagrams in
Appendix B. We present the possible one-loop models and the quantum number
assignments of mediators for the $0\nu\beta\beta$ decay operator N4 with the
diagram N-1-1-1 in Appendix C.

##  2 Generate genuine topologies, diagrams and models

The diagrammatic decomposition approach [22, 23] allows us to systematically
identify the possible ultraviolet (UV) completion of any effective high
dimensional operator $\mathcal{O}_{EF}$ in the SM effective theory, and it has
been used to systematically classify the Majorana neutrino mass models
generated at one-loop level [24], two-loop level [25] and three-loop level
[26]. The classification of Dirac neutrino mass models have also been
performed in this approach [27, 28, 29, 30]. The general procedure of
decomposing a non-renormalizable operator $\mathcal{O}_{EF}$ at $L-$ loop
level is as follows: ① Generate the $L-$loop connected topology with three-
point and four-point vertices, and the number of the external lines should be
exactly equal to the number of the fields contained in the considered
operators. The topologies with tadpoles and self-energies are generally
discarded. ② Assign the fields involved in the concerned operators to the
external lines, and specify the Lorentz nature (spinor or scalar) of each
internal line. Thus the topologies are promoted to diagrams. For the four-
point vertex, the corresponding interaction is renormalizable if and only if
all the four relevant legs are scalars, while the trilinear vertices can be
the fermion-fermion-scalar and scalar-scalar-scalar interactions. We discard
all the diagrams with non-renormalizable interactions. ③ Construct the UV
complete models by specifying the quantum numbers of the internal fields, and
each vertex should be invariant under the SM gauge symmetry.

Although the resulting UV models obtained in this way can generate the
concerned effective operator $\mathcal{O}_{EF}$ at $L-$loop level, one should
guarantee that it is the leading order contribution. It is called a genuine
model if there is no lower order contribution to $\mathcal{O}_{EF}$ otherwise
the model is non-genuine. The topologies associated with the genuine models is
called genuine topologies. If all the UV models generated from a topology are
non-genuine, the corresponding topology is non-genuine. In a similar fashion,
one can define genuine diagrams and non-genuine diagrams. A diagram is non-
genuine if it is possible to compress one or more loops into a renormalizable
vertex [26], in this case, lower order contributions will be allowed. In
general, any internal loop (or loops) with three external legs can be
compressed into a three point vertex regardless of the field content of a
model, thus the associated topology is non-genuine. It is only possible to
compress a internal loop (or loops) with four external legs into a
renormalizable vertex if all external lines are scalars. In the present work,
we shall apply the above procedures to study the one-loop decomposition of the
neutrinoless double decay operators in Eq. (3).

###  2.1  Generate topologies

All the dimension 9 effective $0\nu\beta\beta$ decay operators involve six
quark and lepton fields. We first use our own code to generate the connected
one-loop topologies with six external legs. All the topologies with tadpoles
and self-energy are excluded, because their loop integrals are divergent and
thus lower order counter term is required in a consistent renormalization
scheme. After excluding the topologies which contain an internal loop with
three legs, we find there are totally 43 topologies shown in figure 1 and
figure 14. Only the six one-loop topologies in figure 1 can lead to genuine
diagrams (models). The topologies in figure 14 need non-renormalizable
interactions to accommodate six external fermion lines. For the sake of
completeness, we show the non-genuine topologies which contain an internal
loop with three legs in figure 15. We see that all these non-genuine
topologies can be regarded as extensions of the tree-level $0\nu\beta\beta$
decay models, where one of the vertices is generated at one loop. As regards
the tree-level decomposition of the $0\nu\beta\beta$ decay operators, there
are only two renormalizable topologies [16] which are displayed in figure 1.
Generally three mediators are needed in the tree-level UV complete models and
their quantum numbers are unambiguously fixed. We also provide the tree-level
decomposition in the present work, since these results are essential when
determining the genuineness of a one-loop $0\nu\beta\beta$ decay model.

Figure 1: Genuine topologies which can lead to tree-level and one-loop genuine
models for the short-range part of the $0\nu\beta\beta$ decay amplitude.

###  2.2  Generate diagrams

All the six external legs for $0\nu\beta\beta$ decay are quark and lepton
fields. We proceed to insert fermion or scalar into each internal line, and
Lorentz invariance requires that each vertex must contain even number of
fermion. The UV complete models are assumed to be renormalizable, thus the
interaction vertices should be dimension-four three-point and four-point
vertices, and consequently the interaction terms in the Lagrangian are of the
form $\overline{F}_{1}F_{2}S_{1}$, $S_{1}S_{2}S_{3}$ or
$S_{1}S_{2}S_{3}S_{4}$, where $F_{i}$ and $S_{i}$ denote a generic fermion and
scalar fields respectively. A topology is non-renormalizable if the diagrams
for all possible fermionic and scalar insertions contain non-renormalizable
vertex. We show the non-renormalizable topologies in figure 14 for
completeness. Some diagrams arising from a topology are possibly superfluous.
One needs to consider all possible permutations of fermion and scalar lines
and compare the couplings of the vertices. If the couplings of two diagrams
match with each other, they are exactly the same diagrams. For the genuine
topologies in figure 1, we show the corresponding independent diagrams in
figure 2. Both topologies N-1-4 and N-1-5 give rise to two diagrams, while
there is only one diagram for each of the remaining four topologies N-1-1,
N-1-2, N-1-3 and N-1-6. Because all the external legs of the short-range
$0\nu\beta\beta$ decay are fermions, the four-point interaction
$S_{1}S_{2}S_{3}S_{4}$ is only involved in the diagram N-1-1-1 and the
interaction vertices are of the forms $\overline{F}_{1}F_{2}S_{1}$ and
$S_{1}S_{2}S_{3}$ in other cases, as can be seen from figure 2. In particular,
it is impossible to compress the loop into a renormalizable vertex, hence we
expect all these diagrams in figure 2 could give genuine models. We also show
the tree-level diagrams which lead to $0\nu\beta\beta$ decay in figure 2,
since they are essential when determining whether a loop model is genuine.
There are only two diagrams at tree level, and the Lorentz natures of all
mediators are completely fixed.

Figure 2: List of genuine diagrams arising from the topologies in figure 1,
where the solid and dashed lines stand for fermion and scalar fields
respectively.

###  2.3 Generate models

The next step is to generate models for the diagrams displayed in figure 2,
and we need to specify the $SU(3)_{C}$, $SU(2)_{L}$ and $U(1)_{Y}$ quantum
numbers to each line in these diagrams. Firstly we attach the fields in the
$0\nu\beta\beta$ decay operators to external lines of the diagram,
subsequently gauge invariance allows us to determine the quantum numbers of
each internal line. Notice that there are only finite possibilities for the
quantum numbers of the internal fermions and scalars in the tree-level models.
However, any one-loop diagram in figure 2 leads to an infinite number of
genuine models. For illustration, we shall impose an upper limit on the
dimensions of the $SU(3)_{C}$ and $SU(2)_{L}$ representations assigned to the
mediator fields. The details of how to generate models are given in the
following subsections.

####  2.3.1 Attach external fields

The eleven $0\nu\beta\beta$ decay operators in Eq. (3) can be divided into six
categories as listed in table 1, and the operators of each category involve
the same fields as well as the same UV completion. The vector fields should be
gauge bosons of certain gauge symmetry and their masses are generated by the
spontaneous breaking of symmetry. Thus it would be necessary to extend the SM
gauge group if vectors are mediators. Therefore we focus on the scalar and
fermion mediated cases in the present work, the cases of vectors as mediators
would be not considered. The SM lepton and quark fields are chiral fields in
weak basis, Lorentz invariance requires vector mediators for certain
attachment of external fields, and these cases would be dropped. As an
example, we shall consider the one-loop UV completion of the $0\nu\beta\beta$
decay operator N4 based on the topology N-1-1. Let us consider the diagram in
figure 3, there are six possible ways to attach the fields of the operator N4
to the external lines of the vertex A. If the two external legs are assigned
to be $u_{R}$, $d_{R}$ or $e_{R}$, $d_{R}$, the non-vanishing fermion
bilinears are $\bar{u}_{R}\gamma^{\mu}d_{R}$ and
$\bar{e}_{R}\gamma^{\mu}d_{R}$ respectively, consequently the internal leg has
to be a vector field $X_{\mu}$ to form a Lorentz invariant interaction
$\bar{u}_{R}\gamma^{\mu}d_{R}X_{\mu}$ and
$\bar{e}_{R}\gamma^{\mu}d_{R}X_{\mu}$. As a result, these two cases will be
discarded. Since we can freely attach the fields to the external lines,
therefore there are generally multiple external line structures for a given
diagram and some of them are superfluous. The redundant assignments of
external lines must base on the same diagram. In order to find out the
redundant ones, we need to consider the possible permutations of vertices and
compare the interlinkage between vertices. If several kinds of attachments
have the same structure of interlinkage, they will be identified as one.

Notation |  $0\nu\beta\beta$ decay operators |  External fields
---|---|---
N1 | $\mathcal{O}^{0\nu}_{1},~{}\mathcal{O}^{0\nu}_{2}$ | $\overline{Q},Q,\bar{u}_{R},d_{R},\overline{\ell},\ell^{c}$
N2 | $\mathcal{O}^{0\nu}_{3},~{}\mathcal{O}^{0\nu}_{4}$ | $Q,Q,\bar{u}_{R},\bar{u}_{R},\overline{\ell},\ell^{c}$
N3 | $\mathcal{O}^{0\nu}_{5},~{}\mathcal{O}^{0\nu}_{6}$ | $\overline{Q},\overline{Q},d_{R},d_{R},\overline{\ell},\ell^{c}$
N4 | $\mathcal{O}^{0\nu}_{7}$ | $\bar{u}_{R},\bar{u}_{R},d_{R},d_{R},\bar{e}_{R},e_{R}^{c}$
N5 | $\mathcal{O}^{0\nu}_{8},~{}\mathcal{O}^{0\nu}_{9}$ | $\overline{u}_{R},\overline{Q},d_{R},d_{R},\overline{\ell},e_{R}^{c}$
N6 | $\mathcal{O}^{0\nu}_{10},~{}\mathcal{O}^{0\nu}_{11}$ | $\bar{u}_{R},\bar{u}_{R},Q,d_{R},\overline{\ell},e_{R}^{c}$
Table 1: Summary of fields involved in the effective $0\nu\beta\beta$ decay
operators of dimension 9. Figure 3:  The possible ways of attaching the fields
of the operator N4 to the external legs of the vertex A in diagram N-1-1-1,
where the dashed lines and wavy lines on the right side denote scalar and
vector fields respectively. Notice that Lorentz invariance fixes the mediator
to be scalar or vector because the external legs are fermions with definite
chirality.  Figure 4:  The possible ways to attach the fields of the operator
N4 to the diagram N-1-1-1.

####  2.3.2 Determine $U(1)_{Y}$ charge

Once the external legs are specified, the hypercharges fulfilling of the
internal lines can be determined by imposing the hypercharge conservation
vertex by vertex. Note that the hypercharge of a field is related to its
electric charge via the Gell-Mann-Nishijima formula $Q=T_{3}+Y$, where $T_{3}$
is the third component of the weak isospin. For example, the equations of
hypercharge conservation for the diagram N4-1-1-1-1 read as

$\displaystyle\qquad\qquad
Y_{F_{1}}-Y_{d_{R}}-Y_{S_{1}}=0,~{}~{}~{}~{}Y_{F_{1}}+Y_{S_{2}}+Y_{d_{R}}=0,$
$\displaystyle
Y_{S_{1}}+Y_{S_{2}}+Y_{S_{3}}+Y_{S_{4}}=0,~{}~{}~{}Y_{e_{R}}+Y_{u_{R}}-Y_{S_{1}}=0,~{}~{}~{}Y_{e_{R}}+Y_{u_{R}}-Y_{S_{4}}=0\,,$
(4)

with $Y_{e_{R}}=-1$, $Y_{u_{R}}=\frac{2}{3}$ and $Y_{d_{R}}=-\frac{1}{3}$. The
solution to the above equations is given by

$\displaystyle
Y_{F_{1}}=-\frac{1}{3}+\alpha,~{}~{}~{}Y_{S_{1}}=\alpha,~{}~{}Y_{S_{2}}=-\frac{2}{3}-\alpha,~{}~{}~{}Y_{S_{3}}=Y_{S_{4}}=-\frac{1}{3}\,,$
(5)

where $\alpha$ is arbitrary real parameter. A definite value of $\alpha$
should be taken in a concrete model. In a similar way, we can determine the
hypercharges for the other diagrams in figure 4 as follows,

$\displaystyle\text{N4-1-1-1-2}~{}:~{}Y_{F_{1}}=-\frac{2}{3}+\alpha,~{}~{}Y_{S_{1}}=\alpha,~{}~{}Y_{S_{2}}=-\frac{1}{3}-\alpha,~{}~{}Y_{S_{3}}=-\frac{1}{3},~{}~{}Y_{S_{4}}=\frac{2}{3}\,,$
$\displaystyle\text{N4-1-1-1-3}~{}:~{}Y_{F_{1}}=1+\alpha,~{}~{}Y_{S_{1}}=\alpha,~{}~{}Y_{S_{2}}=-2-\alpha,~{}~{}Y_{S_{3}}=\frac{4}{3},~{}~{}Y_{S_{4}}=\frac{2}{3}\,,$
$\displaystyle\text{N4-1-1-1-4}~{}:~{}Y_{F_{1}}=-\frac{1}{3}+\alpha,~{}~{}Y_{S_{1}}=\alpha,~{}~{}Y_{S_{2}}=\frac{2}{3}-\alpha,~{}~{}Y_{S_{3}}=-2,~{}~{}Y_{S_{4}}=\frac{4}{3}\,,$
$\displaystyle\text{N4-1-1-1-5}~{}:~{}Y_{F_{1}}=-\frac{2}{3}+\alpha,~{}~{}Y_{S_{1}}=\alpha,~{}~{}Y_{S_{2}}=\frac{4}{3}-\alpha,~{}~{}Y_{S_{3}}=-2,~{}~{}Y_{S_{4}}=\frac{2}{3}\,.$
(6)

We see that a hypercharge flow encoded in the parameter $\alpha$ in the loop
is allowed.

####  2.3.3 Assign $SU(2)_{L}$ transformation

We proceed to assign the $SU(2)_{L}$ transformations to the internal lines,
and they are denoted by the dimensions of the $SU(2)_{L}$ representations. For
the concerned renormalizable interactions, the invariance under $SU(2)_{L}$
gauge symmetry requires the following conditions should be fulfilled,

$\displaystyle\overline{F}_{1}F_{2}S_{1}~{}$ $\displaystyle:$
$\displaystyle~{}n_{F_{1}}\otimes n_{F_{2}}\otimes n_{S_{1}}\supset 1\,,$
$\displaystyle S_{1}S_{2}S_{3}~{}$ $\displaystyle:$
$\displaystyle~{}n_{S_{1}}\otimes n_{S_{2}}\otimes n_{S_{3}}\supset 1\,,$
$\displaystyle S_{1}S_{2}S_{3}S_{4}~{}$ $\displaystyle:$
$\displaystyle~{}n_{S_{1}}\otimes n_{S_{2}}\otimes n_{S_{3}}\otimes
n_{S_{4}}\supset 1\,,$ (7)

where $n_{F_{i}}$ and $n_{S_{i}}$ denote the $SU(2)_{L}$ representations to
which the fields $F_{i}$ and $S_{i}$ are assigned. Notice that each $SU(2)$
representation and its complex conjugate are equivalent. Moreover, it is
necessary to check that some $0\nu\beta\beta$ decay operators can really be
generated after integrating out the mediators. For example, let us consider
the models in figure 5, the external legs are
$Q,Q,\bar{u}_{R},\bar{u}_{R},\overline{\ell},\ell^{c}$ which coincide with
those of the operators $\mathcal{O}^{0\nu}_{3}$ and $\mathcal{O}^{0\nu}_{4}$.
In the left panel, the renormalizable interaction $\bar{\ell}\ell^{c}S$ at the
vertex $A$ is vanishing because the $SU(2)$ singlet contraction of two
doublets is antisymmetric. The model in the right panel of figure 5 leads to
the effective operator
$\bar{e}_{L}\nu^{c}_{L}\bar{u}_{R}u_{L}\bar{u}_{R}d_{L}$ which can not induce
the $0\nu\beta\beta$ decay process, although its interaction vertices are
invariant under SM gauge symmetry. As a consequence, we need to explicitly
expand the $SU(2)_{L}$ index of the interaction at each vertex, contract the
internal lines and see whether the $0\nu\beta\beta$ decay operator can be
produced. For illustration, we shall consider the $SU(2)_{L}$ singlet, doublet
and triplet assignments for the messenger fields in the present work, and
higher dimensional representations can be discussed in the same fashion.

Figure 5:  Both models shown in panels (a) and (b) fulfill the $SU(2)_{L}$
invariance, where the numbers 1 and 2 denote the $SU(2)$ transformations of
the fields. For the diagram (a), there is a vanishing point interaction
$\bar{\ell}\ell^{c}S$ at the vertex $A$. The model in panel (b) doesn’t
contribute to the $0\nu\beta\beta$ decay.

The SM quantum number assignments for the internal fields are independent of
their Lorentz property, therefore this step can be performed at the topology
level. In table 2, we list the possible $SU(2)_{L}$ assignments of fields for
the topology N-1-1. For the set of operators N1, N2 and N3, four of the
external fields are $SU(2)_{L}$ doublets and the other two are singlets,
consequently there are four different assignments of external legs, as shown
in table 2. The external lines of the operator N4 are all singlets under
$SU(2)_{L}$. The fields of the operators N5 and N6 include two doublets and
four singlets of $SU(2)$, and there are four different assignments of external
fields as well. The possible $SU(2)_{L}$ quantum numbers of the internal
fields can be fixed by the invariance condition in Eq. (7), and it is
necessary to check that the $0\nu\beta\beta$ decay operators can really be
produced, as explained above.

Figure 6: Field assignments and color flow for the topology N-1-1, where $E_{i}$ and $I_{i}$ denote the external and internal fields respectively.  Operator | $E_{1}$ | $E_{2}$ | $E_{3}$ | $E_{4}$ | $E_{5}$ | $E_{6}$ | $I_{1}$ | $I_{2}$ | $I_{3}$ | $I_{4}$ | $I_{5}$ | Operator | $E_{1}$ | $E_{2}$ | $E_{3}$ | $E_{4}$ | $E_{5}$ | $E_{6}$ | $I_{1}$ | $I_{2}$ | $I_{3}$ | $I_{4}$ | $I_{5}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
N1,N2,N3 | $1$ | $2$ | $2$ | $2$ | $1$ | $2$ | $2$ | $1$ | $2$ | $2$ | $1$ | N1,N2,N3 | $2$ | $2$ | $1$ | $1$ | $2$ | $2$ | $2$ | $1$ | $1$ | $1$ | $1$
$2$ | $1$ | $2$ | $2$ | $3$ | $2$ | $1$ | $3$ | $3$ | $1$
$1$ | $2$ | $1$ | $2$ | $1$ | $1$ | $2$ | $2$ | $1$ | $1$
$1$ | $2$ | $1$ | $2$ | $3$ | $1$ | $2$ | $2$ | $3$ | $1$
$3$ | $2$ | $3$ | $2$ | $1$ | $3$ | $2$ | $2$ | $1$ | $1$
$3$ | $2$ | $3$ | $2$ | $3$ | $3$ | $2$ | $2$ | $3$ | $1$
$2$ | $3$ | $2$ | $2$ | $1$ | $2$ | $3$ | $3$ | $1$ | $1$
$2$ | $3$ | $2$ | $2$ | $3$ | $2$ | $3$ | $3$ | $3$ | $1$
$1$ | $1$ | $2$ | $2$ | $2$ | $2$ | $1$ | $1$ | $1$ | $1$ | $1$ | $2$ | $2$ | $1$ | $2$ | $2$ | $1$ | $2$ | $1$ | $1$ | $2$ | $2$
$1$ | $1$ | $1$ | $3$ | $3$ | $2$ | $1$ | $3$ | $2$ | $2$
$2$ | $2$ | $2$ | $1$ | $1$ | $1$ | $2$ | $2$ | $2$ | $2$
$2$ | $2$ | $2$ | $3$ | $1$ | $3$ | $2$ | $2$ | $2$ | $2$
$2$ | $2$ | $2$ | $3$ | $3$ | $2$ | $3$ | $3$ | $2$ | $2$
$3$ | $3$ | $3$ | $1$ | $1$ | N4 | $1$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$
$3$ | $3$ | $3$ | $3$ | $1$ | $2$ | $2$ | $2$ | $1$ | $1$
$3$ | $3$ | $3$ | $3$ | $3$ | $3$ | $3$ | $3$ | $1$ | $1$
N5,N6 | $1$ | $1$ | $1$ | $1$ | $2$ | $2$ | $1$ | $1$ | $1$ | $1$ | $1$ | N5,N6 | $2$ | $1$ | $1$ | $1$ | $2$ | $1$ | $2$ | $2$ | $1$ | $2$ | $1$
$2$ | $2$ | $2$ | $1$ | $1$ | $1$ | $1$ | $2$ | $2$ | $1$
$2$ | $2$ | $2$ | $3$ | $1$ | $3$ | $3$ | $2$ | $2$ | $1$
$3$ | $3$ | $3$ | $1$ | $1$ | $2$ | $2$ | $3$ | $2$ | $1$
$3$ | $3$ | $3$ | $3$ | $1$ | $2$ | $2$ | $1$ | $1$ | $1$ | $1$ | $2$ | $1$ | $1$ | $1$ | $1$
$1$ | $1$ | $1$ | $2$ | $2$ | $1$ | $1$ | $1$ | $1$ | $2$ | $2$ | $1$ | $2$ | $2$ | $1$ | $1$
$2$ | $2$ | $2$ | $2$ | $2$ | $3$ | $2$ | $2$ | $1$ | $1$
$3$ | $3$ | $3$ | $2$ | $2$ | $2$ | $3$ | $3$ | $1$ | $1$
Table 2: Assign the $SU(2)_{L}$ quantum numbers to fields of the topology
N-1-1 in figure 1. The labels $E_{i}$ and $I_{i}$ of the external and internal
lines are showed in figure 6.

####  2.3.4 Assign $SU(3)_{C}$ transformation

$E_{1}$ | $E_{2}$ | $E_{3}$ | $E_{4}$ | $E_{5}$ | $E_{6}$ | $I_{1}$ | $I_{2}$ | $I_{3}$ | $I_{4}$ | $I_{5}$ | $E_{1}$ | $E_{2}$ | $E_{3}$ | $E_{4}$ | $E_{5}$ | $E_{6}$ | $I_{1}$ | $I_{2}$ | $I_{3}$ | $I_{4}$ | $I_{5}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$1$ | $\bar{3}$ | $3$ | $\bar{3}$ | $1$ | $3$ | $1$ | $3$ | $1$ | $\bar{3}$ | $1$ | $1$ | $3$ | $\bar{3}$ | $3$ | $\bar{3}$ | $1$ | $1$ | $\bar{3}$ | $1$ | $3$ | $1$
$3$ | $1$ | $3$ | $\bar{3}$ | $1$ | $3$ | $3$ | $3$ | $3$ | $1$
$\bar{3}$ | $\bar{3}$ | $\bar{3}$ | $\bar{3}$ | $1$ | $3$ | $\bar{6}$ | $3$ | $3$ | $1$
$\bar{3}$ | $6$ | $\bar{3}$ | $\bar{3}$ | $1$ | $\bar{3}$ | $1$ | $\bar{3}$ | $3$ | $1$
$6$ | $\bar{3}$ | $6$ | $\bar{3}$ | $1$ | $\bar{6}$ | $3$ | $\bar{6}$ | $3$ | $1$
$\bar{3}$ | $\bar{3}$ | $3$ | $1$ | $1$ | $3$ | $\bar{3}$ | $\bar{3}$ | $1$ | $\bar{3}$ | $\bar{3}$ | $1$ | $1$ | $3$ | $\bar{3}$ | $\bar{3}$ | $3$ | $1$ | $1$ | $1$ | $1$ | $1$
$\bar{3}$ | $6$ | $1$ | $\bar{3}$ | $\bar{3}$ | $3$ | $\bar{3}$ | $3$ | $1$ | $1$
$1$ | $3$ | $3$ | $\bar{3}$ | $\bar{3}$ | $6$ | $\bar{6}$ | $6$ | $1$ | $1$
$1$ | $3$ | $1$ | $3$ | $\bar{3}$ | $\bar{3}$ | $1$ | $\bar{3}$ | $1$ | $\bar{3}$ | $\bar{3}$ | $1$ | $\bar{3}$ | $3$ | $3$ | $\bar{3}$ | $1$ | $1$ | $3$ | $1$ | $3$ | $3$
$1$ | $\bar{3}$ | $1$ | $6$ | $\bar{3}$ | $1$ | $3$ | $1$ | $3$ | $\bar{6}$
$3$ | $3$ | $3$ | $\bar{3}$ | $\bar{3}$ | $3$ | $1$ | $3$ | $3$ | $3$
$3$ | $3$ | $3$ | $6$ | $\bar{3}$ | $3$ | $1$ | $3$ | $3$ | $\bar{6}$
$3$ | $\bar{6}$ | $3$ | $\bar{3}$ | $\bar{3}$ | $\bar{3}$ | $\bar{3}$ | $\bar{3}$ | $3$ | $3$
$3$ | $\bar{6}$ | $3$ | $6$ | $\bar{3}$ | $\bar{3}$ | $\bar{3}$ | $\bar{3}$ | $3$ | $\bar{6}$
$\bar{3}$ | $1$ | $\bar{3}$ | $\bar{3}$ | $\bar{3}$ | $\bar{3}$ | $6$ | $\bar{3}$ | $3$ | $3$
$\bar{3}$ | $1$ | $\bar{3}$ | $6$ | $\bar{3}$ | $\bar{3}$ | $6$ | $\bar{3}$ | $3$ | $\bar{6}$
$\bar{6}$ | $3$ | $\bar{6}$ | $\bar{3}$ | $\bar{3}$ | $6$ | $\bar{3}$ | $6$ | $3$ | $3$
$\bar{6}$ | $3$ | $\bar{6}$ | $6$ | $\bar{3}$ | $6$ | $\bar{3}$ | $6$ | $3$ | $\bar{6}$
$\bar{3}$ | $3$ | $1$ | $3$ | $\bar{3}$ | $1$ | $\bar{3}$ | $1$ | $1$ | $3$ | $\bar{3}$ | $3$ | $3$ | $1$ | $\bar{3}$ | $\bar{3}$ | $1$ | $\bar{3}$ | $1$ | $3$ | $3$ | $3$
$1$ | $\bar{3}$ | $3$ | $3$ | $\bar{3}$ | $\bar{3}$ | $1$ | $\bar{6}$ | $3$ | $3$
$3$ | $3$ | $\bar{3}$ | $3$ | $\bar{3}$ | $1$ | $\bar{3}$ | $\bar{3}$ | $3$ | $3$
$3$ | $\bar{6}$ | $\bar{3}$ | $3$ | $\bar{3}$ | $\bar{3}$ | $\bar{3}$ | $3$ | $3$ | $1$ | $1$ | $\bar{3}$ | $\bar{3}$ | $1$ | $1$ | $3$
$\bar{6}$ | $3$ | $\bar{3}$ | $3$ | $\bar{3}$ | $\bar{3}$ | $6$ | $1$ | $1$ | $\bar{6}$
$3$ | $3$ | $6$ | $3$ | $\bar{3}$ | $1$ | $3$ | $3$ | $1$ | $3$
$3$ | $\bar{6}$ | $6$ | $3$ | $\bar{3}$ | $1$ | $3$ | $3$ | $1$ | $\bar{6}$
$1$ | $1$ | $3$ | $3$ | $\bar{3}$ | $\bar{3}$ | $1$ | $1$ | $1$ | $\bar{3}$ | $3$ | $\bar{3}$ | $3$ | $\bar{3}$ | $3$ | $1$ | $1$ | $\bar{3}$ | $1$ | $1$ | $1$ | $1$
$1$ | $1$ | $1$ | $6$ | $\bar{6}$ | $1$ | $\bar{3}$ | $3$ | $1$ | $1$
$3$ | $\bar{3}$ | $3$ | $\bar{3}$ | $3$ | $3$ | $3$ | $\bar{3}$ | $1$ | $1$
$3$ | $\bar{3}$ | $3$ | $\bar{3}$ | $\bar{6}$ | $\bar{6}$ | $3$ | $\bar{3}$ | $1$ | $1$
$3$ | $\bar{3}$ | $3$ | $6$ | $3$ | $3$ | $\bar{6}$ | $6$ | $1$ | $1$
$3$ | $\bar{3}$ | $3$ | $6$ | $\bar{6}$ | $3$ | $3$ | $\bar{3}$ | $\bar{3}$ | $1$ | $1$ | $\bar{3}$ | $1$ | $3$ | $1$ | $\bar{3}$
$6$ | $\bar{6}$ | $6$ | $\bar{3}$ | $3$ | $\bar{3}$ | $1$ | $\bar{6}$ | $1$ | $6$
$6$ | $\bar{6}$ | $6$ | $\bar{3}$ | $\bar{6}$ | $1$ | $\bar{3}$ | $\bar{3}$ | $1$ | $\bar{3}$
$6$ | $\bar{6}$ | $6$ | $6$ | $3$ | $1$ | $\bar{3}$ | $\bar{3}$ | $1$ | $6$
$6$ | $\bar{6}$ | $6$ | $6$ | $\bar{6}$ |  |  |  |  |  |  |  |  |  |  |
Table 3:  Assign the $SU(3)_{C}$ quantum numbers to fields of the topology
N-1-1 in figure 1. The labels $E_{i}$ and $I_{i}$ of the external and internal
fields are showed in figure 6.

Then we discuss the possible transformation properties of the internal fields
under $SU(3)_{C}$, we consider the $SU(3)$ representation up to dimension 8,
namely, $1$, 3, $\bar{3}$, 6, $\bar{6}$ and $8$, high dimensional
representations can be discussed analogously. All the $0\nu\beta\beta$ decay
operators in Eq. (3) involve two quarks, two anti-quarks and two lepton fields
which are in the irreducible representations 3, $\bar{3}$ and $1$ of the
$SU(3)_{C}$ group respectively. We assign $SU(3)_{C}$ quantum numbers to each
internal line, then consider whether $SU(3)_{C}$ invariant contractions can be
formed at the vertices by using the tensor product reduction rules such as

$\displaystyle 3\otimes 3=\bar{3}\oplus 6,~{}~{}~{}3\otimes\bar{3}=1\oplus
8,~{}~{}~{}3\otimes 6=8\otimes
10,~{}~{}~{}3\otimes\bar{6}=\bar{3}\oplus\overline{15}\,,$ $\displaystyle
3\otimes 8=3\oplus\bar{6}\oplus 15,~{}~{}~{}3\otimes 10=15\oplus
15,~{}~{}~{}3\otimes\overline{10}=\overline{6}\oplus\overline{24}\,,$
$\displaystyle 6\otimes 6=\bar{6}\oplus 15\oplus
15,~{}~{}~{}6\otimes\bar{6}=1\oplus 8\oplus 27\,,~{}~{}~{}6\otimes
8=\bar{3}\oplus 6\oplus\overline{15}\oplus 24,$ $\displaystyle 6\otimes
10=\overline{15}\oplus 21\oplus
24,~{}~{}~{}6\otimes\overline{10}=\bar{3}\oplus\overline{15}\oplus\overline{42}\,,~{}~{}~{}8\otimes
8=1\oplus 8\oplus 8\oplus 10\oplus\overline{10}\oplus 27\,,$ $\displaystyle
8\otimes 10=8\oplus 10\oplus 27\oplus 35,~{}~{}~{}10\otimes
10=\overline{10}\oplus 27\oplus 28\oplus
35,~{}~{}~{}10\otimes\overline{10}=1\oplus 8\oplus 27\oplus 64\,.$ (8)

We shall use the Mathematica package GroupMath [31] to determine the
$SU(3)_{C}$ gauge invariance of the interaction vertices. The renormalizable
interaction compatible with $SU(3)_{C}$ gauge symmetry could be vanishing, in
particular, when the repeated fields are involved. We take the interaction
displayed in figure 7 as an example, which is of the form
$\overline{d_{R}^{c}}d_{R}S$, with $S$ being a color triplet. Expanding the
color index of the down quark explicitly, we obtain

$\displaystyle\overline{d_{R}^{c}}d_{R}S$ $\displaystyle=$
$\displaystyle\overline{d_{Rr}^{c}}d_{Rg}S_{b}-\overline{d_{Rg}^{c}}d_{Rr}S_{b}-\overline{d_{Rr}^{c}}d_{Rb}S_{g}+\overline{d_{Rb}^{c}}d_{Rr}S_{g}+\overline{d_{Rg}^{c}}d_{Rb}S_{r}-\overline{d_{Rb}^{c}}d_{Rg}S_{r}$
(9) $\displaystyle=$
$\displaystyle\overline{d_{Rr}^{c}}d_{Rg}S_{b}-\overline{d_{Rr}^{c}}d_{Rg}S_{b}-\overline{d_{Rr}^{c}}d_{Rb}S_{g}+\overline{d_{Rr}^{c}}d_{Rb}S_{g}+\overline{d_{Rg}^{c}}d_{Rb}S_{r}-\overline{d_{Rg}^{c}}d_{Rb}S_{r}$
$\displaystyle=$ $\displaystyle 0\,,$

where $r$, $g$, $b$ denote the three different colors of quark fields. Hence
it is necessary to explicitly expand the $SU(3)_{C}$ indices of each
interaction vertex to see whether point interactions are vanishing. The
$SU(3)_{C}$ quantum numbers of the internal fields can be discussed at the
topology level as well. We list the possible color charge assignments of the
topology N-1-1 in table 3, where the $SU(3)_{C}$ representations up to
dimension 6 are considered due to limited space and the results including
octet assignment are given in the attachment [32]. The labels of the fields
and the color flow of the topology N-1-1 are shown in figure 6.

As explained above, we can start from a diagram in figure 2, attach the fields
of the $0\nu\beta\beta$ decay operators to the external legs, subsequently
determine the possible $U(1)_{Y}$, $SU(2)_{L}$ and $SU(3)_{C}$ quantum numbers
of the internal fields by imposing gauge invariance at each vertex, then the
one-loop UV complete models for the effective $0\nu\beta\beta$ decay operators
can be obtained. The different models are named as Ni-a-b-c-d-e-f, where Ni
with i=1, 2, 3, 4, 5, 6 denotes the set of $0\nu\beta\beta$ decay operators,
the second index “a” refer to tree (a=0) and one-loop (a=1) models, the
symbols “b” and “c” stand for the topology and diagram respectively, the index
“d” indicates the attachment of external fields, and the last two indices “e”
and “f” denote the assignments of the $SU(2)_{L}$ and $SU(3)_{C}$
respectively. For instance, combining figure 4 with the SM electroweak-sector
quantum number assignments in table 2 and table 3, we can obtain all the
models arising from the diagram N-1-1-1. The possible UV completions for the
operator N4 are listed in table 6, table 7, table 8, table 9 and table 10.

Figure 7:  An example of $SU(3)_{C}$ invariant vertex which is vanishing
because of antisymmetric color contractions.

####  2.3.5 Determine genuineness

Numerous one-loop models for short-range $0\nu\beta\beta$ decays can be
generated through the algorithm described in previous sections, in particular
all the tree-level $0\nu\beta\beta$ decay models are generated. A one-loop
model is the dominant contribution if and only if the new messenger fields are
insufficient to mediate the more important tree-level diagram, then such model
would be called genuine. In other words, the fields which generate the tree-
level $0\nu\beta\beta$ decays should be absent in a genuine one-loop model. If
the lower order contributions can not be forbidden without extra symmetry, the
model would be non-genuine. We determine the genuineness of each model by
comparing its field content with that of tree-level models one by one. Since
the quantum numbers of the internal particles in the tree-level
$0\nu\beta\beta$ decay models are unambiguously fixed, genuineness of a one-
loop model generally excludes certain value of the hypercharge parameter
$\alpha$. We take the model N1-1-1-1-2-1-1 for illustration, see figure 8, we
have defined the notation $C_{L}^{Y}$ for the quantum numbers in the diagram,
where $C$ refers to the $SU(3)_{C}$ transformation, $L$ refers to the
$SU(2)_{L}$ transformation, and $Y$ stands for $U(1)_{Y}$ charge. When the
hypercharge parameter $\alpha=\pm 1$, the mediators $S_{2}$ and $S_{3}$ and
the associated interactions allow to generate the tree-level $0\nu\beta\beta$
decay models N1-0-1-1-2-1-1, N2-0-1-1-1-1-1, N3-0-1-1-3-1-1 which are the more
important ones, as shown in figure 8. Therefore the genuineness of the one-
loop model N1-1-1-1-2-1-1 requires $\alpha\neq\pm 1$. The condition of
genuineness has been considered for each possible one-loop decomposition of
the $0\nu\beta\beta$ decay operators, see the attachment for the full results
[32].

Figure 8:  One example of a non-genuine model N1-1-1-1-2-1-1, the mediators
$S_{2}$ and $S_{3}$ will also lead to tree-level contributions N1-0-1-1-2-1-1,
N2-0-1-1-1-1-1, N3-0-1-1-3-1-1 when the hypercharge parameter $\alpha=\pm 1$.
The quantum numbers are given in the notation $C_{L}^{Y}$, where $C$ refers to
the $SU(3)_{C}$ transformation, $L$ refers to the $SU(2)_{L}$ transformation,
and $Y$ stands for $U(1)_{Y}$ charge.

##  3 Neutrino masses in $0\nu\beta\beta$ decay models

We proceed to discuss the neutrino masses which are generated in the above
models of short-range $0\nu\beta\beta$ decay. The famous black box theorem
establishes that the models which can lead to $0\nu\beta\beta$ decay must
produce Majorana neutrino masses at or below four-loop order [2]. The four-
loop diagram guaranteed by the black box theorem produces a very tiny neutrino
mass of order $m_{\nu}\sim 10^{-24}$ eV which is too small to explain the
current neutrino experimental data [33]. However, generally the neutrino
masses are generated at a lower loop order in a specific $0\nu\beta\beta$
decay model so that the correct order of neutrino masses could be obtained. It
is well known that the exchange of a light neutrino between two SM charged
current interactions could lead to $0\nu\beta\beta$ decay, and it is usually
called mass mechanism which is characterized by the effective neutrino mass
$m_{\beta\beta}$ in Eq. (2). For models with tree and one-loop neutrino
masses, the mass mechanism is generally expected to dominate unless the model
parameters are fine-tuned. The short-range contribution could be comparable
with the mass mechanism for models with two-loop and three-loop neutrino
masses [34]. One can construct neutrino mass diagrams by using the SM fields
and the mediators which appear in one-loop renormalizable $0\nu\beta\beta$
decay models. In order to identify the neutrino mass diagrams, one need to
compare the messenger fields of the $0\nu\beta\beta$ decay models with those
of neutrino mass models. In the present work, we are concerned with the
neutrino mass models up to two-loop level which can be generated from
$0\nu\beta\beta$ decay models.

Figure 9:  The neutrino mass models generated from the $0\nu\beta\beta$ decay
model N1-1-1-1-1-1-1. The one-loop neutrino mass diagram W-1-2-2-1-1-1 can be
generated by the mediators of N1-1-1-1-1-1-1 when the hypercharge
$\alpha=-7/6,5/6$, and the two-loop neutrino mass diagram W-2-1-2-1-11-7 can
always be generated.

We take model N1-1-1-1-1-1-1 for example. The new mediators could possibly
give rise to the one-loop neutrino mass diagram W-1-2-2-1-1-1, as shown in
figure 9. By examining the $SU(3)_{C}$ and $SU(2)_{L}$ quantum numbers of the
mediators, one could possibly identify the fields in the mass model
W-1-2-2-1-1-1 with those of the $0\nu\beta\beta$ decay model as follows,

$S_{1}^{\prime}\sim S_{4}^{(\dagger)}\,\quad S_{2}^{\prime}\sim
S_{2}^{(\dagger)},H\,\quad F_{1}^{\prime}\sim F^{(c)},\ell_{L}^{(c)}\,\quad
F_{2}^{\prime}\sim e_{R}^{(c)}\,.$ (10)

Furthermore, the identification of the hypercharge requires the following
relations should be fulfilled,

$\displaystyle
S_{1}^{\prime}:\alpha_{1}=-\frac{1}{6}-\alpha~{}~{}~{}\text{or}~{}~{}~{}\alpha_{1}=\frac{1}{6}+\alpha\,,$
$\displaystyle
S_{2}^{\prime}:\frac{1}{2}-\alpha_{1}=\frac{1}{2}~{}~{}~{}\text{or}~{}~{}~{}\frac{1}{2}-\alpha_{1}=-\frac{1}{2}\,,$
$\displaystyle
F_{1}^{\prime}:\frac{1}{2}-\alpha_{1}=-\frac{1}{3}+\alpha~{}~{}\text{or}~{}~{}\frac{1}{2}-\alpha_{1}=\frac{1}{3}-\alpha~{}~{}\text{or}~{}~{}\frac{1}{2}-\alpha_{1}=\frac{1}{2}~{}~{}\text{or}~{}~{}\frac{1}{2}-\alpha_{1}=-\frac{1}{2}\,,$
$\displaystyle
F_{2}^{\prime}:\alpha_{1}=-1~{}~{}~{}\text{or}~{}~{}~{}\alpha_{1}=1\,.$ (11)

The solutions to these equations are given by

$\left\\{\begin{array}[]{cl}\alpha&=-7/6,\\\
\alpha_{1}&=1,\end{array}\right.\quad\text{or}\quad\left\\{\begin{array}[]{cl}\alpha&=5/6,\\\
\alpha_{1}&=1.\end{array}\right.$ (12)

Hence the interactions in the $0\nu\beta\beta$ decay model N1-1-1-1-1-1-1 lead
us to the one-loop neutrino mass model W-1-2-2-1-1-1 when $\alpha=-7/6$ or
$\alpha=5/6$. Furthermore, by matching with the mediators of the one-loop
decomposition of the Weinberg operator [24], we find that the absence of one-
loop neutrino masses excludes the values of $\alpha=$$-7/6,-1/6,5/6$. Then we
proceed to show that the two-loop neutrino mass diagram W-2-1-2-1-11-7 given
in figure 9 can be constructed without introducing new interactions. All the
mediators of the model W-2-1-2-1-11-7 should be contained in the model
N1-1-1-1-1-1-1, their transformations under $SU(3)_{C}\times SU(2)_{L}$
implies the following possible identifications,

$S_{1}^{\prime\prime}\sim S_{1}^{\dagger}\,,~{}~{}S_{2}^{\prime\prime}\sim
S_{3}^{\dagger}\,,~{}~{}S_{3}^{\prime\prime}\sim
S_{4}^{(\dagger)}\,,~{}~{}F_{1}^{\prime\prime}\sim
Q_{L}\,,~{}~{}F_{2}^{\prime\prime}\sim
u_{R},d_{R}\,,~{}~{}F_{3}^{\prime\prime}\sim F^{(c)},\ell_{L}^{(c)}\,.$ (13)

The equality of hypercharges imposes the following constraints:

$\displaystyle S_{1}^{\prime\prime}:\beta_{2}=\frac{1}{3}\,,$ $\displaystyle
S_{2}^{\prime\prime}:\alpha_{2}=-\alpha\,,$ $\displaystyle
S_{3}^{\prime\prime}:\frac{1}{2}-\alpha_{2}-\beta_{2}=-\frac{1}{6}-\alpha\quad\text{or}\quad\frac{1}{2}-\alpha_{2}-\beta_{2}=\frac{1}{6}+\alpha\,,$
$\displaystyle F_{1}^{\prime\prime}:\frac{1}{2}-\beta_{2}=\frac{1}{6}\,,$
$\displaystyle
F_{2}^{\prime\prime}:-\beta_{2}=\frac{2}{3}\quad\text{or}\quad-\beta_{2}=-\frac{1}{3}\,,$
$\displaystyle
F_{3}^{\prime\prime}:\alpha_{2}+\beta_{2}=-\frac{1}{3}+\alpha~{}~{}\text{or}~{}~{}\alpha_{2}+\beta_{2}=\frac{1}{3}-\alpha~{}~{}\text{or}~{}~{}\alpha_{2}+\beta_{2}=-\frac{1}{2}~{}~{}\text{or}~{}~{}\alpha_{2}+\beta_{2}=\frac{1}{2}\,,$
(14)

whose solution is

$\alpha=-\alpha_{2},~{}~{}~{}\beta_{2}=1/3\,.$ (15)

This implies that the $0\nu\beta\beta$ decay model N1-1-1-1-1-1-1 always leads
to the two-loop neutrino mass diagram W-2-1-2-1-11-7 independent of the value
of $\alpha$. After applying the method to all the mass models up to two-loop
level, one can find that the model N1-1-1-1-1-1-1 can generate neutrino masses
when $\alpha=-7/6,-1/6,1/3,5/6,-\alpha_{2}$. In a similar way, we consider the
neutrino mass diagrams for each $0\nu\beta\beta$ decay model.

##  4  An example model of one-loop $0\nu\beta\beta$ decay

Figure 10:  An example of $0\nu\beta\beta$ decay model N1-1-6-1-4-1-1 (left panel) with the hypercharge parameter $\alpha=0$, which is a UV completion of the operator N1 based on the diagram N-1-6-1. The neutrino masses arise from the two-loop diagram W-2-5-4-1-3-10 (right panel) at leading order, and it is mediated by the messenger fields of the model N1-1-6-1-4-1-1.  | N1-1-6-1-4-1-1
---|---
Fields | $F_{1}$ | $F_{2}$ | $F_{3}$ | $S_{1}$ | $S_{2}$ | $S_{3}$
$C_{L}^{Y}$ | $1_{2}^{-1/6}$ | $3_{2}^{-1/2}$ | $3_{1}^{1/3}$ | $1_{1}^{-1/3}$ | $\bar{3}_{1}^{0}$ | $1_{2}^{-1/6}$
Table 4:  The mediators of the $0\nu\beta\beta$ decay model shown in the left
panel of figure 10 as well as their transformation properties under the SM
gauge group. The quantum numbers are given in the notation $C_{L}^{Y}$, where
$C$ refers to the $SU(3)_{C}$ transformation, $L$ refers to the $SU(2)_{L}$
transformation, and $Y$ stands for $U(1)_{Y}$ charge.

In this section, we present an example of model N1-1-6-1-4-1-1 which is a
decomposition of the $0\nu\beta\beta$ decay operator N1 based on the diagram
N-1-6-1 in figure 2. We add six new fields to the SM particle content in this
model, and their quantum number assignments are listed in table 4. The
interaction Lagrangian for these fields should respect the SM gauge symmetry,
and it is of the following form

$\displaystyle\mathcal{L}_{\text{int}}$ $\displaystyle=$ $\displaystyle
y_{I\alpha}\overline{Q_{L}}_{,\alpha}i\tau^{2}F_{1}^{c}S_{2}^{*}+y_{II\alpha}\overline{\ell}_{\alpha}F_{1}S_{1}+y_{III\alpha}\overline{u_{R}}_{,\alpha}F_{3}S_{1}^{*}+y_{IV\alpha}S_{3}^{\dagger}\overline{F_{3}}Q_{L,\alpha}$
(16)
$\displaystyle+y_{V\alpha}\overline{F_{2}}d_{R,\alpha}S_{3}+y_{VI\alpha}\overline{\ell}_{\alpha}F_{2}S_{2}+\sum_{i=1}^{3}m_{F_{i}}\overline{F_{i}}F_{i}+\mathrm{H.c.}\,,$

where $\tau^{2}$ is the Pauli matrix, and $i\tau^{2}=\epsilon$ is the
2-component antisymmetric tensor. This Lagrangian can generate the one-loop
diagram for $0\nu\beta\beta$ decay shown in figure 10. The gauge invariant
scalar potential terms involving new scalars are given by

$V\supset\left({\mu}S_{1}S_{3}^{T}i\tau^{2}H+\mathrm{H.c.}\right)+\sum_{i=1}^{3}{\lambda}_{S_{i}H}S_{i}^{\dagger}S_{i}H^{\dagger}H+\sum_{i,j=1}^{3}{\lambda}_{S_{i}S_{j}}S_{i}^{\dagger}S_{i}S_{j}^{\dagger}S_{j}+\sum_{i=1}^{3}m_{S_{i}}^{2}S_{i}^{\dagger}S_{i}\,.$
(17)

We denote $S_{3}=(S_{3,1},S_{3,2})^{T}$, the first term in Eq. (17) can induce
mixing between the two states $S_{1}$ and $S_{3,1}$ after electroweak symmetry
breaking. The bilinear terms of $S_{1}$ and $S_{3,1}$ read as

$\mu
v_{H}S_{1}S_{3,1}+\mu^{*}v_{H}S^{\dagger}_{1}S^{\dagger}_{3,1}+m_{S_{1}}^{2}S_{1}^{\dagger}S_{1}+m_{S_{3}}^{2}S_{3,1}^{\dagger}S_{3,1}=\left(S_{1}^{\dagger}\quad
S_{3,1}\right)\left(\begin{array}[]{cc}m_{S_{1}}^{2}{}&\mu^{*}\,v_{H}\\\
{\mu}\,v_{H}{}&m_{S_{3}}^{2}\end{array}\right)\left(\begin{array}[]{c}S_{1}\\\
S_{3,1}^{\dagger}\end{array}\right)\,,$ (18)

where $v_{H}$ is the vacuum expectation value of the Higgs field with $\langle
H^{0}\rangle=v_{H}$, we have redefined $m_{S_{1}}^{2}$ and $m_{S_{3}}^{2}$ to
absorb the contributions of ${\lambda}_{S_{1}H}v^{2}_{H}$ and
${\lambda}_{S_{3}H}v^{2}_{H}$ respectively. Assuming the mixing parameter
$\mu$ to be real, the above mass matrix can be diagonalized by a two-
dimensional rotation matrix

$\left(\begin{array}[]{cc}\cos\theta{}&\sin\theta\\\
-\sin\theta{}&\cos\theta\end{array}\right)\left(\begin{array}[]{cc}m_{S_{1}}^{2}{}&\mu\,v_{H}\\\
{\mu}\,v_{H}{}&m_{S_{3}}^{2}\end{array}\right)\left(\begin{array}[]{cc}\cos\theta{}&-\sin\theta\\\
\sin\theta{}&\cos\theta\end{array}\right)=\text{diag}(m^{2}_{S_{n1}},m^{2}_{S_{n2}})\,,$
(19)

where

$\displaystyle\tan 2\theta$ $\displaystyle=$
$\displaystyle\frac{2{\mu}\,v_{H}}{m_{S_{1}}^{2}-m_{S_{3}}^{2}}\,,$
$\displaystyle m^{2}_{S_{n1}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[m_{S_{1}}^{2}+m_{S_{3}}^{2}+\sqrt{(m_{S_{1}}^{2}-m_{S_{3}}^{2})^{2}+4\mu^{2}v^{2}_{H}}\right]\,,$
$\displaystyle m^{2}_{S_{n2}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[m_{S_{1}}^{2}+m_{S_{3}}^{2}-\sqrt{(m_{S_{1}}^{2}-m_{S_{3}}^{2})^{2}+4\mu^{2}v^{2}_{H}}\right]\,.$
(20)

Accordingly the mass eigenstates are related to the interaction basis by

$\left(\begin{array}[]{c}S_{n1}\\\
S_{n2}\end{array}\right)=\left(\begin{array}[]{cc}\cos\theta{}&\sin\theta\\\
-\sin\theta{}&\cos\theta\end{array}\right)\left(\begin{array}[]{c}S_{1}\\\
S_{3,1}^{*}\end{array}\right)\,.$ (21)

Rotating the interaction in Eq. (16) to the scalar mass eigenstate basis, the
terms relevant to $0\nu\beta\beta$ decay are given by

$\displaystyle\mathcal{L}_{\text{int}}=$ $\displaystyle
y_{I1}\overline{u_{L}}F_{1,2}^{c}S_{2}^{*}+y_{IIe}\overline{e_{L}}F_{1,2}(S_{n1}\cos\theta-
S_{n2}\sin\theta)+y_{III1}\overline{u_{R}}F_{3}(S_{n1}^{*}\cos\theta-
S_{n2}^{*}\sin\theta)$
$\displaystyle+y_{IV1}\overline{F_{3}}d_{L}S_{3,2}^{*}+y_{V1}\overline{F_{2,2}}d_{R}S_{3,2}+y_{VIe}\overline{e_{L}}F_{2,2}S_{2}+\mathrm{H.c.}+\ldots\,.$
(22)

The one-loop diagram in figure 10 can be straightforwardly calculated, and we
find the following effective operator is generated

$\frac{G_{F}^{2}\cos^{2}\theta_{C}}{2m_{P}}\epsilon_{1}^{LRL}J_{L}J_{R}j_{L}\,,$
(23)

where $J_{R,L}=\bar{u}(1\pm\gamma_{5})d$, $j_{L}=\bar{e}(1+\gamma_{5})e^{c}$
are defined in Eq. (38), and the coefficient $\epsilon_{1}^{LRL}$ is given by

$\displaystyle\epsilon_{1}^{LRL}=\frac{m_{F_{2}}m_{F_{3}}m_{P}}{128\pi^{2}m^{7}_{F_{1}}G_{F}^{2}\cos^{2}\theta_{C}}y_{I1}y_{IIe}y_{III1}y_{IV1}y_{V1}y_{VIe}\Bigg{[}\cos^{2}\theta\,\mathcal{F}\left(1,\frac{m_{F_{2}}^{2}}{m_{F_{1}}^{2}},\frac{m_{F_{3}}^{2}}{m_{F_{1}}^{2}},\frac{m_{S_{n1}}^{2}}{m_{F_{1}}^{2}},\frac{m_{S_{2}}^{2}}{m_{F_{1}}^{2}},\frac{m_{S_{3}}^{2}}{m_{F_{1}}^{2}}\right)$
$\displaystyle\hskip
32.52127pt+\sin^{2}\theta\,\mathcal{F}\left(1,\frac{m_{F_{2}}^{2}}{m_{F_{1}}^{2}},\frac{m_{F_{3}}^{2}}{m_{F_{1}}^{2}},\frac{m_{S_{n2}}^{2}}{m_{F_{1}}^{2}},\frac{m_{S_{2}}^{2}}{m_{F_{1}}^{2}},\frac{m_{S_{3}}^{2}}{m_{F_{1}}^{2}}\right)\Bigg{]}\,.$
(24)

The function $\mathcal{F}(x_{1},x_{2},x_{3},x_{4},x_{5},x_{6})$ is the loop
integral and it is given by

$\displaystyle\mathcal{F}(x_{1},x_{2},x_{3},x_{4},x_{5},x_{6})$
$\displaystyle=$ $\displaystyle\frac{1}{i\pi^{2}}\int
d^{4}p\frac{1}{\left(p^{2}-x_{1}\right)\left(p^{2}-x_{2}\right)\left(p^{2}-x_{3}\right)\left(p^{2}-x_{4}\right)\left(p^{2}-x_{5}\right)\left(p^{2}-x_{6}\right)}$
(25) $\displaystyle=$ $\displaystyle\frac{x_{1}\left(1-\ln
x_{1}\right)}{\left(x_{1}-x_{2}\right)\left(x_{1}-x_{3}\right)\left(x_{1}-x_{4}\right)\left(x_{1}-x_{5}\right)\left(x_{1}-x_{6}\right)}$
$\displaystyle+\frac{x_{2}\left(1-\ln
x_{2}\right)}{\left(x_{2}-x_{1}\right)\left(x_{2}-x_{3}\right)\left(x_{2}-x_{4}\right)\left(x_{2}-x_{5}\right)\left(x_{2}-x_{6}\right)}$
$\displaystyle+\frac{x_{3}\left(1-\ln
x_{3}\right)}{\left(x_{3}-x_{1}\right)\left(x_{3}-x_{2}\right)\left(x_{3}-x_{4}\right)\left(x_{3}-x_{5}\right)\left(x_{3}-x_{6}\right)}$
$\displaystyle+\frac{x_{4}\left(1-\ln
x_{4}\right)}{\left(x_{4}-x_{1}\right)\left(x_{4}-x_{2}\right)\left(x_{4}-x_{3}\right)\left(x_{4}-x_{5}\right)\left(x_{4}-x_{6}\right)}$
$\displaystyle+\frac{x_{5}\left(1-\ln
x_{5}\right)}{\left(x_{5}-x_{1}\right)\left(x_{5}-x_{2}\right)\left(x_{5}-x_{3}\right)\left(x_{5}-x_{4}\right)\left(x_{5}-x_{6}\right)}$
$\displaystyle+\frac{x_{6}\left(1-\ln
x_{6}\right)}{\left(x_{6}-x_{1}\right)\left(x_{6}-x_{2}\right)\left(x_{6}-x_{3}\right)\left(x_{6}-x_{4}\right)\left(x_{6}-x_{5}\right)}\,.$

In the limit
$m_{F_{1}}=m_{F_{2}}=m_{F_{3}}=m_{S_{1}}=m_{S_{2}}=m_{S_{3}}\equiv\bar{m}$, we
have $x_{1}=x_{2}=x_{3}=x_{5}=x_{6}=1$,
$x_{4}=m_{S_{n1}}^{2}/m_{F_{1}}^{2}=1+\mu v_{H}/\bar{m}^{2}$ or
$x_{4}=m_{S_{n2}}^{2}/m_{F_{1}}^{2}=1-\mu v_{H}/\bar{m}^{2}$, Eq. (25) is not
valid anymore. In this case, we can use the following formula instead,

$\mathcal{F}(1,1,1,x_{4},1,1)=\frac{3+10x_{4}-18x_{4}^{2}+6x_{4}^{3}-x_{4}^{4}+12x_{4}\ln
x_{4}}{12\left(1-x_{4}\right)^{5}}\,.$ (26)

At present, the most stringent experimental limits on the $0\nu\beta\beta$
decay half-life of Germanium and Xenon are $T_{1/2}(^{76}\text{Ge})>1.8\times
10^{26}$ yr [35] and $T_{1/2}(^{136}\text{Xe})>1.07\times 10^{26}$ yr [3]
respectively at $90\%$ confidence level, which lead to the upper limits on the
short-range coupling $|\epsilon_{1}^{LRL}|<2.90\times 10^{-10}$ and
$|\epsilon_{1}^{LRL}|<2.50\times 10^{-10}$ respectively [36] if only the
short-range mechanism is considered.

###  4.1 Prediction for neutrino masses

We have checked that for this model the leading order neutrino masses arise
from a two-loop model W-2-5-4-1-3-10, as shown in figure 10. The relevant
interaction Lagrangian is given by

$\displaystyle\mathcal{L}^{\nu}_{\text{int}}$ $\displaystyle=$ $\displaystyle
y_{I\alpha}\overline{Q_{L}}_{,\alpha}i\tau^{2}F_{1}^{c}S_{2}^{*}+y_{II\alpha}\overline{\ell}_{\alpha}F_{1}S_{1}+y_{d,\alpha\beta}\overline{Q_{L}}_{,\alpha}d_{R,\beta}H$
(27)
$\displaystyle+y_{V\alpha}\overline{F_{2}}d_{R,\alpha}S_{3}+y_{VI\alpha}\overline{\ell}_{\alpha}F_{2}S_{2}+{\mu}S_{1}S_{3}^{T}i\tau^{2}H+\mathrm{H.c.}\,.$

The lepton number is violated by the simultaneous presence of the above terms.
After the electroweak symmetry breaking, the above interaction can be written
in the mass eigenstate basis as follows,

$\displaystyle\mathcal{L}^{\nu}_{\text{int}}$ $\displaystyle=$ $\displaystyle
y_{I\alpha}\overline{d_{L}}_{,\alpha}F^{c}_{1,1}S_{2}^{*}+y_{II\alpha}\overline{\nu_{L}}_{,\alpha}F_{1,1}(S_{n1}\cos\theta-
S_{n2}\sin\theta)$ (28)
$\displaystyle+y_{V\alpha}\overline{F_{2}}_{,1}d_{R,\alpha}(S_{n1}^{*}\sin\theta+S_{n2}^{*}\cos\theta)+y_{VI\alpha}\overline{\nu_{L}}_{,\alpha}F_{2,1}S_{2}+\mathrm{H.c.}+\ldots\,,$

where the unitary rotations of the down type quarks to mass eigenstates have
been absorbed into the couplings $y_{I\alpha}$ and $y_{V\alpha}$. The neutrino
mass matrix arises from the two-loop diagrams of figure 11 in the mass basis
and it can be straightforwardly calculated as

$\displaystyle(m_{\nu})_{\alpha\beta}$ $\displaystyle=$
$\displaystyle\frac{3}{512\pi^{4}}\sin
2\theta\,(y^{*}_{II\alpha}y^{*}_{VI\beta}+y^{*}_{II\beta}y^{*}_{VI\alpha})y^{*}_{Ik}y^{*}_{Vk}\frac{m_{F_{1}}m_{F_{2}}m_{d_{k}}}{m^{2}_{S_{2}}}\bigg{[}\mathcal{I}\left(\frac{m_{F_{1}}^{2}}{m^{2}_{S_{2}}},\frac{m_{S_{n2}}^{2}}{m^{2}_{S_{2}}},\frac{m_{F_{2}}^{2}}{m^{2}_{S_{2}}},1,\frac{m_{d_{k}}^{2}}{m^{2}_{S_{2}}}\right)$
(29)
$\displaystyle-\mathcal{I}\left(\frac{m_{F_{1}}^{2}}{m^{2}_{S_{2}}},\frac{m_{S_{n1}}^{2}}{m^{2}_{S_{2}}},\frac{m_{F_{2}}^{2}}{m^{2}_{S_{2}}},1,\frac{m_{d_{k}}^{2}}{m^{2}_{S_{2}}}\right)\bigg{]}\,.$

Here the loop integral
$\mathcal{I}\left(r_{a},t_{\alpha},r_{b},t_{\beta},r_{c}\right)$ is

$\displaystyle\mathcal{I}\left(r_{a},t_{\alpha},r_{b},t_{\beta},r_{c}\right)=\left(\frac{1}{i\pi^{2}}\right)^{2}\int
d^{4}k\int
d^{4}q\frac{1}{\left(k^{2}-r_{a}\right)\left(k^{2}-t_{\alpha}\right)\left(q^{2}-r_{b}\right)\left(q^{2}-t_{\beta}\right)\left[\left(q+k\right)^{2}-r_{c}\right]}$
$\displaystyle\qquad\qquad\qquad=\frac{1}{r_{a}-t_{\alpha}}\frac{1}{r_{b}-t_{\beta}}\left[T_{3}(r_{a},r_{b},r_{c})-T_{3}(r_{a},t_{\beta},r_{c})-T_{3}(t_{\alpha},r_{b},r_{c})+T_{3}(t_{\alpha},t_{\beta},r_{c})\right]\,,$
(30)

with

$\displaystyle
T_{3}(x_{1},x_{2},x_{3})=\left(\frac{1}{i\pi^{2}}\right)^{2}\int d^{4}k\int
d^{4}q\frac{1}{\left(k^{2}-x_{1}\right)\left(q^{2}-x_{2}\right)\left[\left(q+k\right)^{2}-x_{3}\right]}\,.$
(31)

The momentum space integration $T_{3}(x_{1},x_{2},x_{3})$ can be performed
through the introduction of Feynman parameters, and it has been calculated in
Refs. [37, 38, 39, 40]. The analytical expression of
$\mathcal{I}\left(r_{a},t_{\alpha},r_{b},t_{\beta},r_{c}\right)$ reads as

$\displaystyle\mathcal{I}\left(r_{a},t_{\alpha},r_{b},t_{\beta},r_{c}\right)$
$\displaystyle=$
$\displaystyle\frac{r_{c}}{(r_{a}-t_{\alpha})(r_{b}-t_{\beta})}\big{[}g(r_{a}/r_{c},r_{b}/r_{c})-g(r_{a}/r_{c},t_{\beta}/r_{c})$
(32)
$\displaystyle-g(t_{\alpha}/r_{c},r_{b}/r_{c})+g(t_{\alpha}/r_{c},t_{\beta}/r_{c})\Big{]}\,,$

where

$\displaystyle g(s,t)$ $\displaystyle=$ $\displaystyle\frac{s}{2}\ln s\ln
t+\sum_{\pm}\pm\frac{s(1-s)+3st+2(1-t)x_{\pm}}{2\omega}$ (33)
$\displaystyle\times\left[\mathrm{Li}_{2}\left(\frac{x_{\pm}}{x_{\pm}-s}\right)-\mathrm{Li}_{2}\left(\frac{x_{\pm}-s}{x_{\pm}}\right)+\mathrm{Li}_{2}\left(\frac{t-1}{x_{\pm}}\right)-\mathrm{Li}_{2}\left(\frac{t-1}{x_{\pm}-s}\right)\right]$

with $\mathrm{Li}_{2}(x)$ being the dilogarithm function

$\mathrm{Li}_{2}(x)=-\int_{0}^{x}\frac{\ln(1-y)}{y}dy$ (34)

and

$x_{\pm}=\frac{1}{2}(-1+s+t\pm\omega)\,,\quad\omega=\sqrt{1+s^{2}+t^{2}-2(s+t+st)}~{}\,.$
(35)

The neutrino mass in Eq. (29) is dominated by the contribution of the bottom
quark because the hierarchical mass pattern $m_{b}\gg m_{s}\gg m_{d}$. Only
one generation for every new fermions and scalars are introduced in our model,
therefore the neutrino mass matrix of Eq. (29) is of rank two so that the
lightest neutrino is massless. If one allows for more than one generation of
the new fermions and scalars, the corresponding neutrino mass matrix can be
read out analogously and the mass of the lightest neutrino could be non-zero.

Figure 11:  The diagrams for the neutrino mass generation in the model
W-2-5-4-1-3-10, which is mediated by the fields of the $0\nu\beta\beta$ decay
model N1-1-6-1-4-1-1, where $d_{k}$ with $k=1,2,3$ stand for the down quark,
strange quark and bottom quark respectively.

###  4.2 Half-life time of $0\nu\beta\beta$ decay

As shown in previous subsection, the mediators of the $0\nu\beta\beta$ decay
model can generate non-vanishing light neutrino masses at two-loop level.
Therefore both the short-range operator in Eq. (23) and the mass mechanism of
the light neutrino exchange contribute to the $0\nu\beta\beta$ decay, and the
two contributions should be added coherently. Using the general formula of Eq.
(40) for the inverse half-life, we have

$\displaystyle T_{1/2}^{-1}$ $\displaystyle=$ $\displaystyle
G_{11+}^{(0)}\absolutevalue{\frac{m_{\beta\beta}}{m_{e}}\mathcal{M}_{\nu}+\epsilon_{1}^{LRL}\mathcal{M}^{LR}_{1}}^{2}$
(36) $\displaystyle=$ $\displaystyle
G_{11+}^{(0)}\left[\absolutevalue{\mathcal{M}_{\nu}}^{2}\frac{\absolutevalue{m_{\beta\beta}}^{2}}{m^{2}_{e}}+\absolutevalue{\mathcal{M}^{LR}_{1}}^{2}\absolutevalue{\epsilon_{1}^{LRL}}^{2}+2\absolutevalue{\frac{m_{\beta\beta}}{m_{e}}\mathcal{M}_{\nu}}\absolutevalue{\epsilon_{1}^{LRL}\mathcal{M}^{LR}_{1}}\cos(\alpha-\beta)\right]\,,$

where $\alpha$ and $\beta$ denote the complex phases of $m_{\beta\beta}$ and
$\epsilon_{1}^{LRL}$ respectively, note that the nuclear matrix elements are
conventionally defined to be real. We see that the interference between the
short-range contribution and mass mechanism is constructive for
$\alpha-\beta=0$ and destructive for $\alpha-\beta=\pi$. It is remarkable that
the two contributions would exactly cancel each other if
$\epsilon_{1}^{LRL}=-m_{\beta\beta}\mathcal{M}_{\nu}/(m_{e}\mathcal{M}^{LR}_{1})$.

We show the constraints on the effective Majorana neutrino mass
$\absolutevalue{m_{\beta\beta}}$ and the short-range coupling
$|\epsilon_{1}^{LRL}|$ in figure 12, where the values of phase space factor
and nuclear matrix elements are adopted from Ref. [36]. The highlighted areas
denote the allowed regions by the current limits and future sensitivities of
the $0\nu\beta\beta$ decay half-life of the isotopes ${}^{76}\text{Ge}$ and
${}^{136}\text{Xe}$, where the phase $\alpha-\beta$ freely varies between $0$
to $\pi$. The half-life sensitivity would be increased by a factor of
approximately one hundred in upcoming generation of tonne-scale experiments,
therefore the constraint on parameter space would be improved considerably, as
can be seen from figure 12. Moreover, we notice that the constraint imposed by
${}^{136}\text{Xe}$ is a bit more stringent than that of ${}^{76}\text{Ge}$.
Cosmology has played an important role in constraining neutrino properties,
and the analysis of cosmological observables can provide an upper bound on a
total neutrino mass. Within the standard minimal $\Lambda$CDM model, the
present cosmological bound on neutrino mass is $\sum_{i}m_{\nu_{i}}<0.12$ eV
for Planck TT, TE, EE+lowE+lensing+BAO [41]. Taking into account the values of
the neutrino mass squared differences and mixing angles measured by the
neutrino oscillation experiments [42], one can obtain that the effective
Majorana neutrino mass is in the region: $|m_{\beta\beta}|\leq 31.12$ meV for
normal ordering (NO) and $18.62\,\text{meV}\leq|m_{\beta\beta}|\leq 51.14$ meV
for inverted ordering (IO). With this constraint the current bounds
$T_{1/2}(^{76}\text{Ge})>1.8\times 10^{26}$ yr [35] and
$T_{1/2}(^{136}\text{Xe})>1.07\times 10^{26}$ yr [3] require the coupling
$|\epsilon_{1}^{LRL}|$ smaller than $4\times 10^{-10}$. It is remarkable that
the forthcoming tonne-scale searches for $0\nu\beta\beta$ decays would
constrain the parameter space into a narrow band where the short-range
mechanism is comparable with the mass mechanism.

We proceed to discuss the relative size of the short-range mechanism and mass
mechanism to the decay rate. We take
$m_{S_{1}}=m_{S_{2}}=m_{S_{3}}=m_{F_{1}}=m_{F_{2}}=m_{F_{3}}=1$ TeV and the
coupling constants $y_{Ik}=y_{III1}=y_{IIe}=y_{IV1}=y_{Vk}=y_{VIe}\equiv
y_{eff}$ to be real. We display the contour region of
$\log_{10}\big{(}\big{|}\frac{\epsilon^{\text{LRL}}_{1}\mathcal{M}^{\text{LR}}_{1}}{\epsilon_{\nu}\mathcal{M}_{\nu}}\big{|}\big{)}$
in the plane $\mu$ versus $y_{eff}$ in figure 13, where $\epsilon_{\nu}\equiv
m_{\beta\beta}/m_{e}$. The parameter $\mu$ dictates the mass splitting of
$m^{2}_{S_{n1}}$ and $m^{2}_{S_{n2}}$, as shown in Eq. (20). If $\mu$ is quite
small with $\mu/\text{TeV}\ll 1$, $m^{2}_{S_{n1}}$ and $m^{2}_{S_{n2}}$ tend
to be degenerate so that cancellation between two terms of the light neutrino
mass matrix in Eq. (29) occurs. As a result, the mass mechanism is generally
sub-dominant to the short-range contribution for $\mu/\text{TeV}\ll 1$.
Furthermore, from Eq. (24) and Eq. (29) we see $\epsilon_{1}^{LRL}\propto
Y^{6}_{eff}$ and $(m_{\nu})_{\alpha\beta}\propto Y^{4}_{eff}$. Therefore the
mass mechanism is dominant in the region $y_{eff}\ll 1$ while the short-range
mechanism is dominant in the region $y_{eff}>1$.

Figure 12:  Constraints on the effective neutrino mass $|m_{\beta\beta}|$ and
the short-range coupling $\epsilon_{1}^{LRL}$ for normal ordering and inverted
ordering neutrino masses. The black dashed lines from top to bottom correspond
to
$\absolutevalue{\epsilon_{1}^{LRL}\mathcal{M}^{LR}_{1}}/\absolutevalue{\epsilon_{\nu}\mathcal{M}_{\nu}}=10,\;1,\;0.1$
respectively. The highlighted regions represent the parameter space allowed by
the current bounds $T_{1/2}(^{76}\text{Ge})>1.8\times 10^{26}$ yr [35],
$T_{1/2}(^{136}\text{Xe})>1.07\times 10^{26}$ yr [3] and the future
sensitivities $T_{1/2}(^{76}\text{Ge})>1.3\times 10^{28}$ yr [4],
$T_{1/2}(^{136}\text{Xe})>1.35\times 10^{28}$ yr [5]. The gray bands denote
the regions excluded by the cosmological data $\sum_{i=1}^{3}m_{\nu_{i}}<0.12$
eV at $95\%$ confidence level from the Planck Collaboration [41], and the
oscillation parameters are fixed to the current best fit values [42].  Figure
13:  The relative size of the short-range mechanism and mass mechanism to the
decay rate in the plane $\mu$ versus the couplings $y_{eff}$. The mass
mechanism is dominant in the blue region while the short-range mechanism is
dominant in the dark red region. The shaded regions denote the contours of the
current limit and future sensitivities in ${}^{136}\text{Xe}$ (left) and
${}^{76}\text{Ge}$ (right) $0\nu\beta\beta$ decay search, where the phase
$\alpha-\beta$ freely varies between $0$ and $\pi$.

##  5 Conclusion

The $0\nu\beta\beta$ decay is the most sensitive experiments to establish that
the massive neutrinos are Majorana or Dirac particles. If the signal of
$0\nu\beta\beta$ decay is observed in the future, lepton number would be
violated by two units and neutrinos must be Majorana particles. The rare
$0\nu\beta\beta$ decay process might also provide important information on the
neutrino mass spectrum, the Majorana phases and the absolute scale of neutrino
masses. On the other hand, if neutrinos are Majorana particles,
$0\nu\beta\beta$ decay can be mediated by exchange of a light neutrino between
two charged current vertices, and it is the so-called mass mechanism.
Additional contributions from short-range mechanism is also possible, where
the exchanged particles are much heavier than the $0\nu\beta\beta$ decay
energy scale typically 100 MeV, and they can be represented by non-
renormalizable contact interactions with six external fermions.

In the present work, we have systematically studied the possible one-loop UV
completions of dimension-9 short-range $0\nu\beta\beta$ decay operators in Eq.
(3). The general procedure to be followed involve three steps: firstly
generate one-loop topologies with six legs, then promote topology to diagram
by specifying the Lorentz nature of fields, finally attach fields to external
lines and determine the SM quantum numbers of internal particles by gauge
invariance. Notice that there are infinite possibilities for the quantum
numbers of the internal fermions and scalars, thus a given diagram can lead to
an infinite number of UV completions. To be more specific, excluding the
topologies with tadpole and self-energy, there are only 48 one-loop connected
topologies with three-point and four-point vertices and six external lines, as
shown in figure 1, figure 14 and figure 15. Only 11 of them in figure 1 and
figure 15 can accommodate 6 external fermion lines using renormalizable
interactions. The 5 topologies in figure 15 contain a one-loop sub-diagram
which can be compressed into a three-point vertex so that the more important
tree-level contributions can not be forbidden. As a result, it is sufficient
to only consider the 6 one-loop topologies in figure 1 which lead to 8 one-
loop genuine diagrams shown in figure 2. One can determine the transformations
of the internal fields under SM gauge group by imposing gauge invariance at
each vertex, the $SU(3)_{C}$, $SU(2)_{L}$ representations and hypercharge
$U(1)_{Y}$ assignments can be unambiguously fixed in tree-level realizations
while there are an infinite number of UV completions at one-loop level. The
complete list of the models together with the internal fields and their
quantum numbers are given in the supplementary file [32], where the
$SU(2)_{L}$ and $SU(3)_{C}$ representations with dimensions less than 4 and 10
respectively are considered.

In order to guarantee the dominance of the one-loop contribution to
$0\nu\beta\beta$ decay, the mediators of the tree diagrams should be absent in
the one-loop models. As a consequence, certain values of the hypercharge
parameter $\alpha$ are excluded. The $0\nu\beta\beta$ decay models violate
lepton number, thus they can also generate Majorana neutrino masses, as is
required by the black box theorem [2]. Therefore both short-range contribution
and mass mechanism are present in these models. We have discussed the Majorana
neutrino masses up to two-loop level which can be generated from the mediators
of the $0\nu\beta\beta$ decay models. Furthermore, we give a concrete one-loop
$0\nu\beta\beta$ decay model in which the neutrino masses are generated at
2-loop level. The predictions for the neutrino mass matrix and the half-life
of $0\nu\beta\beta$ decay are presented. The short-range contribution can be
as important as the mass mechanism in some region of parameter space.

There are many experiments that are searching for $0\nu\beta\beta$ decay or
are in various stages of planning and construction over the world. We expect
it would become an important probe for physics beyond the SM. If the
observation of $0\nu\beta\beta$ decay experiments is in conflict with neutrino
mass ordering measured by neutrino oscillation experiments and the neutrino
mass bounds from $\beta$ decay and cosmology in the future, the new
contributions beyond the mass mechanism would be necessary. Then our results
can be used to construct short-range $0\nu\beta\beta$ decay models.
Furthermore, generally a few messenger fields are required in the one-loop
$0\nu\beta\beta$ decay models, it is interesting to study the phenomenological
implications of these new fields in collider and cosmology, and this is left
for future work.

## Acknowledgements

PTC and GJD are supported by the National Natural Science Foundation of China
under Grant Nos 11975224, 11835013 and 11947301 and the Key Research Program
of the Chinese Academy of Sciences under Grant NO. XDPB15. CYY is supported in
part by the Grants No. NSFC-11975130, No. NSFC-12035008, No. NSFC-12047533, by
the National Key Research and Development Program of China under Grant No.
2017YFA0402200 and the China Post-doctoral Science Foundation under Grant No.
2018M641621.

## Appendix

##  Appendix A  $0\nu\beta\beta$ decay operators below the electroweak scale
and half-life

The SM gauge symmetry $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}$ is
spontaneous broken down to $SU(3)_{C}\times U(1)_{EM}$. The leading operators
which contribute to $0\nu\beta\beta$ decay at short distance involve four
quark fields and two charged leptons, they respect $SU(3)_{C}\times U(1)_{EM}$
symmetry and can be written as the product of three fermion currents [6, 36]

$\displaystyle(\mathcal{O}_{1})_{\\{LL\\}L}=J_{L}J_{L}j_{L},$
$\displaystyle(\mathcal{O}_{1})_{\\{LL\\}R}=J_{L}J_{L}j_{R}\,,$
$\displaystyle(\mathcal{O}_{1})_{\\{RR\\}L}=J_{R}J_{R}j_{L},$
$\displaystyle(\mathcal{O}_{1})_{\\{RR\\}R}=J_{R}J_{R}j_{R}\,,$
$\displaystyle(\mathcal{O}_{1})_{\\{LR\\}L}=J_{L}J_{R}j_{L},$
$\displaystyle(\mathcal{O}_{1})_{\\{LR\\}R}=J_{L}J_{R}j_{R}\,,$
$\displaystyle(\mathcal{O}_{2})_{\\{LL\\}L}=J^{\mu\nu}_{L}J_{L\mu\nu}j_{L},$
$\displaystyle(\mathcal{O}_{2})_{\\{LL\\}R}=J^{\mu\nu}_{L}J_{L\mu\nu}j_{R}\,,$
$\displaystyle(\mathcal{O}_{2})_{\\{RR\\}L}=J^{\mu\nu}_{R}J_{R\mu\nu}j_{L},$
$\displaystyle(\mathcal{O}_{2})_{\\{RR\\}R}=J^{\mu\nu}_{R}J_{R\mu\nu}j_{R}\,,$
$\displaystyle(\mathcal{O}_{3})_{\\{LL\\}L}=J^{\mu}_{L}J_{L\mu}j_{L},$
$\displaystyle(\mathcal{O}_{3})_{\\{LL\\}R}=J^{\mu}_{L}J_{L\mu}j_{R}\,,$
$\displaystyle(\mathcal{O}_{3})_{\\{RR\\}L}=J^{\mu}_{R}J_{R\mu}j_{L},$
$\displaystyle(\mathcal{O}_{3})_{\\{RR\\}R}=J^{\mu}_{R}J_{R\mu}j_{R}\,,$
$\displaystyle(\mathcal{O}_{3})_{\\{LR\\}L}=J^{\mu}_{L}J_{R\mu}j_{L},$
$\displaystyle(\mathcal{O}_{3})_{\\{LR\\}R}=J^{\mu}_{L}J_{R\mu}j_{R}\,,$
$\displaystyle(\mathcal{O}_{4})_{LL}=J_{L}^{\mu}J_{L\mu\nu}j^{\nu},$
$\displaystyle(\mathcal{O}_{4})_{RR}=J_{R}^{\mu}J_{R\mu\nu}j^{\nu}\,,$
$\displaystyle(\mathcal{O}_{4})_{LR}=J_{L}^{\mu}J_{R\mu\nu}j^{\nu},$
$\displaystyle(\mathcal{O}_{4})_{RL}=J_{R}^{\mu}J_{L\mu\nu}j^{\nu}\,,$
$\displaystyle(\mathcal{O}_{5})_{LL}=J_{L}J_{L}^{\mu}j_{\mu},$
$\displaystyle(\mathcal{O}_{5})_{RR}=J_{R}J_{R}^{\mu}j_{\mu}\,,$
$\displaystyle(\mathcal{O}_{5})_{LR}=J_{L}J_{R}^{\mu}j_{\mu},$
$\displaystyle(\mathcal{O}_{5})_{RL}=J_{R}J_{L}^{\mu}j_{\mu}\,.$ (37)

Here the normalization factor $\frac{G^{2}_{F}\cos^{2}\theta_{C}}{2m_{P}}$ has
been omitted, where $m_{P}$ is the proton mass and $\theta_{C}$ is the Cabibbo
angle. The factor $\cos^{2}\theta_{C}$ is included so that the phase space
factor of short-range contribution can be defined in the same way as that for
the standard light neutrino exchange [36]. The hadronic current $J$ and
leptonic current $j$ are bilinears of quarks and electrons,

$\displaystyle J_{R,L}$ $\displaystyle=$
$\displaystyle\bar{u}_{a}(1\pm\gamma_{5})d_{a},\quad
J_{R,L}^{\mu}=\bar{u}_{a}\gamma^{\mu}(1\pm\gamma_{5})d_{a},\quad
J_{R,L}^{\mu\nu}=\bar{u}_{a}\sigma^{\mu\nu}(1\pm\gamma_{5})d_{a},$ (38)
$\displaystyle\qquad\qquad j_{R,L}=\bar{e}(1\mp\gamma_{5})e^{c},\quad
j^{\mu}=\bar{e}\gamma^{\mu}\gamma_{5}e^{c}\,,$

where $a=1,2,3$ is $SU(3)_{C}$ color indices and
$\sigma^{\mu\nu}=\frac{i}{2}[\gamma^{\mu},\gamma^{\nu}]$. The chirality
assignment in the charged lepton current $j_{R,L}$ is flipped, and the
subscripts $L$ and $R$ are associated with $1+\gamma_{5}$ and $1-\gamma_{5}$
respectively. Thus the operators $j_{R}$ and $j_{L}$ describes the creation of
two right-handed and left-handed electrons respectively. Notice that using the
Fierz identity one can show that both operators
$(\mathcal{O}_{2})_{\\{LR\\}L}=J^{\mu\nu}_{L}J_{R\mu\nu}j_{L}$ and
$(\mathcal{O}_{2})_{\\{LR\\}R}=J^{\mu\nu}_{L}J_{R\mu\nu}j_{R}$ are vanishing.
Hence the most general short-range interaction Lagrangian responsible for
$0\nu\beta\beta$ decay can be written as

$\mathcal{L}_{\text{short}}=\frac{G^{2}_{F}\cos^{2}\theta_{C}}{2m_{P}}\sum_{\chi_{1},\chi_{2},\chi}\left[\sum^{3}_{i=1}\epsilon_{i}^{\chi_{1}\chi_{2}\chi}(\mathcal{O}_{i})_{\\{\chi_{1}\chi_{2}\\}\chi}+\sum^{5}_{i=4}\epsilon_{i}^{\chi_{1}\chi_{2}}(\mathcal{O}_{i})_{\chi_{1}\chi_{2}}\right]+\mathrm{H.c.}\,,$
(39)

where the indices $\chi_{1},\chi_{2},\chi$ are $L$ or $R$ and they denote the
chiralities of the quark and electron currents in the order of appearance in
Eq. (37). Under the presence of both standard mass mechanism and the short-
range mechanism, the inverse $0\nu\beta\beta$ decay half-life time can be
expressed as [36]

$\displaystyle T_{1/2}^{-1}$ $\displaystyle=$ $\displaystyle
G_{11+}^{(0)}\absolutevalue{\sum_{I=1}^{3}\epsilon_{I}^{L}\mathcal{M}_{I}+\epsilon_{\nu}\mathcal{M}_{\nu}}^{2}+G_{11+}^{(0)}\absolutevalue{\sum_{I=1}^{3}\epsilon_{I}^{R}\mathcal{M}_{I}}^{2}+G_{66}^{(0)}\absolutevalue{\sum_{I=4}^{5}\epsilon_{I}\mathcal{M}_{I}}^{2}$
(40) $\displaystyle+G_{11-}^{(0)}\times
2\real\left[\left(\sum_{I=1}^{3}\epsilon_{I}^{L}\mathcal{M}_{I}+\epsilon_{\nu}\mathcal{M}_{\nu}\right)\left(\sum_{I=1}^{3}\epsilon_{I}^{R}\mathcal{M}_{I}\right)^{*}\right]$
$\displaystyle+G_{16}^{(0)}\times
2\real\left[\left(\sum_{I=1}^{3}\epsilon_{I}^{L}\mathcal{M}_{I}-\sum_{I=1}^{3}\epsilon_{I}^{R}\mathcal{M}_{I}+\epsilon_{\nu}\mathcal{M}_{\nu}\right)\left(\sum_{I=4}^{5}\epsilon_{I}\mathcal{M}_{I}\right)^{*}\right]\,,$

where $\epsilon_{\nu}\equiv m_{\beta\beta}/m_{e}$, the summation is over the
different types of short-range contributions including different chiralities
with $I=(i,\chi_{1}\chi_{2})$, $i=1,\ldots,5$ and
$\chi_{1},\chi_{2}\in\left\\{L,R\right\\}$. The notations $\epsilon_{I}^{L}$
and $\epsilon_{I}^{R}$ stand for the coefficients
$\epsilon^{\chi_{1}\chi_{2}L}_{i}$ and $\epsilon^{\chi_{1}\chi_{2}R}_{i}$
respectively when $i=1,2,3$. Moreover, $G_{11+}^{(0)}$, $G_{66}^{(0)}$,
$G_{11-}^{(0)}$, $G_{16}^{(0)}$ are the phase space factors,
$\mathcal{M}_{\nu}$ and $\mathcal{M}_{I}$ are the nuclear matrix elements for
the mass mechanism and the short-range operators respectively, their numerical
values within the interacting boson model are listed in Ref. [36].

As shown in Eq. (3), there are eleven $0\nu\beta\beta$ decay operators
invariant under the SM gauge symmetry $SU(3)_{C}\times SU(2)_{L}\times
U(1)_{Y}$ and twenty-four $SU(3)_{C}\times U(1)_{EM}$ invariant operators at
dimension $d=9$. The matching between the $0\nu\beta\beta$ decay operators at
the electroweak scale is listed in table 5. Multiplying the dimension-9
operators in Eq. (3) by $H^{\dagger}H$ and $(H^{\dagger}H)^{2}$, one can
obtain the remaining thirteen low-energy dimension-9 operators after
electroweak symmetry breaking. Following the general procedures listed in
section 2, the possible UV completions of the above effective operators in Eq.
(37) can be straightforwardly found out.

Above the electroweak scale | Below the electroweak scale
---|---
$\mathcal{O}^{0\nu}_{1}$ | $\frac{1}{8}(\mathcal{O}_{3})_{\\{LR\\}L}$
$\mathcal{O}^{0\nu}_{2}$ | $-\frac{1}{4}(\mathcal{O}_{1})_{\\{LR\\}L}-\frac{1}{24}(\mathcal{O}_{3})_{\\{LR\\}L}$
$\mathcal{O}^{0\nu}_{3}$ | $\frac{1}{8}(\mathcal{O}_{1})_{\\{LL\\}L}$
$\mathcal{O}^{0\nu}_{4}$ | $-\frac{5}{24}(\mathcal{O}_{1})_{\\{LL\\}L}-\frac{1}{32}(\mathcal{O}_{2})_{\\{LL\\}L}$
$\mathcal{O}^{0\nu}_{5}$ | $\frac{1}{8}(\mathcal{O}_{1})_{\\{RR\\}L}$
$\mathcal{O}^{0\nu}_{6}$ | $-\frac{5}{24}(\mathcal{O}_{1})_{\\{RR\\}L}-\frac{1}{32}(\mathcal{O}_{2})_{\\{RR\\}L}$
$\mathcal{O}^{0\nu}_{7}$ | $\frac{1}{8}(\mathcal{O}_{3})_{\\{RR\\}R}$
$\mathcal{O}^{0\nu}_{8}$ | $\frac{1}{8}(\mathcal{O}_{5})_{RR}$
$\mathcal{O}^{0\nu}_{9}$ | $\frac{i}{8}(\mathcal{O}_{4})_{RR}-\frac{5}{24}(\mathcal{O}_{5})_{RR}$
$\mathcal{O}^{0\nu}_{10}$ | $\frac{1}{8}(\mathcal{O}_{5})_{RL}$
$\mathcal{O}^{0\nu}_{11}$ | $-\frac{i}{8}(\mathcal{O}_{4})_{RL}-\frac{5}{24}(\mathcal{O}_{5})_{RL}$
Table 5:  Matching between the $0\nu\beta\beta$ decay operators above the
electroweak scale and those below the electroweak scale.

##  Appendix B Non-renormalizable topologies and finite non-genuine diagrams

For completeness, we list the topologies which always lead to the non-
renormalizable diagrams in figure 14 as well as the non-genuine topologies in
figure 15. The non-genuine model arising from the topologies in figure 15 can
be regarded as extensions of the tree-level $0\nu\beta\beta$ decay models, and
one of the vertices is generated at one-loop level, as illustrated in figure
16 and figure 17.

Figure 14: The non-renormalizable tree-level and one-loop topologies. Figure
15: The one-loop topologies leading to non-genuine diagrams. Figure 16:  Non-
genuine one-loop diagrams which correspond to the one-loop generation of one
of the couplings in the tree-level diagram N-0-1-1(see figure 2). The diagrams
in the box are divergent.  Figure 17:  Non-genuine one-loop diagrams which
correspond to the one-loop generation of one of the couplings in the tree-
level diagram N-0-2-1 (see figure 2). The diagrams in the box are divergent.

##  Appendix C Possible UV completions for the $0\nu\beta\beta$ decay operator
N4 based on the diagram N-1-1-1

In this appendix, we present the possible one-loop models which generate the
$0\nu\beta\beta$ decay operator N4, and the models are based on the diagram
N-1-1-1. The SM quantum numbers of the messenger fields are listed in tables
6-10.

Model Name | $U(1)_{Y}$ | $SU(2)_{L}$ | $SU(3)_{C}$
---|---|---|---
$F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$ | $F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$ | $F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$
N4-1-1-1-1-1-1 | $-\frac{1}{3}+\alpha$ | $\alpha$ | $\frac{2}{3}-\alpha$ | $-\frac{1}{3}$ | $-\frac{1}{3}$ | $1$ | $1$ | $1$ | $1$ | $1$ | $\bar{3}$ | $3$ | $1$ | $3$ | $3$
N4-1-1-1-1-1-2 | $\bar{3}$ | $\bar{6}$ | $1$ | $3$ | $3$
N4-1-1-1-1-1-3 | $3$ | $8$ | $3$ | $3$ | $3$
N4-1-1-1-1-1-4 | $\bar{6}$ | $8$ | $3$ | $3$ | $3$
N4-1-1-1-1-1-5 | $1$ | $\bar{3}$ | $\bar{3}$ | $3$ | $3$
N4-1-1-1-1-1-6 | $8$ | $\bar{3}$ | $\bar{3}$ | $3$ | $3$
N4-1-1-1-1-1-7 | $8$ | $6$ | $\bar{3}$ | $3$ | $3$
N4-1-1-1-1-1-8 | $3$ | $8$ | $\bar{6}$ | $3$ | $3$
N4-1-1-1-1-1-9 | $8$ | $6$ | $6$ | $3$ | $3$
N4-1-1-1-1-2-1 | $2$ | $2$ | $2$ | $1$ | $1$ | $\bar{3}$ | $3$ | $1$ | $3$ | $3$
N4-1-1-1-1-2-2 | $\bar{3}$ | $\bar{6}$ | $1$ | $3$ | $3$
N4-1-1-1-1-2-3 | $3$ | $8$ | $3$ | $3$ | $3$
N4-1-1-1-1-2-4 | $\bar{6}$ | $8$ | $3$ | $3$ | $3$
N4-1-1-1-1-2-5 | $1$ | $\bar{3}$ | $\bar{3}$ | $3$ | $3$
N4-1-1-1-1-2-6 | $8$ | $\bar{3}$ | $\bar{3}$ | $3$ | $3$
N4-1-1-1-1-2-7 | $8$ | $6$ | $\bar{3}$ | $3$ | $3$
N4-1-1-1-1-2-8 | $3$ | $8$ | $\bar{6}$ | $3$ | $3$
N4-1-1-1-1-2-9 | $8$ | $6$ | $6$ | $3$ | $3$
N4-1-1-1-1-3-1 | $3$ | $3$ | $3$ | $1$ | $1$ | $\bar{3}$ | $3$ | $1$ | $3$ | $3$
N4-1-1-1-1-3-2 | $\bar{3}$ | $\bar{6}$ | $1$ | $3$ | $3$
N4-1-1-1-1-3-3 | $3$ | $8$ | $3$ | $3$ | $3$
N4-1-1-1-1-3-4 | $\bar{6}$ | $8$ | $3$ | $3$ | $3$
N4-1-1-1-1-3-5 | $1$ | $\bar{3}$ | $\bar{3}$ | $3$ | $3$
N4-1-1-1-1-3-6 | $8$ | $\bar{3}$ | $\bar{3}$ | $3$ | $3$
N4-1-1-1-1-3-7 | $8$ | $6$ | $\bar{3}$ | $3$ | $3$
N4-1-1-1-1-3-8 | $3$ | $8$ | $\bar{6}$ | $3$ | $3$
N4-1-1-1-1-3-9 | $8$ | $6$ | $6$ | $3$ | $3$
Table 6:  The one-loop models for the $0\nu\beta\beta$ decay operator N4 arising from the diagram N4-1-1-1-1 shown in figure 4.  Model Name | $U(1)_{Y}$ | $SU(2)_{L}$ | $SU(3)_{C}$
---|---|---|---
$F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$ | $F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$ | $F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$
N4-1-1-1-2-1-1 | $-\frac{2}{3}+\alpha$ | $\alpha$ | $-\frac{1}{3}-\alpha$ | $-\frac{1}{3}$ | $\frac{2}{3}$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$ | $3$ | $1$ | $3$ | $\bar{6}$
N4-1-1-1-2-1-2 | $\bar{3}$ | $1$ | $3$ | $3$ | $\bar{6}$
N4-1-1-1-2-1-3 | $\bar{3}$ | $8$ | $3$ | $3$ | $\bar{6}$
N4-1-1-1-2-1-4 | $3$ | $\bar{3}$ | $\bar{3}$ | $3$ | $\bar{6}$
N4-1-1-1-2-1-5 | $3$ | $6$ | $\bar{3}$ | $3$ | $\bar{6}$
N4-1-1-1-2-1-6 | $6$ | $8$ | $\bar{6}$ | $3$ | $\bar{6}$
N4-1-1-1-2-1-7 | $\bar{6}$ | $\bar{3}$ | $6$ | $3$ | $\bar{6}$
N4-1-1-1-2-1-8 | $8$ | $3$ | $8$ | $3$ | $\bar{6}$
N4-1-1-1-2-1-9 | $8$ | $\bar{6}$ | $8$ | $3$ | $\bar{6}$
N4-1-1-1-2-2-1 | $2$ | $2$ | $2$ | $1$ | $1$ | $1$ | $3$ | $1$ | $3$ | $\bar{6}$
N4-1-1-1-2-2-2 | $\bar{3}$ | $1$ | $3$ | $3$ | $\bar{6}$
N4-1-1-1-2-2-3 | $\bar{3}$ | $8$ | $3$ | $3$ | $\bar{6}$
N4-1-1-1-2-2-4 | $3$ | $\bar{3}$ | $\bar{3}$ | $3$ | $\bar{6}$
N4-1-1-1-2-2-5 | $3$ | $6$ | $\bar{3}$ | $3$ | $\bar{6}$
N4-1-1-1-2-2-6 | $6$ | $8$ | $\bar{6}$ | $3$ | $\bar{6}$
N4-1-1-1-2-2-7 | $\bar{6}$ | $\bar{3}$ | $6$ | $3$ | $\bar{6}$
N4-1-1-1-2-2-8 | $8$ | $3$ | $8$ | $3$ | $\bar{6}$
N4-1-1-1-2-2-9 | $8$ | $\bar{6}$ | $8$ | $3$ | $\bar{6}$
N4-1-1-1-2-3-1 | $3$ | $3$ | $3$ | $1$ | $1$ | $1$ | $3$ | $1$ | $3$ | $\bar{6}$
N4-1-1-1-2-3-2 | $\bar{3}$ | $1$ | $3$ | $3$ | $\bar{6}$
N4-1-1-1-2-3-3 | $\bar{3}$ | $8$ | $3$ | $3$ | $\bar{6}$
N4-1-1-1-2-3-4 | $3$ | $\bar{3}$ | $\bar{3}$ | $3$ | $\bar{6}$
N4-1-1-1-2-3-5 | $3$ | $6$ | $\bar{3}$ | $3$ | $\bar{6}$
N4-1-1-1-2-3-6 | $6$ | $8$ | $\bar{6}$ | $3$ | $\bar{6}$
N4-1-1-1-2-3-7 | $\bar{6}$ | $\bar{3}$ | $6$ | $3$ | $\bar{6}$
N4-1-1-1-2-3-8 | $8$ | $3$ | $8$ | $3$ | $\bar{6}$
N4-1-1-1-2-3-9 | $8$ | $\bar{6}$ | $8$ | $3$ | $\bar{6}$
Table 7:  The one-loop models for the $0\nu\beta\beta$ decay operator N4 arising from the diagram N4-1-1-1-2 shown in figure 4. Model Name | $U(1)_{Y}$ | $SU(2)_{L}$ | $SU(3)_{C}$
---|---|---|---
$F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$ | $F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$ | $F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$
N4-1-1-1-3-1-1 | $1+\alpha$ | $\alpha$ | $-2-\alpha$ | $\frac{4}{3}$ | $\frac{2}{3}$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$ | $6$ | $\bar{6}$
N4-1-1-1-3-1-2 | $\bar{3}$ | $\bar{3}$ | $3$ | $6$ | $\bar{6}$
N4-1-1-1-3-1-3 | $6$ | $6$ | $\bar{6}$ | $6$ | $\bar{6}$
N4-1-1-1-3-1-4 | $8$ | $8$ | $8$ | $6$ | $\bar{6}$
N4-1-1-1-3-2-1 | $2$ | $2$ | $2$ | $1$ | $1$ | $1$ | $1$ | $1$ | $6$ | $\bar{6}$
N4-1-1-1-3-2-2 | $\bar{3}$ | $\bar{3}$ | $3$ | $6$ | $\bar{6}$
N4-1-1-1-3-2-3 | $6$ | $6$ | $\bar{6}$ | $6$ | $\bar{6}$
N4-1-1-1-3-2-4 | $8$ | $8$ | $8$ | $6$ | $\bar{6}$
N4-1-1-1-3-3-1 | $3$ | $3$ | $3$ | $1$ | $1$ | $1$ | $1$ | $1$ | $6$ | $\bar{6}$
N4-1-1-1-3-3-2 | $\bar{3}$ | $\bar{3}$ | $3$ | $6$ | $\bar{6}$
N4-1-1-1-3-3-3 | $6$ | $6$ | $\bar{6}$ | $6$ | $\bar{6}$
N4-1-1-1-3-3-4 | $8$ | $8$ | $8$ | $6$ | $\bar{6}$
Table 8:  The one-loop models for the $0\nu\beta\beta$ decay operator N4 arising from the diagram N4-1-1-1-3 shown in figure 4. Model Name | $U(1)_{Y}$ | $SU(2)_{L}$ | $SU(3)_{C}$
---|---|---|---
$F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$ | $F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$ | $F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$
N4-1-1-1-4-1-1 | $-\frac{1}{3}+\alpha$ | $\alpha$ | $\frac{2}{3}-\alpha$ | $-2$ | $\frac{4}{3}$ | $1$ | $1$ | $1$ | $1$ | $1$ | $\bar{3}$ | $\bar{6}$ | $1$ | $1$ | $6$
N4-1-1-1-4-1-2 | $3$ | $8$ | $3$ | $1$ | $6$
N4-1-1-1-4-1-3 | $\bar{6}$ | $8$ | $3$ | $1$ | $6$
N4-1-1-1-4-1-4 | $1$ | $\bar{3}$ | $\bar{3}$ | $1$ | $6$
N4-1-1-1-4-1-5 | $8$ | $\bar{3}$ | $\bar{3}$ | $1$ | $6$
N4-1-1-1-4-1-6 | $3$ | $8$ | $\bar{6}$ | $1$ | $6$
N4-1-1-1-4-1-7 | $8$ | $6$ | $6$ | $1$ | $6$
N4-1-1-1-4-2-1 | $2$ | $2$ | $2$ | $1$ | $1$ | $\bar{3}$ | $\bar{6}$ | $1$ | $1$ | $6$
N4-1-1-1-4-2-2 | $3$ | $8$ | $3$ | $1$ | $6$
N4-1-1-1-4-2-3 | $\bar{6}$ | $8$ | $3$ | $1$ | $6$
N4-1-1-1-4-2-4 | $1$ | $\bar{3}$ | $\bar{3}$ | $1$ | $6$
N4-1-1-1-4-2-5 | $8$ | $\bar{3}$ | $\bar{3}$ | $1$ | $6$
N4-1-1-1-4-2-6 | $3$ | $8$ | $\bar{6}$ | $1$ | $6$
N4-1-1-1-4-2-7 | $8$ | $6$ | $6$ | $1$ | $6$
N4-1-1-1-4-3-1 | $3$ | $3$ | $3$ | $1$ | $1$ | $\bar{3}$ | $\bar{6}$ | $1$ | $1$ | $6$
N4-1-1-1-4-3-2 | $3$ | $8$ | $3$ | $1$ | $6$
N4-1-1-1-4-3-3 | $\bar{6}$ | $8$ | $3$ | $1$ | $6$
N4-1-1-1-4-3-4 | $1$ | $\bar{3}$ | $\bar{3}$ | $1$ | $6$
N4-1-1-1-4-3-5 | $8$ | $\bar{3}$ | $\bar{3}$ | $1$ | $6$
N4-1-1-1-4-3-6 | $3$ | $8$ | $\bar{6}$ | $1$ | $6$
N4-1-1-1-4-3-7 | $8$ | $6$ | $6$ | $1$ | $6$
Table 9:  The one-loop models for the $0\nu\beta\beta$ decay operator N4 arising from the diagram N4-1-1-1-4 shown in figure 4. Model Name | $U(1)_{Y}$ | $SU(2)_{L}$ | $SU(3)_{C}$
---|---|---|---
$F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$ | $F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$ | $F_{1}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$
N4-1-1-1-5-1-1 | $-\frac{2}{3}+\alpha$ | $\alpha$ | $\frac{4}{3}-\alpha$ | $-2$ | $\frac{2}{3}$ | $1$ | $1$ | $1$ | $1$ | $1$ | $3$ | $6$ | $1$ | $1$ | $\bar{6}$
N4-1-1-1-5-1-2 | $1$ | $3$ | $3$ | $1$ | $\bar{6}$
N4-1-1-1-5-1-3 | $8$ | $3$ | $3$ | $1$ | $\bar{6}$
N4-1-1-1-5-1-4 | $\bar{3}$ | $8$ | $\bar{3}$ | $1$ | $\bar{6}$
N4-1-1-1-5-1-5 | $6$ | $8$ | $\bar{3}$ | $1$ | $\bar{6}$
N4-1-1-1-5-1-6 | $8$ | $\bar{6}$ | $\bar{6}$ | $1$ | $\bar{6}$
N4-1-1-1-5-1-7 | $\bar{3}$ | $8$ | $6$ | $1$ | $\bar{6}$
N4-1-1-1-5-2-1 | $2$ | $2$ | $2$ | $1$ | $1$ | $3$ | $6$ | $1$ | $1$ | $\bar{6}$
N4-1-1-1-5-2-2 | $1$ | $3$ | $3$ | $1$ | $\bar{6}$
N4-1-1-1-5-2-3 | $8$ | $3$ | $3$ | $1$ | $\bar{6}$
N4-1-1-1-5-2-4 | $\bar{3}$ | $8$ | $\bar{3}$ | $1$ | $\bar{6}$
N4-1-1-1-5-2-5 | $6$ | $8$ | $\bar{3}$ | $1$ | $\bar{6}$
N4-1-1-1-5-2-6 | $8$ | $\bar{6}$ | $\bar{6}$ | $1$ | $\bar{6}$
N4-1-1-1-5-2-7 | $\bar{3}$ | $8$ | $6$ | $1$ | $\bar{6}$
N4-1-1-1-5-3-1 | $3$ | $3$ | $3$ | $1$ | $1$ | $3$ | $6$ | $1$ | $1$ | $\bar{6}$
N4-1-1-1-5-3-2 | $1$ | $3$ | $3$ | $1$ | $\bar{6}$
N4-1-1-1-5-3-3 | $8$ | $3$ | $3$ | $1$ | $\bar{6}$
N4-1-1-1-5-3-4 | $\bar{3}$ | $8$ | $\bar{3}$ | $1$ | $\bar{6}$
N4-1-1-1-5-3-5 | $6$ | $8$ | $\bar{3}$ | $1$ | $\bar{6}$
N4-1-1-1-5-3-6 | $8$ | $\bar{6}$ | $\bar{6}$ | $1$ | $\bar{6}$
N4-1-1-1-5-3-7 | $\bar{3}$ | $8$ | $6$ | $1$ | $\bar{6}$
Table 10: The one-loop models for the $0\nu\beta\beta$ decay operator N4
arising from the diagram N4-1-1-1-5 shown in figure 4.

## References

  * [1] Particle Data Group Collaboration, P. A. Zyla et al., “Review of Particle Physics,” PTEP 2020 no. 8, (2020) 083C01.
  * [2] J. Schechter and J. W. F. Valle, “Neutrinoless Double beta Decay in SU(2) x U(1) Theories,” Phys. Rev. D 25 (1982) 2951.
  * [3] KamLAND-Zen Collaboration, A. Gando et al., “Search for Majorana Neutrinos near the Inverted Mass Hierarchy Region with KamLAND-Zen,” Phys. Rev. Lett. 117 no. 8, (2016) 082503, arXiv:1605.02889 [hep-ex]. [Addendum: Phys.Rev.Lett. 117, 109903 (2016)].
  * [4] LEGEND Collaboration, N. Abgrall et al., “The Large Enriched Germanium Experiment for Neutrinoless $\beta\beta$ Decay: LEGEND-1000 Preconceptual Design Report,” arXiv:2107.11462 [physics.ins-det].
  * [5] nEXO Collaboration, G. Adhikari et al., “nEXO: Neutrinoless double beta decay search beyond $10^{28}$ year half-life sensitivity,” arXiv:2106.16243 [nucl-ex].
  * [6] H. Päs, M. Hirsch, H. V. Klapdor-Kleingrothaus, and S. G. Kovalenko, “A superformula for neutrinoless double beta decay II: The short range part,” Physics Letters B 498 no. 1-2, (Jan., 2001) 35–39, arXiv:hep-ph/0008182.
  * [7] H. Päs, M. Hirsch, H. V. Klapdor-Kleingrothaus, and S. G. Kovalenko, “Towards a superformula for neutrinoless double beta decay,” Physics Letters B 453 no. 3-4, (May, 1999) 194–198, arXiv:hep-ph/9804374.
  * [8] S. R. Elliott and P. Vogel, “Double beta decay,” Ann. Rev. Nucl. Part. Sci. 52 (2002) 115–151, arXiv:hep-ph/0202264.
  * [9] W. Rodejohann, “Neutrino-less Double Beta Decay and Particle Physics,” International Journal of Modern Physics E 20 no. 09, (Sept., 2011) 1833–1930, arXiv:1106.1334.
  * [10] F. F. Deppisch, M. Hirsch, and H. Päs, “Neutrinoless Double Beta Decay and Physics Beyond the Standard Model,” Journal of Physics G: Nuclear and Particle Physics 39 no. 12, (Dec., 2012) 124007, arXiv:1208.0727.
  * [11] S. Dell’Oro, S. Marcocci, M. Viel, and F. Vissani, “Neutrinoless double beta decay: 2015 review,” Adv. High Energy Phys. 2016 (2016) 2162659, arXiv:1601.07512 [hep-ph].
  * [12] M. J. Dolinski, A. W. P. Poon, and W. Rodejohann, “Neutrinoless Double-Beta Decay: Status and Prospects,” Ann. Rev. Nucl. Part. Sci. 69 (2019) 219–251, arXiv:1902.04097 [nucl-ex].
  * [13] M. L. Graesser, “An electroweak basis for neutrinoless double $\beta$ decay,” Journal of High Energy Physics 2017 no. 8, (Aug., 2017) 99, arXiv:1606.04549.
  * [14] Y. Liao and X.-D. Ma, “An explicit construction of the dimension-9 operator basis in the standard model effective field theory,” JHEP 11 (2020) 152, arXiv:2007.08125 [hep-ph].
  * [15] H.-L. Li, Z. Ren, M.-L. Xiao, J.-H. Yu, and Y.-H. Zheng, “Complete set of dimension-nine operators in the standard model effective field theory,” Phys. Rev. D 104 no. 1, (2021) 015025, arXiv:2007.07899 [hep-ph].
  * [16] F. Bonnet, M. Hirsch, T. Ota, and W. Winter, “Systematic decomposition of the neutrinoless double beta decay operator,” Journal of High Energy Physics 2013 no. 3, (Mar., 2013) 55, arXiv:1212.3045.
  * [17] J. C. Helo, M. Hirsch, and T. Ota, “Long-range contributions to double beta decay revisited,” JHEP 06 (2016) 006, arXiv:1602.03362 [hep-ph].
  * [18] M. Gustafsson, J. M. No, and M. A. Rivera, “Radiative neutrino mass generation linked to neutrino mixing and 0$\nu$$\beta$$\beta$-decay predictions,” Phys. Rev. D 90 no. 1, (2014) 013012, arXiv:1402.0515 [hep-ph].
  * [19] Z. Liu and P.-H. Gu, “Extending two Higgs doublet models for two-loop neutrino mass generation and one-loop neutrinoless double beta decay,” Nucl. Phys. B 915 (2017) 206–223, arXiv:1611.02094 [hep-ph].
  * [20] J. Alcaide, D. Das, and A. Santamaria, “A model of neutrino mass and dark matter with large neutrinoless double beta decay,” JHEP 04 (2017) 049, arXiv:1701.01402 [hep-ph].
  * [21] W. Rodejohann and X.-J. Xu, “Loop-enhanced rate of neutrinoless double beta decay,” JHEP 11 (2019) 029, arXiv:1907.12478 [hep-ph].
  * [22] S. Antusch, J. P. Baumann, and E. Fernandez-Martinez, “Non-Standard Neutrino Interactions with Matter from Physics Beyond the Standard Model,” Nucl. Phys. B 810 (2009) 369–388, arXiv:0807.1003 [hep-ph].
  * [23] M. B. Gavela, D. Hernandez, T. Ota, and W. Winter, “Large gauge invariant non-standard neutrino interactions,” Phys. Rev. D 79 (2009) 013007, arXiv:0809.3451 [hep-ph].
  * [24] F. Bonnet, M. Hirsch, T. Ota, and W. Winter, “Systematic study of the d=5 Weinberg operator at one-loop order,” JHEP 07 (2012) 153, arXiv:1204.5862 [hep-ph].
  * [25] D. Aristizabal Sierra, A. Degee, L. Dorame, and M. Hirsch, “Systematic classification of two-loop realizations of the Weinberg operator,” JHEP 03 (2015) 040, arXiv:1411.7038 [hep-ph].
  * [26] R. Cepedello, R. M. Fonseca, and M. Hirsch, “Systematic classification of three-loop realizations of the Weinberg operator,” JHEP 10 (2018) 197, arXiv:1807.00629 [hep-ph]. [Erratum: JHEP 06, 034 (2019)].
  * [27] C.-Y. Yao and G.-J. Ding, “Systematic Study of One-Loop Dirac Neutrino Masses and Viable Dark Matter Candidates,” Phys. Rev. D 96 no. 9, (2017) 095004, arXiv:1707.09786 [hep-ph]. [Erratum: Phys.Rev.D 98, 039901 (2018)].
  * [28] C.-Y. Yao and G.-J. Ding, “Systematic analysis of Dirac neutrino masses from a dimension five operator,” Phys. Rev. D 97 no. 9, (2018) 095042, arXiv:1802.05231 [hep-ph].
  * [29] S. Centelles Chuliá, R. Cepedello, E. Peinado, and R. Srivastava, “Systematic classification of two loop $d$ = 4 Dirac neutrino mass models and the Diracness-dark matter stability connection,” JHEP 10 (2019) 093, arXiv:1907.08630 [hep-ph].
  * [30] S. Jana, P. Vishnu, and S. Saad, “Minimal realizations of Dirac neutrino mass from generic one-loop and two-loop topologies at $d=5$,” JCAP 04 (2020) 018, arXiv:1910.09537 [hep-ph].
  * [31] R. M. Fonseca, “GroupMath: A Mathematica package for group theory calculations,” arXiv:2011.01764 [hep-th].
  * [32] P.-T. Chen, G.-J. Ding, and C.-Y. Yao, Supplementary material: Decomposition of $d=9$ short-range $0\nu\beta\beta$ decay operators at one-loop level, 2021 (accessed Octomber 29, 2021).  http://staff.ustc.edu.cn/~dinggj/supplementary_materials/0nbb.zip.
  * [33] M. Duerr, M. Lindner, and A. Merle, “On the Quantitative Impact of the Schechter-Valle Theorem,” JHEP 06 (2011) 091, arXiv:1105.0901 [hep-ph].
  * [34] J. C. Helo, M. Hirsch, T. Ota, and F. A. P. dos Santos, “Double beta decay and neutrino mass models,” Journal of High Energy Physics 2015 no. 5, (May, 2015) 92, arXiv:1502.05188.
  * [35] GERDA Collaboration, M. Agostini et al., “Final Results of GERDA on the Search for Neutrinoless Double-$\beta$ Decay,” Phys. Rev. Lett. 125 no. 25, (2020) 252502, arXiv:2009.06079 [nucl-ex].
  * [36] 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].
  * [37] J. van der Bij and M. J. G. Veltman, “Two Loop Large Higgs Mass Correction to the rho Parameter,” Nucl. Phys. B 231 (1984) 205–234.
  * [38] K. L. McDonald and B. H. J. McKellar, “Evaluating the two loop diagram responsible for neutrino mass in Babu’s model,” arXiv:hep-ph/0309270.
  * [39] P. W. Angel, Y. Cai, N. L. Rodd, M. A. Schmidt, and R. R. Volkas, “Testable two-loop radiative neutrino mass model based on an $LLQd^{c}Qd^{c}$ effective operator,” JHEP 10 (2013) 118, arXiv:1308.0463 [hep-ph]. [Erratum: JHEP 11, 092 (2014)].
  * [40] A. Freitas, “Three-loop vacuum integrals with arbitrary masses,” JHEP 11 (2016) 145, arXiv:1609.09159 [hep-ph].
  * [41] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [42] I. Esteban, M. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, “The fate of hints: updated global analysis of three-flavor neutrino oscillations,” JHEP 09 (2020) 178, arXiv:2007.14792 [hep-ph].