# Nuclear $0\nu 2\beta$ decays in $B-L$ symmetric SUSY model and in TeV scale
left-right symmetric model

Jin-Lei Yang1,2<EMAIL_ADDRESS>Chao-Hsi Chang1,2<EMAIL_ADDRESS>Tai-
Fu Feng1,3<EMAIL_ADDRESS>1CAS Key Laboratory of Theoretical Physics,
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190,
China
2School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing 100049, China
3Department of Physics, Hebei University, Key Laboratory of High-precision
Computation and Application of Quantum Field Theory of Hebei Province,
Baoding, 071002, China

###### Abstract

In this paper we take B-L supersymmetric standard model (B-LSSM) and TeV scale
left-right symmetric model (LRSM) as two representations of the two kinds of
new physics models to study the nuclear neutrinoless double beta decays ($0\nu
2\beta$) so as to see the senses onto these two kinds of models when the
decays are taken into account additionally. Within the parameter spaces
allowed by all the exist experimental data, the decay half-life of the nucleus
76Ge and 136Xe, $T^{0\nu}_{1/2}$(76Ge, 136Xe), is precisely calculated and the
results are presented properly. Based on the numerical results, we conclude
that the room of LRSM type models for the foreseeable future experimental
observations on the decays is greater than that of B-LSSM type models.

neutrinoless double beta decay, QCD correction, B-LSSM, LRSM

##  I Introduction

Tiny but nonzero neutrino masses explaining neutrino oscillation experiments
Tanabashi:2018oca  are unambiguous evidences of new physics (NP) beyond the
standard model (SM). It is because that in SM there are only left-handed
neutrinos, so that the neutrinos can acquire neither Dirac masses nor Majorana
masses. Hence to explore any mechanism inducing the tiny neutrino masses, as
well as relevant phenomenology, is an important direction to search for NP.
The simple extension of SM is to introduce three right-handed neutrinos in
singlet of the gauge group SU(2) additionally, where the neutrinos acquire
Dirac masses, and to fit the neutrino oscillation and nuclear decay
experiments as well as astronomy observations, the corresponding Yukawa
couplings of Higgs to the neutrinos are requested so tiny as $\lesssim
10^{-12}$, that is quite unnatural.

However neutrino(s) may acquire masses naturally by introducing Majorana mass
terms in extended SMs. Once a Majorana mass term is introduced, certain
interesting physics arise. One of the consequences is that the lepton-number
violation (LNV) processes, e.g. the nuclear neutrinoless double beta decays
($0\nu 2\beta$) etc may occur. Of them, the $0\nu 2\beta$ decays are specially
interesting, because they may tell us sensitively some on the nature, whether
Dirac Dirac:1925jy  or Majorana Majorana:1937vz , of the neutrinos. When the
decays $0\nu 2\beta$ are observed in experiments, most likely, the neutrinos
contain Majorana components. Thus studying $0\nu 2\beta$ decays is attracting
special attentions.

Nowadays there are several experiments running to observe the $0\nu 2\beta$
decays, and the most stringent experimental bounds on the processes are
obtained by GERDA Agostini:2017iyd ; Guiseppe:2017rcr  and KamLAND-Zen
KamLAND-Zen:2016pfg ; Albert:2017owj . They adopt suitable approaches and
nucleus such as 76Ge and 136Xe respectively. Now the latest experimental lower
bound on the decay half-life given by GERDA experiments is
$T_{1/2}^{0\nu}>1.8\times 10^{26}$ years ($90\%$ C.L.) for nucleus 76Ge
Agostini:2020xta , and in the near future the sensitivity can reach to
$10^{28}$ years talk . For the nucleus 136Xe, the most stringent lower bound
on the decay half-life is $T_{1/2}^{0\nu}>1.07\times 10^{26}$ years ($90\%$
C.L.) given by KamLAND-Zen KamLAND-Zen:2016pfg , and the corresponding future
sensitivity can reach to $2.4\times 10^{27}$ years Agostini:2020adk .
Moreover, underground experiments PANDAX, CDEX etc, which are originally
designed for searching for WIMP dark matter, are also planning to seek the
$0\nu 2\beta$ decays i.e. they may observe the $0\nu 2\beta$ decays with the
sensitivity which at least will set a fresh lower bound.

In literature, there are a lot of theoretical analyses on the $0\nu 2\beta$
decays. The analyses are carried out generally by dividing the estimation of
the $0\nu 2\beta$ decays into three ‘factors’: one is at quark level to
evaluate the amplitude for the ‘core’ process $d+d\to u+u+e+e$ of the decays;
the second one is, from quark level to nucleon level, to involve the quark
process into the relevant nucleon one i.e. the ‘initial’ quarks $d$ and $d$
involve into the two neutrons in the initial nucleus and the ‘final’ quarks
$u$ and $u$ involve into the two protons in the final nucleus; the third one
is, from nucleon level to nucleus level, the relevant nucleons involve in the
initial nucleus and the final nucleus properly. For the ‘core’ process, in
Ref. Pas:1999fc  a general Lorentz-invariant effective Lagrangian is
constructed by dimension-9 operators, and in Ref. Gonzalez:2015ady ;
Liao:2020roy  the QCD corrections to all of these dimension-9 operators are
calculated. In Ref. Pas:2000vn  the short-range effects at quark level are
considered, the analyses of the decay rates in the SM effective field theory
are presented in Refs. Cirigliano:2017djv ; Cirigliano:2018hja ;
Cirigliano:2018yza , the $0\nu 2\beta$ decay rates are derived in Refs.
Dekens:2020ttz , the corresponding nuclear matrix elements (NME) and phase-
space factors (PSF) for the second and the third factors of the decays i.e.
from quark level to nucleon and nucleus levels, are considered in Refs.
Barea:2015kwa ; Suhonen:1991sk ; Caurier:2007xz ; Simkovic:2007vu ;
Menendez:2008jp ; Barea:2009zza ; Rodriguez:2010mn ; Suhonen:2012wd ;
Barea:2013bz ; Graf:2018ozy ; Deppisch:2020ztt , some theoretical predictions
on the $0\nu 2\beta$ for certain models are presented in Refs. Deppisch:2012nb
; LopezPavon:2012zg ; Awasthi:2013ff ; Asaka:2016zib ; Pascoli:2013fiz , and
the theoretical analyses on the decays are reviewed in Refs. DellOro:2016tmg ;
Engel:2016xgb ; Dolinski:2019nrj .

In this work, we are investigating the constrains from the $0\nu 2\beta$
decays for the $B-L$ supersymmetric model (B-LSSM) and for the TeV scale left-
right symmetric model (LRSM) comparatively. It is because that the two models
are typical: both have LNV source but the mechanisms which give rise to the
Majorana mass terms are different Dev:2016dja ; Patra:2015bga ; Mitra:2016kov
; Maiezza:2016ybz ; Khalil:2008ps ; Elsayed:2011de ; Khalil:2015naa ;
DelleRose:2017ukx ; Yang:2018fvw ; Yang:2018utw ; Yang:2018guw ; Yang:2019aao
; Yang:2020ebs ; Yang:2020bmh . In the B-LSSM, the tiny neutrino masses are
acquired naturally through the so-called type-I seesaw mechanism which is
proposed firstly by Weinberg Weinberg:1979sa . In the LRSM Mohapatra:2006gs ;
Mohapatra:2005wg , the tiny neutrino masses are acquired by both of type-I and
type-II seesaw mechanisms, in addition, the new right-handed gauge bosons
$W^{\pm}_{R}$ is introduced in this model, then both of left-handed and right-
handed currents cause the $0\nu 2\beta$ decays Mohapatra:1979ia ;
Mohapatra:1981pm ; Picciotto:1982qe ; Hirsch:1996qw ; Tello:2010am ;
Chakrabortty:2012mh ; Barry:2013xxa ; Dev:2013vxa ; Dev:2013oxa ;
Huang:2013kma ; Dev:2014xea ; Ahmed:2019vum ; Borah:2015ufa ; Ge:2015yqa ;
Awasthi:2015ota ; Bambhaniya:2015ipg ; Pritimita:2016fgr ; Deppisch:2017vne .
As results, the computations of the decays are much more complicated in the
LRSM than those in the B-LSSM. Hence these two models, being representatives
of NP models, are typical for the $0\nu 2\beta$ decays, and one may learn the
mechanisms in the models well via analyzing the $0\nu 2\beta$ decays
comparatively.

In the study here, we will mainly focus on the first ‘factor’ about the quark
level i.e. the core process, which relates to the applied specific model
closely. We will evaluate the Wilson coefficients of the operators relevant to
the core process $d+d\to u+u+e+e$ etc on the models, whereas the estimation of
the other two ‘factors’, i.e. to evaluate ‘NME’ and ‘PSF’ etc, we will follow
the literatures Graf:2018ozy ; Deppisch:2012nb ; LopezPavon:2012zg ;
Awasthi:2013ff ; Asaka:2016zib . Respect to the ‘core’ process $d+d\to
u+u+e+e$, all of the contributions in the B-LSSM can be deduced directly quite
well, while the contributions in LRSM cannot be so. As shown in Ref.
Kotila:2021xgw , the calculations in the LRSM are much complicated and the
interference effects are quite hard to be considered well. In this work, a new
approximation, i.e. the momenta of the two involved quarks inside the initial
or final nuclei is tried to be set equal, is made so as may reduce all
contributions in the LRSM quite similar to the case of B-LSSM. Then the
calculations in LRSM are simplified quite a lot and under the approximation
the interference effects can be treated comparatively well. Finally, for
comparison, we also present the results obtained by the traditional method
Dev:2014xea .

The paper is organized as follows: In Sec. II, for B-LSSM, the seesaw
mechanisms which give rise to the tiny neutrino masses, the heavy neutral
leptons as well, the relevant interactions etc, the calculations of the $0\nu
2\beta$ decay half-lifes of the nuclei are given. Similarly, in Sec. III, for
LRSM, the seesaw mechanisms which give rise to the tiny neutrino masses, as
well the heavy neutral leptons, the relevant interactions and the calculations
of the $0\nu 2\beta$ decay half-lifes of the nuclei are given. In Sec. IV.1
and IV.2 the numerical results for B-LSSM and LRSM are presented respectively.
Finally, in Sec. V brief discussions and conclusions are given. In the
appendix, the needed QCD corrections to the effective Lagrangian which
contains the dimension-9 operators in the region from the energy scale
$\mu\simeq M_{W}$ to the energy scale $\mu\simeq 1.0\;$GeV are collected.

##  II The B-LSSM for $0\nu 2\beta$ decays

In the B-LSSM, the local gauge group is $SU(3)_{C}\bigotimes
SU(2)_{L}\bigotimes U(1)_{Y}\bigotimes U(1)_{B-L}$, where $B$, $L$ denote the
baryon number and lepton number respectively, and the details about the gauge
fields, their breaking, extra lepton and Higgs fields etc can be found in
Refs. Khalil:2008ps ; Elsayed:2011de ; Khalil:2015naa ; DelleRose:2017ukx ;
Yang:2018fvw ; Yang:2018utw ; Yang:2018guw ; Yang:2019aao ; Yang:2020ebs ;
Yang:2020bmh . In this model, tiny neutrino masses are acquired by the so-
called type-I seesaw mechanism, when the $U(1)_{B-L}$ symmetry is broken
spontaneously by the two $U(1)_{B-L}$ singlet scalars (Higgs). The mass matrix
for neutrinos and neutral heavy leptons in the model can be expressed as

$\displaystyle\left(\begin{array}[]{c}0,\;\;\;\;\;\;\;M_{D}^{\rm T}\\\
M_{D},\;\;\;\;M_{R}\end{array}\right),$ (3)

and the mass matrix can be diagonalized in terms of a unitary matrix $U_{\nu}$
as follows:

$\displaystyle U_{\nu}^{\rm
T}\left(\begin{array}[]{c}0,\;\;\;\;\;\;\;M_{D}^{\rm T}\\\
M_{D},\;\;\;\;M_{R}\end{array}\right)U_{\nu}=\left(\begin{array}[]{c}\hat{m}_{\nu},\;\;\;\;0\\\
0,\;\;\;\;\hat{M}_{N}\end{array}\right),$ (8)

where the neutrino masses $\hat{m}_{\nu}={\rm
diag}(m_{\nu_{1}},m_{\nu_{2}},m_{\nu_{3}})$, the masses of the heavy neutral
leptons $\hat{M}_{N}={\rm diag}(M_{N_{1}},M_{N_{2}},M_{N_{3}})$ and $U_{\nu}$
is a matrix of $6\times 6$ which can be rewritten as

$\displaystyle U_{\nu}=\left(\begin{array}[]{c}U\quad S\\\ T\quad
V\end{array}\right),$ (11)

where $U,S,T,V$ are matrices of $3\times 3$.

The interactions, being applied later on, in the model are

$\displaystyle\mathcal{L}_{I}=\frac{{\rm
i}g_{2}}{\sqrt{2}}\sum_{j=1}^{3}\Big{[}U_{ij}\bar{e}_{i}\gamma^{\mu}P_{L}\nu_{j}W_{L,\mu}^{-}+S_{ij}\bar{e}_{i}\gamma^{\mu}P_{L}N_{j}W_{L,\mu}^{-}+h.c\Big{]},$
(12)

where $\nu,\;N$ are the four-component fermion fields of the light and heavy
neutral leptons respectively.

Figure 1: The dominant Feynman diagrams for the $0\nu 2\beta$ decays in the
B-LSSM Model. (a) The contributions from the heavy neutral lepton exchanges,
(b) The contributions from the light neutrino exchangs.

The Feynman diagrams which response for the dominated contributions to the
$0\nu 2\beta$ decays in the B-LSSM are plotted in Fig. 1. Note here that the
contributions from charged Higgs exchange(s) are ignored safely as they are
highly suppressed by the charged Higgs masses and Yukawa couplings, thus in
Fig. 1 the charged Higgs exchange do not appear at all. To evaluate the
contributions corresponding to the Feynman diagrams for the $0\nu 2\beta$
decays, and to consider the roles of the neutral leptons (light neutrinos and
heavy neutral leptons) in the Feynman diagrams, the useful formulae are
collected below so as to deal with the neutrino propagator sandwiched by
various chiral project operators $P_{L,R}=\frac{1}{2}(1\mp\gamma_{5})$:

$\displaystyle
P_{L}\frac{k\\!\\!\\!/+m}{k^{2}-m^{2}}P_{L}=\frac{m}{k^{2}-m^{2}}P_{L},\;\;P_{R}\frac{k\\!\\!\\!/+m}{k^{2}-m^{2}}P_{R}=\frac{m}{k^{2}-m^{2}}P_{R},$
$\displaystyle
P_{L}\frac{k\\!\\!\\!/+m}{k^{2}-m^{2}}P_{R}=\frac{k\\!\\!\\!/}{k^{2}-m^{2}}P_{R},\;\;P_{R}\frac{k\\!\\!\\!/+m}{k^{2}-m^{2}}P_{L}=\frac{k\\!\\!\\!/}{k^{2}-m^{2}}P_{L},$
(13)

and the propagator

$\displaystyle\frac{k\\!\\!\\!/+m}{k^{2}-m^{2}}\simeq\left\\{\begin{array}[]{l}-\frac{1}{m}\;,\;\;\;m^{2}\gg
k^{2}\\\ \frac{k\\!\\!\\!/+m}{k^{2}}\;,\;\;\;m^{2}\ll k^{2}\\\
\end{array}\right.,$ (16)

Relating to the exchanges of the heavy neutral leptons (the virtual neutral
lepton momentum $k$ has $|k|\simeq 0.10\;{\rm GeV}\ll M_{N_{i}}$) for the
decays, from Fig. 1 (a) the effective Lagrangian at the energy scale
$\mu\simeq M_{W_{L}}$ can be read out as

$\displaystyle\frac{2m_{p}}{G_{F}^{2}\cos\theta_{C}^{2}}\mathcal{L}_{\rm
eff}^{\rm
BL}(N)=\sum_{i}\frac{2m_{p}}{M_{N_{i}}}(S_{1i})^{2}[4(\bar{u}\gamma_{\mu}P_{L}d)(\bar{u}\gamma^{\mu}P_{L}d)\bar{e}P_{R}e^{c}]\equiv
C_{3R}^{LL}(N)\mathcal{O}_{3R}^{LL},$ $\displaystyle
C_{3R}^{LL}(N)=\sum_{i}\frac{2m_{p}}{M_{N_{i}}}(S_{1i})^{2},\;\mathcal{O}_{3Z}^{XY}=8(\bar{u}\gamma_{\mu}P_{X}d)(\bar{u}\gamma^{\mu}P_{Y}d)\bar{e}P_{Z}e^{c},$
(17)

where $X,Y,Z=L,R$, $\theta_{C}$ is the Cabibbo angle, $m_{p}$ is proton mass
introduced for normalization of the effective Lagrangian, and $S_{1i}$ is the
matrix elements in Eq. (11).

Since the nuclear $0\nu 2\beta$ decays take place at the energy scale about
$\mu\approx 0.10\;$GeV, obviously we need to consider the QCD corrections for
the effective Lagrangian obtained at the energy scale $\mu\simeq M_{W_{L}}$ in
Eq. (17) i.e. to evolute the effective Lagrangian in terms of renormalization
group equation (RGE) method from the energy scale $\mu\simeq M_{W_{L}}$ to
that $\mu\simeq 1.0\;$GeV first, where the corrections are in perturbative QCD
(pQCD) region. Thus, for completeness, the QCD corrections to all of the
possible dimension-9 operators which may contribute to the nuclear $0\nu
2\beta$ decays, are calculated by the RGE method, and the details of
computations are collected in the appendix. Whereas the QCD corrections in the
energy scale region $\mu\simeq 1.0\;$GeV$\sim\mu\simeq 0.10\;$GeV, being in
non-perturbative QCD region, we take into account them by inputting the
experimental measurements for the relevant current matrix elements of
nucleons, which emerge when calculating the amplitude based on the effective
Lagrangian at $\mu\simeq 0.10\;$GeV.

For the decays when considering contributions from the neutrino
($m_{\nu_{i}}\ll|k|$) exchanges as described by the Feynman diagram Fig. 1
(b), and the interferences between the light neutrinos’ and heavy neutral
leptons’ contributions, to derive the effective Lagrangian for the neutrino
contributions is better at the energy scale $\mu\simeq 1.0\;{\rm GeV}$ as
heavy neutral leptons too. At this energy scale the effective Lagrangian can
be written down according to the Feynman diagram Fig. 1 (b) as below:

$\displaystyle\frac{2m_{p}}{G_{F}^{2}\cos\theta_{C}^{2}}\mathcal{L}_{\rm
eff}^{\rm
BL}(\nu)=\frac{m_{\nu_{i}}}{m_{e}}(U_{1i})^{2}\cdot\frac{2m_{p}m_{e}}{-k^{2}}\mathcal{O}_{3R}^{LL}\equiv
C_{3R}^{LL}(\nu)\frac{2m_{p}m_{e}}{-k^{2}}\mathcal{O}_{3R}^{LL},$
$\displaystyle C_{3R}^{LL}(\nu)=\frac{m_{\nu_{i}}}{m_{e}}(U_{1i})^{2}.$ (18)

Since the light neutrino exchange is of long range, the QCD corrections in the
region $\mu\simeq 1.0\;$GeV$\sim\mu\simeq 0.10\;$GeV to the coefficients in
Eq. (18) may be involved just as that in the case done as in the above for the
heavy neutral lepton case, via inputting the experimental measurements for the
relevant current matrix elements of nucleons, which emerge when calculating
the matrix elements in the amplitudes at $\mu\approx 0.10\;$GeV.

To evaluate the half-life $T_{1/2}^{0\nu}\equiv\ln 2/\Gamma$ of the $0\nu
2\beta$ decays, the contributions from the heavy neutral leptons (Fig 1 (a))
and those from the neutrinos (Fig 1 (b)) should be summed up for the
amplitudes. The half-life $T_{1/2}^{0\nu}\equiv\ln 2/\Gamma$ of the $0\nu
2\beta$ decays can be written as Deppisch:2020ztt

$\displaystyle\frac{1}{T_{1/2}^{0\nu}}=G^{0\nu}|M^{0\nu}|^{2}\Big{|}\frac{m_{ee}^{\rm
BL}}{m_{e}}\Big{|}^{2},$ (19)

where $G^{0\nu}=2.36\times 10^{-15}\;(14.56\times 10^{-15})\;{\rm years}^{-1}$
Deppisch:2020ztt  for ${}^{76}{\rm Ge}\;(^{136}{\rm Xe})$ is the PSF,
$M^{0\nu}=-6.64\pm 1.06$ $(-3.60\pm 0.58)$ Barea:2015kwa ; Deppisch:2020ztt
for ${}^{76}{\rm Ge}\;(^{136}{\rm Xe})$111The NMEs adopted here are obtained
within the framework of the microscopic interacting meson model for nuclei
Deppisch:2020ztt , and the uncertainties for NMEs calculations from the
various nuclear structure models are quite wild, here we remind only that the
NMEs obtained by various approaches are varying by a factor of $(2\sim 3)$
roughly Kotila:2021gwk . is the NME corresponding to the long range
contributions which is defined as

$\displaystyle
M^{0\nu}\equiv\langle{\mathcal{O}}_{F}^{+}|\frac{2m_{p}m_{e}}{-k^{2}}[4(\bar{u}\gamma_{\mu}P_{L}d)(\bar{u}\gamma^{\mu}P_{L}d)]|\mathcal{O}_{I}^{+}\rangle$
(20)

with $|\mathcal{O}_{I}^{+}\rangle$, $\langle\mathcal{O}_{F}^{+}|$ denoting the
initial and final nuclear states respectively. Note that in Eq. (20) that the
factor $\frac{2m_{p}m_{e}}{-k^{2}}$ in Eq. (18) is absorbed into the so-called
‘neutrino potential’ which is used to compute the long range NME. And

$\displaystyle m_{ee}^{\rm BL}\equiv
U_{3}^{XX}C_{3R}^{LL}(N)\frac{M_{3}^{XX}(N)}{M^{0\nu}}+C_{3R}^{LL}(\nu),$ (21)

where $U_{3}^{XX}$ is the QCD running factor from $\mu\simeq M_{W_{L}}$ to
$\mu\simeq 1.0\;$GeV (the numerical result of $U_{3}^{XX}$ can be found in Eq.
(137)), $M_{3}^{XX}(N)=-200\pm 56\;(-111\pm 31.08)$ Barea:2015kwa ;
Deppisch:2020ztt  for ${}^{76}{\rm Ge}\;(^{136}{\rm Xe})$ is the NME
corresponding to short range contributions which is defined as

$\displaystyle
M_{3}^{XX}(N)\equiv\langle{\mathcal{O}}_{F}^{+}|[4(\bar{u}\gamma_{\mu}P_{X}d)(\bar{u}\gamma^{\mu}P_{X}d)]|\mathcal{O}_{I}^{+}\rangle.$
(22)

##  III The LRSM

For the model LRSM, the gauge fields are $SU(3)_{C}\bigotimes
SU(2)_{L}\bigotimes SU(2)_{R}\bigotimes U(1)_{B-L}$, and the details about the
gauge fields and their breaking can be found in Refs. Dev:2016dja ;
Patra:2015bga ; Mitra:2016kov ; Maiezza:2016ybz . In this model, the tiny
neutrino masses are obtained by both of type-I and type-II seesaw mechanisms
due to introducing the right-handed neutral leptons and two triplet Higgs
(scalars) accordingly. In the model, the mass matrix for the neutral leptons
generally is written as

$\displaystyle\left(\begin{array}[]{c}M_{L},\;\;\;\;M_{D}^{\rm T}\\\
M_{D},\;\;\;\;M_{R}\end{array}\right)$ (25)

and the mass matrix Eq. (25)222The matrix Eq. (25) with $M_{L}=0$ indicates
the masses are acquired by type-I seesaw mechanism; it with $M_{D}=0$
indicates the masses are acquired by type-II seesaw mechanism; it in general
feature indicates the masses are acquired by type-I+II seesaw mechanism. can
be diagonalized in terms of a unitary matrix $U_{\nu}$, whereas the matrix
$U_{\nu}$ can be expressed similarly as that in the case of the B-LSSM Eq.
(11).

For the model LRSM, if the left-right symmetry is not broken manifestly but
spontaneously, i.e. $g_{L}=g_{R}\equiv g_{2}$ and as one consequence, the mass
terms of $W$ bosons can be written as

$\displaystyle\mathcal{L}_{M_{W}}=\frac{g_{2}^{2}}{4}\left(\begin{array}[]{c}W_{L}^{+},\;\;W_{R}^{+}\end{array}\right)\left(\begin{array}[]{c}v_{1}^{2}+v_{2}^{2}+2v_{L}^{2},\quad
2v_{1}v_{2}\\\ 2v_{1}v_{2},\quad
v_{1}^{2}+v_{2}^{2}+2v_{R}^{2}\end{array}\right)\left(\begin{array}[]{c}W_{L}^{-}\\\
W_{R}^{-}\end{array}\right),$ (31)

where $v_{1},v_{2},v_{L},v_{R}$ ($v_{L}\ll v_{R}$) are the VEVs of new scalars
(Higgs) in the LRSM. Then the physical masses of the $W$ bosons can be
obtained Deppisch:2017vne

$\displaystyle
M_{W_{1}}\simeq\frac{g_{2}}{2}(v_{1}^{2}+v_{2}^{2})^{1/2},\;\;\;\;\;M_{W_{2}}\simeq\frac{g_{2}}{\sqrt{2}}v_{R}.$
(32)

The mass eigenstates $W_{1,2}^{\pm}$ are related to the interaction
eigenstates $W_{L,R}^{\pm}$ by $\zeta$

$\displaystyle\left(\begin{array}[]{c}W_{1}^{\pm}\\\
W_{2}^{\pm}\end{array}\right)=\left(\begin{array}[]{c}\cos\zeta,\;\;\sin\zeta\\\
-\sin\zeta,\;\;\cos\zeta\end{array}\right)\left(\begin{array}[]{c}W_{L}^{\pm}\\\
W_{R}^{\pm}\end{array}\right),$ (39)

where $\tan 2\zeta=\frac{2v_{1}v_{2}}{v_{R}^{2}-v_{L}^{2}}$.

The interactions, being applied later on, in the model are

$\displaystyle\mathcal{L}_{I}^{LRSM}=\frac{{\rm
i}g_{2}}{\sqrt{2}}\sum_{j=1}^{3}\Big{[}\bar{e}_{i}(\cos\zeta
U_{ij}\gamma^{\mu}P_{L}+\sin\zeta
T^{*}_{ij}\gamma^{\mu}P_{R})\nu_{j}W_{1,\mu}^{-}$
$\displaystyle\qquad\qquad\;\;\;+\bar{e}_{i}(\cos\zeta
T^{*}_{ij}\gamma^{\mu}P_{R}-\sin\zeta
U_{ij}\gamma^{\mu}P_{L})\nu_{j}W_{2,\mu}^{-}$
$\displaystyle\qquad\qquad\;\;\;+\bar{e}_{i}(\cos\zeta
S_{ij}\gamma^{\mu}P_{L}+\sin\zeta
V^{*}_{ij}\gamma^{\mu}P_{R})N_{j}W_{1,\mu}^{-}$
$\displaystyle\qquad\qquad\;\;\;+\bar{e}_{i}(\cos\zeta
V^{*}_{ij}\gamma^{\mu}P_{R}-\sin\zeta
S_{ij}\gamma^{\mu}P_{L})N_{j}W_{2,\mu}^{-}$
$\displaystyle\qquad\qquad\;\;\;+\bar{u}(\cos\zeta\gamma^{\mu}P_{L}+\sin\zeta\gamma^{\mu}P_{R})dW_{1,\mu}^{-}$
$\displaystyle\qquad\qquad\;\;\;+\bar{u}(\cos\zeta\gamma^{\mu}P_{R}-\sin\zeta\gamma^{\mu}P_{L})dW_{2,\mu}^{-}+h.c\Big{]},$
(40)

where the definitions for $U,\;S,\;T,\;V,\;\nu,\;N$ are same as the ones in
the B-LSSM.

In the LRSM, the dominant contributions to the $0\nu 2\beta$ decays are
represented by Feynman diagrams Fig. 2.

Figure 2: The dominant Feynman diagrams for the $0\nu 2\beta$ decays in the
LRSM. (a) The contributions from the heavy neutral lepton exchanges, (b) The
contributions from the light neutrino exchanges.

In comparison with those of the B-LSSM, the contributions from the Higgs
exchanges can be ignored so in Fig. 2 there is no Higgs exchange at all, but
besides the left-handed gauge boson components $W_{L}^{\pm}$, there are right-
handed components $W_{R}^{\pm}$ in $W_{1}^{\pm}$ and $W_{2}^{\pm}$ gauge
bosons. Therefore the situation in determining the effective Lagrangian for
the $0\nu 2\beta$ decays is different and comparatively complicated than that
for the B-LSSM. In the case for the contributions from the heavy neutral
lepton exchanges shown in Fig. 2 (a), considering the fact that the heavy
neutral leptons propagator
$\frac{M_{N_{i}}+k\\!\\!\\!/}{k^{2}-M_{N_{i}}^{2}}\approx\frac{-1}{M_{N_{i}}}$
($M_{N_{i}}\geq M_{W_{1}}$) and $\sin\zeta,S_{1i},T_{1i}\ll 1.0$ in the
decays, by using Eqs. (13, 16) the effective Lagrangian at the energy scale
$\mu\simeq M_{W_{1}}$ ($W_{1}$ is the lighter one boson between $W_{1,2}$) can
be written down as follows

$\displaystyle\frac{2m_{p}}{G_{F}^{2}\cos\theta_{C}^{2}}\mathcal{L}_{\rm
eff}^{\rm
LR}(N)=C_{3R}^{LL}(N)\mathcal{O}_{3R}^{LL}+C_{3L}^{LL}(N)\mathcal{O}_{3L}^{LL}+C_{3L}^{RL}(N)\mathcal{O}_{3L}^{RL}+C_{3L}^{RR}(N)\mathcal{O}_{3L}^{RR}\,;$
$\displaystyle C_{3R}^{LL}(N)=\frac{2m_{p}}{M_{N_{i}}}\cos^{4}\zeta
S^{2}_{1i},\;\;C_{3L}^{LL}(N)=\frac{2m_{p}}{M_{N_{i}}}\cos^{2}\zeta\sin^{2}\zeta
V^{*2}_{1i}.$ $\displaystyle
C_{3L}^{RL}(N)=\frac{2m_{p}}{M_{N_{i}}}\cos^{3}\zeta\sin\zeta
V^{*2}_{1i}(\frac{M_{W_{L}}}{M_{W_{R}}})^{2},\;\;C_{3L}^{RR}(N)=\frac{2m_{p}}{M_{N_{i}}}\cos^{4}\zeta
V^{*2}_{1i}(\frac{M_{W_{L}}}{M_{W_{R}}})^{4}.$ (41)

In literatures, due to small $\zeta$, the contributions corresponding to
$C_{3L}^{RL}(N)$, $C_{3L}^{LL}(N)$ in Eq. (41) are neglected. However, since
$\tan 2\zeta=\frac{2v_{1}v_{2}}{v_{R}^{2}-v_{L}^{2}}$, i.e. $\zeta\approx
xM_{W_{1}}^{2}/M_{W_{2}}^{2}$ and when $x\equiv v_{2}/v_{1}>0.02$
Bertolini:2014sua , the terms with $C_{3L}^{RL}(N)$ and $C_{3L}^{LL}(N)$ can
also make essential contributions compared with the terms with
$C_{3L}^{RR}(N)$, $C_{3R}^{LL}(N)$. Thus in this work, when evaluating the
$0\nu 2\beta$ decays we would like to keep the contributions from the terms of
$C_{3L}^{RL}(N)$, $C_{3L}^{LL}(N)$, and consider the QCD corrections to the
effective Lagrangian Eq. (41) in the similar way as that in the B-LSSM.

As the next step, when considering the contributions from the light neutrino
exchanges as Fig. 2 (b), owing to the fact that $W^{\pm}_{1,2}$ contain both
components $W^{\pm}_{L,R}$ in LRSM, according to Eqs.(13,16), the light
neutrino propagators with chiral project operators
$P_{L,R}\frac{m_{\nu_{i}}+k\\!\\!\\!/}{k^{2}-m_{\nu_{i}}^{2}}P_{L,R}$ or
$P_{L,R}\frac{m_{\nu_{i}}+k\\!\\!\\!/}{k^{2}-m_{\nu_{i}}^{2}}P_{R,L}$, as
central factors, finally contribute the factors as
$\frac{m_{\nu_{i}}}{k^{2}}P_{L,R}$ or $\frac{k\\!\\!\\!/}{k^{2}}P_{L,R}$
respectively to the results. Namely the factors
$\frac{k\\!\\!\\!/}{k^{2}}P_{L,R}$ are new and substantial, and in B-LSSM they
do not appear at all. Moreover when the contributions from the ‘higher order’
terms for $\sin\zeta,\;S_{1i},\;T_{1i}$, such as those small terms
proportional to $\sin^{2}\zeta,\;S_{1i}^{2},\;T_{1i}^{2}$ etc, are ignored,
then the effective Lagrangian at the energy scale $\mu\simeq 1.0\;{\rm GeV}$
may be read out from Fig. 2(b) as

$\displaystyle\frac{2m_{p}}{G_{F}^{2}\cos\theta_{C}^{2}}\mathcal{L}_{\rm
eff}^{\rm LR}(\nu)=\cos^{4}\zeta\frac{U_{1i}^{2}m_{\nu
i}}{m_{e}}\frac{2m_{p}m_{e}}{-k^{2}}\mathcal{O}_{3R}^{LL}+$
$\displaystyle\qquad\qquad\qquad\qquad\quad\cos^{3}\zeta\sin\zeta
U_{1i}T_{1i}^{*}[4(\bar{u}\gamma_{\mu}P_{L}d)(\bar{u}\gamma_{\nu}P_{L}d)\bar{e}\gamma^{\mu}\frac{2m_{p}k\\!\\!\\!/}{k^{2}}\gamma^{\nu}e^{c}]+$
$\displaystyle\qquad\qquad\qquad\qquad\quad\cos^{4}\zeta
U_{1i}T_{1i}^{*}\frac{M_{W_{L}}^{2}}{M_{W_{R}}^{2}}[4(\bar{u}\gamma_{\mu}P_{R}d)(\bar{u}\gamma_{\nu}P_{L}d)\bar{e}\gamma^{\mu}\frac{2m_{p}k\\!\\!\\!/}{k^{2}}\gamma^{\nu}P_{R}e^{c}].$
(42)

The QCD corrections in the energy scale region $\mu\simeq 1.0\;$GeV to
$\mu\simeq 0.10\;$GeV, being of non-perturbative QCD, are taken into account
by inputting in the experimental measurements for the relevant current matrix
elements of nucleons, which emerge at the effective Lagrangian at $\mu\simeq
0.10\;$GeV.

In Refs. Mohapatra:1979ia ; Mohapatra:1981pm ; Picciotto:1982qe ;
Hirsch:1996qw ; Tello:2010am ; Chakrabortty:2012mh ; Barry:2013xxa ;
Dev:2013vxa ; Dev:2013oxa ; Huang:2013kma ; Dev:2014xea ; Ahmed:2019vum , the
second term and the third term of Eq. (42) are defined as $\eta$, $\lambda$
respectively. Extracting the factors

$\displaystyle C_{\eta}=\cos^{3}\zeta\sin\zeta
U_{1i}T_{1i}^{*},\;\;C_{\lambda}=\cos^{4}\zeta
U_{1i}T_{1i}^{*}\frac{M_{W_{1}}^{2}}{M_{W_{2}}^{2}},$ (43)

then the operators
$[4(\bar{u}\gamma_{\mu}P_{L}d)(\bar{u}\gamma_{\nu}P_{L}d)\bar{e}\gamma^{\mu}\frac{2m_{p}k\\!\\!\\!/}{k^{2}}\gamma^{\nu}e^{c}]$,
$[4(\bar{u}\gamma_{\mu}P_{R}d)(\bar{u}\gamma_{\nu}P_{L}d)\bar{e}\gamma^{\mu}\frac{2m_{p}k\\!\\!\\!/}{k^{2}}\gamma^{\nu}P_{R}e^{c}]$
are attributed to the calculations of NME and PSF.

Whereas when calculating the NMEs and PSF, the interference effects among the
contributions, especially to consider the contributions from the factors
$\frac{k\\!\\!\\!/}{k^{2}}P_{L,R}$ for the light neutrino exchanges, are
complicated and hard (in literatures to treat them even the Lorentz covariance
is lost Tomoda:1990rs ). In this work to calculate NMEs and PSF, we try to
make an additional approximation on the contributions relevant the factors
$\frac{k\\!\\!\\!/}{k^{2}}P_{L,R}$ for the light neutrino exchanges, which we
call as ‘frozen approximation’. Under the approximation, the momenta of the
two involved quarks inside the initial nucleus and two involved quarks inside
the final nucleus (Fig. 2) are assumed to be equal approximately:

$p_{1}\simeq p_{2}\equiv\bar{p},\;\;k_{1}\simeq k_{2}\equiv\overline{k}.$ (44)

With the ‘frozen approximation’ and ‘the on-shell approximation’ onto the
momenta for the out legs as well, all contributions corresponding to Fig. 2
(b) can be well-deduced and the final results can be collected as

$\displaystyle\frac{2m_{p}}{G_{F}^{2}\cos\theta_{C}^{2}}\mathcal{L}_{\rm
eff}^{\rm
LR}(\nu)=\frac{2m_{p}m_{e}}{-k^{2}}[C_{3R}^{LL}(\nu)\mathcal{O}_{3R}^{LL}+C_{3L}^{LL}(\nu)\mathcal{O}_{3L}^{LL}+C_{3L}^{RL}(\nu)\mathcal{O}_{3L}^{RL}+C_{3R}^{RL}(\nu)\mathcal{O}_{3R}^{RL}+$
$\displaystyle\qquad\qquad\qquad\qquad\quad
C_{5}^{LL}(\nu)\mathcal{O}_{5}^{LL}+C_{5}^{RR}(\nu)\mathcal{O}_{5}^{RR}+C_{5}^{LR}(\nu)\mathcal{O}_{5}^{LR}+C_{5}^{RL}(\nu)\mathcal{O}_{5}^{RL}],$
(45)

where

$\displaystyle\mathcal{O}_{5}^{XY}=4(\bar{u}\gamma_{\mu}P_{X}d)(\bar{u}P_{Y}d)\bar{e}\gamma_{\mu}\gamma^{5}e^{c},$
$\displaystyle C_{3R}^{LL}(\nu)=\frac{1}{m_{e}}\cos^{3}\zeta
U_{1i}(m_{\nu_{i}}\cos\zeta U_{1i}-m_{e}\sin\zeta T^{*}_{1i}),$ $\displaystyle
C_{3L}^{LL}(\nu)=-\cos^{3}\zeta\sin\zeta U_{1i}T^{*}_{1i},$ $\displaystyle
C_{3L}^{RL}(\nu)=C_{3R}^{RL}(\nu)=-\frac{1}{2}\cos^{4}\zeta
U_{1i}T^{*}_{1i}(\frac{M_{W_{1}}}{M_{W_{2}}})^{2},$ $\displaystyle
C_{5}^{RR}(\nu)=-C_{5}^{LL}(\nu)=\frac{m_{u}-m_{d}}{m_{e}}\cos^{3}\zeta\sin\zeta
U_{1i}T^{*}_{1i}-\frac{m_{d}}{m_{e}}\cos^{4}\zeta
U_{1i}T^{*}_{1i}(\frac{M_{W_{1}}}{M_{W_{2}}})^{2},$ $\displaystyle
C_{5}^{LR}(\nu)=-C_{5}^{RL}(\nu)=\frac{m_{d}-m_{u}}{m_{e}}\cos^{3}\zeta\sin\zeta
U_{1i}T^{*}_{1i}-\frac{m_{u}}{m_{e}}\cos^{4}\zeta
U_{1i}T^{*}_{1i}(\frac{M_{W_{1}}}{M_{W_{2}}})^{2}.$ (46)

Then in LRSM the half-life of $0\nu 2\beta$ decays can be written as
Deppisch:2020ztt

$\displaystyle\frac{1}{T_{1/2}^{0\nu}}=G^{0\nu}|M^{0\nu}|^{2}\Big{|}\frac{m_{ee}^{\rm
LR}}{m_{e}}\Big{|}^{2},$ (47)

where

$\displaystyle m_{ee}^{\rm
LR}=m_{e}\Big{\\{}\Big{|}[C_{3L}^{RR}(N)+C_{3L}^{LL}(N)]U_{3}^{XX}\frac{M_{3}^{XX}(N)}{M^{0\nu}}+C_{3L}^{RL}(N)U_{(31)11}^{XY}\frac{M_{3}^{XY}(N)}{M^{0\nu}}$
$\displaystyle\qquad\quad+C_{3L}^{LL}(\nu)+C_{3L}^{RL}(\nu)\frac{M_{3}^{XY}(\nu)}{M^{0\nu}}\Big{|}^{2}+\Big{|}C_{3R}^{LL}(N)U_{3}^{XX}\frac{M_{3}^{XX}(N)}{M^{0\nu}}+C_{3R}^{LL}(\nu)\Big{|}^{2}$
$\displaystyle\qquad\quad+2\frac{G_{11-}^{(0)}}{G^{0\nu}}\Big{[}[C_{3L}^{RR}(N)+C_{3L}^{LL}(N)]U_{3}^{XX}\frac{M_{3}^{XX}(N)}{M^{0\nu}}+C_{3L}^{RL}(N)U_{(31)11}^{XY}\frac{M_{3}^{XY}(N)}{M^{0\nu}}$
$\displaystyle\qquad\quad+C_{3L}^{LL}(\nu)+C_{3L}^{RL}(\nu)\frac{M_{3}^{XY}(\nu)}{M^{0\nu}}\Big{]}\Big{[}C_{3R}^{LL}(N)U_{3}^{XX}\frac{M_{3}^{XX}(N)}{M^{0\nu}}+C_{3R}^{LL}(\nu)\Big{]}\Big{\\}}^{1/2},$
(48)

where $U_{(31)}^{XY}$ is the $2\times 2$ RGE evolution matrix from $\mu\simeq
M_{W_{1}}$ to $\mu\simeq 1.0\;$GeV (the numerical result of $U_{(31)}^{XY}$
can be found in Eq. (137)), $C_{3Z}^{XY}(N,\nu)$ is the coefficient defined in
Eq. (46), $G^{(0)}_{11-}=-0.28\times 10^{-15}\;(-1.197\times 10^{-15})\;{\rm
years}^{-1}$ Deppisch:2020ztt  for ${}^{76}{\rm Ge}\;(^{136}{\rm Xe})$ is the
PSF, and $M_{3}^{XY}(\nu)=4.24\pm 0.68\;(2.17\pm 0.35)$,
$M_{3}^{XY}(N)=99.8\pm 27.94\;(51.2\pm 14.34)$ Barea:2015kwa ;
Deppisch:2020ztt  for ${}^{76}{\rm Ge}\;(^{136}{\rm Xe})$ are the NMEs
corresponding to the exchange of light neutrinos, heavy neutral leptons
respectively which are defined as

$\displaystyle
M_{3}^{XY}(\nu)\equiv\langle{\mathcal{O}}_{F}^{+}|\frac{2m_{p}m_{e}}{-k^{2}}[4(\bar{u}\gamma_{\mu}P_{X}d)(\bar{u}\gamma^{\mu}P_{Y}d)]|\mathcal{O}_{I}^{+}\rangle,$
$\displaystyle
M_{3}^{XY}(N)\equiv\langle{\mathcal{O}}_{F}^{+}|[4(\bar{u}\gamma_{\mu}P_{X}d)(\bar{u}\gamma^{\mu}P_{Y}d)]|\mathcal{O}_{I}^{+}\rangle.$
(49)

Similar to the case of B-LSSM, the factor $\frac{2m_{p}m_{e}}{-k^{2}}$ in Eq.
(45) is absorbed into the ‘neutrino potential’ which is used to compute the
long range NME. In addition, we should note that in Eq. (48) the terms of
$C_{5}^{RR}(\nu),C_{5}^{LL}(\nu),C_{5}^{RL}(\nu),C_{5}^{LR}(\nu)$ in Eq. (46)
do not make contributions to the $0\nu 2\beta$ decays for the reason below.
They make contributions to $0\nu 2\beta$ decays in the form

$\displaystyle
M_{5}^{XX}(\nu)[C_{5}^{RR}(\nu)+C_{5}^{LL}(\nu)]+M_{5}^{XY}(\nu)[C_{5}^{RL}(\nu)+C_{5}^{LR}(\nu)],$
(50)

where

$\displaystyle
M_{5}^{XX}(\nu)\equiv\langle{\mathcal{O}}_{F}^{+}|\frac{2m_{p}m_{e}}{-k^{2}}[4(\bar{u}\gamma_{\mu}P_{X}d)(\bar{u}P_{X}d)]|\mathcal{O}_{I}^{+}\rangle,$
$\displaystyle
M_{5}^{XY}(\nu)\equiv\langle{\mathcal{O}}_{F}^{+}|\frac{2m_{p}m_{e}}{-k^{2}}[4(\bar{u}\gamma_{\mu}P_{X}d)(\bar{u}P_{Y}d)]|\mathcal{O}_{I}^{+}\rangle.$
(51)

Since $C_{5}^{RR}(\nu)=-C_{5}^{LL}(\nu)$, $C_{5}^{RL}(\nu)=-C_{5}^{LR}(\nu)$
as shown in Eq. (46), then Eq. (50) shows the contributions to the decays from
the operators corresponding to the coefficients $C_{5}^{RR}(\nu)$,
$C_{5}^{LL}(\nu)$, $C_{5}^{RL}(\nu)$, $C_{5}^{LR}(\nu)$ are cancelled
completely.

##  IV Numerical results

With the formulas obtained above, now in this section we do the numerical
calculations and present the results on the $0\nu 2\beta$ decays for the the
nuclei ${}^{76}{\rm Ge}$ and ${}^{136}{\rm Xe}$ accordingly. In our numerical
calculations, the parameters are taken as the weak boson mass:
$M_{W_{L}}=80.385$ GeV for B-LSSM and $M_{W_{1}}=80.385$ GeV for LRSM,
$m_{b}=4.65$ GeV for b-quark mass, $m_{c}=1.275$ GeV for c-quark mass,
$\alpha_{em}(m_{Z})=1/128.9$ for the coupling of the electromagnetic
interaction, $\alpha_{s}(m_{Z})=0.118$ for the coupling of the strong
interaction. The constraints from available experimental data, such as the
most stringent upper limit on the sum of neutrino masses by PLANK
Esteban:2018azc  $\sum_{i}m_{\nu_{i}}<0.12\;{\rm eV}$; the neutrino mass-
squared differences obtained via analyzing the solar and atmospheric neutrino
oscillation data at $3\sigma$ deviations Esteban:2018azc

$\displaystyle\Delta m_{12}^{2}\equiv m_{\nu_{2}}^{2}-m_{\nu_{1}}^{2}=(7.4\pm
0.61)\times 10^{-5}\;{\rm eV}^{2},$
$\displaystyle\left\\{\begin{array}[]{l}\Delta m_{13}^{2}\equiv
m_{\nu_{3}}^{2}-m_{\nu_{1}}^{2}\approx(2.526\pm 0.1)\times 10^{-3}\;{\rm
eV}^{2}\qquad\qquad({\rm NH})\\\ \Delta m_{32}^{2}\equiv
m_{\nu_{2}}^{2}-m_{\nu_{3}}^{2}\approx(2.508\pm 0.1)\times 10^{-3}\;{\rm
eV}^{2}\qquad\qquad({\rm IH}),\\\ \end{array}\right.$ (54)

etc, are well-considered. Owing to the fact that the hierarchy of neutrino
masses has not been fixed yet, we take the two possibilities below to carry
out the analyses i.e. the normal hierarchy (NH)
$m_{\nu_{1}}<m_{\nu_{2}}<m_{\nu_{3}}$ and the inverse hierarchy (IH)
$m_{\nu_{3}}<m_{\nu_{1}}<m_{\nu_{2}}$. Moreover, we use the matrix $U$, the
upper-left sub-matrix of the whole matrix $U_{\nu}$ in Eq. (11), being the
Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix Tanabashi:2018oca  to
describe the mixing of the light neutrinos. Then we have

$\displaystyle(U_{1i})^{2}m_{\nu
i}=c_{12}^{2}c_{13}^{2}e^{i\alpha}m_{\nu_{1}}+c_{13}^{2}s_{12}^{2}e^{i\beta}m_{\nu_{2}}+s_{13}^{2}m_{\nu_{3}},$
(55)

where $\alpha,\beta$ are Majorana CP violating phases,
$c_{12,13}\equiv\cos\theta_{12,13},s_{12,13}\equiv\sin\theta_{12,13}$ are the
neutrino oscillation parameters and at $3\sigma$ error level they can be
written as

$\displaystyle s_{12}^{2}=0.3125\pm 0.0375,\;s_{13}^{2}=0.022405\pm 0.001965.$
(56)

The direct searching for the right-handed gauge boson sets the lower bound for
the mass of $W_{2}$ boson as $M_{W_{2}}\gtrsim 4.8\;{\rm TeV}$
Sirunyan:2017yrk ; Sirunyan:2018pom ; Aaboud:2018spl ; Aaboud:2019wfg , and
the $W_{R}$, $W_{L}$ mixing angle as $\zeta\lesssim 7.7\times 10^{-4}$
Dev:2014xea . Finally we have $x\equiv v_{2}/v_{1}>0.02$ Bertolini:2014sua
required by the CP-violation data for the $K$ and $B$ mesons, where $v_{1}$,
$v_{2}$ are the VEVs of the two Higgs particles in the LRSM.

Eq. (21) and Eq. (48) shows $m_{ee}^{LR}$ depends the NME and phase space
factor of the chosen nuclei. However these factors appears in the expression
of $m_{ee}^{LR}$ in the form of some ratios, and these ratios for ${}^{76}{\rm
Ge}$ are equal to those for ${}^{136}{\rm Xe}$ roughly, hence we can use the
same expression of $m_{ee}^{LR}$ for ${}^{76}{\rm Ge}$, ${}^{136}{\rm Xe}$ in
the analyses later on.

###  IV.1 The numerical results for the B-LSSM

Figure 3: $m_{ee}^{\rm BL}$ versus $m_{\nu-{\rm lightest}}$ under scanning the
neutrino mass-squared differences Eq. (54) at $3\sigma$ error level and the
parameter regions Eq. (57). The green (red) points denote the NH (IH) results,
the blue (purple) solid line denotes the experimental constraints from the
$0\nu 2\beta$ decay half-life of ${}^{76}{\rm Ge}$ (${}^{136}{\rm Xe}$), the
blue (purple) dashed line denotes the experimental ability of ${}^{76}{\rm
Ge}$ (${}^{136}{\rm Xe}$) for the next generation of experiments, the orange
solid (dashed) line denotes the constraints from PLANK 2018 for the case of NH
(IH).

According to the above analysis, in B-LSSM the contributions from the heavy
neutral leptons shown in Eq. (17) are highly suppressed ($S_{1i}\approx
10^{-7}$) for the TeV-scale heavy neutral leptons, hence the dominant
contributions to the $0\nu 2\beta$ decays come from the light neutrinos. Then
with the neutrino mixing parameters

$\displaystyle s_{12}^{2}=(0.275\sim 0.35),\;\;s_{13}^{2}=(0.02044\sim
0.02437),\;\;\alpha=(0\sim 2\pi),\;\;\beta=(0\sim 2\pi).$ (57)

and the neutrino mass-squared differences as those in Eq. (54) at $3\sigma$
error level, $m_{ee}^{\rm BL}$ versus $m_{\nu-{\rm lightest}}$ (the lightest
neutrino mass) for the B-LSSM is plotted in Fig. 3, where the green (red)
points denote the NH (IH) results, the blue (purple) solid line denotes the
constraints from the lower $0\nu 2\beta$ decay half-life bound of ${}^{76}{\rm
Ge}$ (${}^{136}{\rm Xe}$), the blue (purple) dashed line denotes the
experimental ability of ${}^{76}{\rm Ge}$ (${}^{136}{\rm Xe}$) for the next
generation of experiments, the orange solid (dashed) line denotes the
constraints from PLANK 2018 for NH (IH) neutrino masses (the meaning for the
red points, green points, blue lines, purple lines, orange lines is also
adopted accordingly in the figures later on). In the figure Fig. 3,
$m_{ee}^{\rm BL}$ is well-below the experimental upper bounds in the cases of
NH and IH, and there is not tighter restriction on the $m_{\nu-{\rm
lightest}}$ than that offered by PLANK. Additionally, the blue dashed line
shows that there is certain opportunity to observe the $0\nu 2\beta$ decays in
the next generation of experiments, whereas if any of the decays is not
observed in the next generation of experiments, then the IH neutrino masses
will be excluded completely by the blue dashed line in Fig. 3.

###  IV.2 The numerical results for the LRSM

In the LRSM owing to the bosons $W_{L}$, $W_{R}$ and their mixing, the
situation is much more complicated than that in the B-LSSM, and both of the
light neutrinos and heavy neutral leptons make substantial contributions to
the $0\nu 2\beta$ decays. For TeV-scale heavy neutral leptons in the LRSM, the
consequences of the type II seesaw dominance are similar to those of the type
I seesaw dominance, while the consequences are very different from the ones of
type I$+$II seesaw dominance, because the light-heavy neutral lepton mixing is
not much suppressed in this case. Hence we do the numerical computations for
the type I and type I$+$II seesaw dominance cases in our analyses. For
simplicity and not losing general features, we assume that there is only one
of heavy neutral leptons to make substantial contributions in the model. It
indicates that there are no off-diagonal elements in the matrix $M_{R}$ in Eq.
(25), i.e. we have $M_{R}\simeq\hat{M}_{N}={\rm
diag}(M_{N_{1}},M_{N_{2}},M_{N_{3}})$. Then under the dominance of either the
type I seesaw or the type I+II seesaw mass mechanisms, we compute $m_{ee}^{\rm
LR}$ numerically in tern.

####  IV.2.1 The results under type I seesaw dominance

Under the type I seesaw dominance and due to the tiny neutrino masses, the
sub-matrix $T$ in Eq. (11) has $T_{1j}\ll 1\;(j=1,2,3)$. It indicates that the
contributions from the terms proportional to $T_{1i}^{*}$ as shown in Eq. (46)
are highly suppressed. Then scanning the parameter space in Eqs. (54, 57) and
in the parameter space

$\displaystyle x=(0.02\sim 0.5),\;\;M_{W_{2}}=(4.8\sim 10.0)\;{\rm
TeV},\;\;M_{N_{1}}=(0.10\sim 3.0)\;{\rm TeV},$ (58)

we plot $m_{ee}^{\rm LR}$ versus $m_{\nu-{\rm lightest}}$ in Fig. 4.

Figure 4: $m_{ee}^{\rm LR}$ versus $m_{\nu-{\rm lightest}}$ with scanning the
parameter space in Eqs. (54, 57, 58).

Comparing Fig. 4 with Fig. 3, the red points show that the range of
$m_{ee}^{\rm LR}$ is similar to the range of $m_{ee}^{\rm BL}$ in the case of
IH neutrino masses, but in the case of NH neutrino masses there are points
with $m_{ee}^{\rm LR}>m_{ee}^{\rm BL}$. For IH neutrino masses, the
contributions to $m_{ee}^{\rm LR}$ are dominated by the terms proportional to
light neutrino masses, hence $m_{ee}^{\rm LR}$ depends on
$x,M_{W_{2}},M_{N_{1}}$ negligibly, that leads to the the range of reds points
is similar to the results of $m_{ee}^{\rm BL}$. Additionally, $m_{ee}^{\rm
LR}>m_{ee}^{\rm BL}$ shown as green points indicates the contributions to
$m_{ee}^{\rm LR}$ can be dominated by heavy neutral leptons for appropriate
values of $x,M_{W_{2}},M_{N_{1}}$.

In order to see the effects of $x,M_{W_{2}},M_{N_{1}}$ clearly, we take
$m_{\nu_{1}}=0.001\;{\rm eV}$ for the NH neutrino masses,
$s_{12},s_{13},\Delta m_{12}^{2},\Delta m_{13}^{2}$ at the corresponding
center values and the CP violation phases $\alpha=\beta=0$. Then taking
$M_{W_{2}}=5.0\;{\rm TeV}$, we plot $m_{ee}^{\rm LR}$ versus $M_{N_{1}}$ in
Fig. 5 (a), where the solid, dashed, dotted lines denote the results for
$x=0.10,\,0.25,\,0.40$ respectively. Similarly, $m_{ee}^{\rm LR}$ versus
$M_{W_{2}}$ for $M_{N_{1}}=0.20\;{\rm TeV}$ is plotted in Fig. 5 (b).

Figure 5: With $m_{\nu_{1}}=0.001\;{\rm eV}$ for NH neutrino masses, (a):
$m_{ee}^{\rm LR}$ versus $M_{N_{1}}$ for $M_{W_{2}}=5.0\;{\rm TeV}$, (b):
$m_{ee}^{\rm LR}$ versus $M_{W_{2}}$ for $M_{N_{1}}=0.2\;{\rm TeV}$. The
solid, dashed, dotted lines denote the obtained $m_{ee}^{\rm LR}$ for
$x=0.10,\,0.25,\,0.40$ respectively.

From Fig. 5, one may see the fact that the obtained $m_{ee}^{\rm LR}$
decreases with the increasing of $M_{N_{1}}$, $M_{W_{2}}$, then $m_{ee}^{\rm
LR}$ approaches to a constant when $M_{N_{1}}$ or $M_{W_{2}}$ becomes large.
According to the definition of $W_{L}$-$W_{R}$ mixing parameter $\zeta$, the
coefficient $C_{3R}^{RL}(N)$ in Eq. (41) increases with increasing
$x\equiv\frac{v_{2}}{v_{1}}$, which leads to that $m_{ee}^{\rm LR}$ increases
with $x$ increasing as shown in the figures. Additionally one may see that
$m_{ee}^{\rm LR}$ depends on the values of $x,\;M_{W_{2}},\;M_{N_{1}}$ mildly
when $M_{N_{1}}$ or $M_{W_{2}}$ is large. It is because that the contributions
from heavy neutral lepton(s) are highly suppressed when its mass or right
handed boson mass becomes large, thus the contributions from the light
neutrinos become dominant and proportional merely to the light neutrino masses
in this case.

####  IV.2.2 The results under the type I+II seesaw dominance

As pointed out above, the mixing parameters of the light neutrinos and the
heavy neutral leptons are not tiny under the type I$+$II seesaw dominance,
hence the Dirac mass matrix $M_{D}$ can also affect the numerical results via
the mixing of the light neutrinos and the heavy neutral leptons. For
simplicity and not losing general feature, we assume that the mass matrix
$M_{D}$ in Eq. (25) is diagonal as $M_{D}={\rm
diag}(M_{D11},M_{D22},M_{D33})$. Taking the assumption that only one
generation of heavy neutral leptons makes substantial contributions, i.e.
$M_{D22}$ and $M_{D33}$ do not affect the results. Scanning the ranges
determined by Eqs. (54, 57) and the parameter space

$\displaystyle x=(0.02\sim 0.5),\,M_{W_{2}}=(4.8\sim 10)\;{\rm
TeV},M_{N_{1}}=(0.10\sim 3.0)\;{\rm TeV},$ $\displaystyle M_{D11}=(0.10\sim
100)\;{\rm MeV},$ (59)

$m_{ee}^{\rm LR}$ versus $m_{\nu-{\rm lightest}}$ as the results is plotted in
Fig. 6.

Figure 6: $m_{ee}^{\rm LR}$ versus $m_{\nu-{\rm lightest}}$ with scanning the
parameter space in Eqs. (54, 57, 59).

Comparing Fig. 6 with Fig. 4 (also Fig. 3 ), a lager $m_{ee}^{\rm LR}$ can be
reached in the case of type I+II seesaw dominance. Because the contributions
are dominated by the terms proportional to $S_{1i}$ or $T_{1i}^{*}$ when
$M_{D11}$ is large (shown in Eqs. (41, 46)), and lager $m_{ee}^{\rm LR}$ is
reached in this case. As indicated in Fig. 6 there are some points which are
smaller than the points in Fig. 4 both for NH and IH neutrino masses, that is
due to the existence of the cancellation effect, i.e. the contributions from
the terms proportional to $M_{D11}$ can be cancelled some amount of the
contributions from the other terms. The parameter $M_{D11}$ represents the
strength of light-heavy neutral lepton mixing, that we can introduce a
familiar parameter $S_{e}^{2}\equiv\sum_{i=1}^{3}|S_{1i}|^{2}$ for the
following analysis to describe the light-heavy neutral lepton mixing.

In order to show the effects of light-heavy mixing parameter $S_{e}^{2}$ and
the cancellation effect well, we take $M_{N_{1}}=0.2\;{\rm TeV}$,
$m_{\nu_{1}}=0.001\;{\rm eV}$ for the NH neutrino masses,
$s_{12}^{2},\,s_{13}^{2},\,\Delta m_{12}^{2},\,\Delta m_{13}^{2}$ with the
center values and the CP violation phases $\alpha=\beta=0$ accordingly. The
effects with the parameters being fixed above are presented a similar behavior
as the case of type I seesaw dominance, hence we do not repeat the study of
the effects on these parameters. With $M_{W_{2}}=5.0\;{\rm TeV}$, in Fig. 7
(a) we plot $m_{ee}^{\rm LR}$ versus $S_{e}^{2}$, and the solid, dashed,
dotted lines denote the results for $x=0.02,\,0.06,\,0.1$ respectively.
Similarly, with $x=0.02$, $m_{ee}^{\rm LR}$ versus $S_{e}^{2}$ is plotted in
Fig. 7 (b), where the solid, dashed, dotted lines denote the results for
$M_{W_{2}}=5.0,\,6.0,\,7.0\;{\rm TeV}$ respectively.

Figure 7: $m_{ee}^{\rm LR}$ versus $S_{e}^{2}$ with $M_{N_{1}}=0.2\;{\rm
TeV}$, $m_{\nu_{1}}=0.001\;{\rm eV}$ for the NH neutrino masses, $s_{12}$,
$s_{13}$, $\Delta m_{12}^{2}$, $\Delta m_{13}^{2}$ at the corresponding center
values and the CP violation phases $\alpha=\beta=0$. (a): the results for
$M_{W_{2}}=5.0\;{\rm TeV}$, and the solid, dashed, dotted lines denote the
results for $x=0.02,\,0.06,\,0.10$ respectively. (b): the results for
$x=0.02$, and the solid, dashed, dotted lines denote the results for
$M_{W_{2}}=5.0,\,6.0,\,7.0\;{\rm TeV}$ respectively.

From Fig. 7(a, b) one may see the fact that with the mixing parameter
$S_{e}^{2}$ increasing, $m_{ee}^{\rm LR}$ decreases to a minimum value and
then increases. In the case of the type I+II seesaw dominance, i.e. $S_{1i}$
and $T_{1i}$ are not too small, the contributions from the light neutrinos are
dominated over those terms which do not depend on the neutrino mass
$m_{\nu_{i}}$ in Eq. (46). When $S_{e}^{2}$ is small, $C_{3L}^{RR}(N)$ in Eq.
(41) plays a dominant role, and as $S_{e}^{2}$ increasing, the contributions
from $C_{3L/R}^{LL}(\nu)$ in Eq. (46) become larger. The minimum values for
$m_{ee}^{\rm LR}$ as shown in Fig. 7 is due to the opposite signs for
$M_{3}^{XX}(N)$ and $M_{3}^{XY}(\nu)$, i.e. it is caused by cancellation of
the contributions from $C_{3L}^{RR}(N)$ and $C_{3L/R}^{RL}(\nu)$. From the
figures the fact can be seen clearly that when $S_{e}^{2}$ is increasing, the
contributions from $C_{3L/R}^{RL}(\nu)$ become dominant which leads to the
increasing of $m_{ee}^{\rm LR}$ as long as $S_{e}^{2}$ becomes large enough.
The results indeed show that with $S_{e}^{2}$ varying, the cancellation takes
place in a proper manner. When the cancellation takes place, the value of
$S_{e}^{2}$ depends on $M_{W_{2}}$ and $x$ explicitly.

To compare the results obtained under the fresh approximation in this work
with the ones obtained by traditional way, we take the results from Ref.
Dev:2014xea  and present $m_{ee}^{\rm LR}$ versus $M_{W_{2}}$ in Fig. 8 (a),
(b) for $s_{e}^{2}=0.25\times 10^{-8},\;10^{-8}$ respectively, where the
solid, dashed lines denote the results obtained under the fresh approximation
in this work and those by the traditions way from in Ref. Dev:2014xea
respectively.

Figure 8: (a): $m_{ee}^{\rm LR}$ versus $M_{W_{2}}$ for $S_{e}^{2}=0.25\times
10^{-8}$. (b): $m_{ee}^{\rm LR}$ versus $M_{W_{2}}$ for $S_{e}^{2}=10^{-8}$.
Both are with $M_{N_{1}}=0.2\;{\rm TeV}$, $m_{\nu_{1}}=0.01\;{\rm eV}$ for the
NH neutrino masses. The solid, dashed lines denote the results obtained under
the new approximation in this work and the traditions method shown in Ref.
Dev:2014xea  respectively.

In Fig. 8, we take $M_{N_{1}}=0.2\;{\rm TeV}$, $m_{\nu_{1}}=0.010\;{\rm eV}$,
$x=0.020$, $s_{12},\,s_{13},\,\Delta m_{12}^{2},\,\Delta m_{13}^{2}$ at the
corresponding center values and the CP violation phases $\alpha=\beta=0$. Note
that in the figures the results obtained in this work coincide well with the
traditional ones when $S_{e}^{2}$ is small or $M_{W_{2}}$ is large in the
chosen parameter space. In addition, the results obtained in this work, being
approximate ones but having the interference effects considered better, are
smaller than the ones obtained in traditional way when $W_{2}$ is not heavy
enough, and the reduction factor depends on the parameters in LRSM completely
and also comes from the uncertainties of NMEs partly. Considering the
difficulties and the uncertainties etc in computing NMEs and the interference
effects among the various contributions by the way in literatures and the
approximation approach in this work, from Fig.8, it seems that the results
with very heavy $W_{2}$ ($M_{W_{2}}\geq 12$TeV) or small $S_{e}^{2}$
$(\lesssim 10^{-9})$ may be convinced more as the results approach to
coinciding with each other.

##  V Summary and Discussions

In the paper we take B-L supersymmetric standard model (B-LSSM) and TeV scale
left-right symmetric model (LRSM) as two representations of two kinds of new
physics models to study the nuclear neutrinoless double beta decays ($0\nu
2\beta$). As stated in Introduction, the calculations, on which are focused
the lights, are those about the factor on the ’core’ factor of the decays:
evaluating the process $d+d\to u+u+e+e$ by considering the effective
Lagrangian containing the operators with the Wilson coefficients or
considering the relevant amplitude etc. Whereas here the estimations of the
other necessary ‘factors’ for the decays, i.e. to evaluate relevant nuclear
matrix element (NME) and phase space factor (PSF) etc, that does not relate to
the specific models, are treated by following literatures.

In the B-LSSM, all of the calculations can be well deduced and the results are
dominated by the neutrino mass terms. However, in the LRSM owing to the
existence of right-handed gauge boson $W_{R}$, the calculations are
complicated, and the interference effects are hard and not well to be
estimated. In this work, a new approximation, i.e. the momenta of the two
involved quarks inside the initial nucleus and inside the final nucleus are
assumed to be equal approximately, is made, then all contributions in the LRSM
can be well reduced and summed up all the contributions. Then the calculations
in LRSM are simplified and the interference effects can be calculated
comparatively easily. To see the consequences of the approximation, we compare
the results obtained in this work with the ones obtained by the traditional
method numerically, and the results coincide with each other well when light-
heavy neutral lepton mixing parameter $S_{e}^{2}$ is small or $M_{W_{2}}$ is
large. For the effective dimension-9 contact interactions in these two models,
the contributions from heavy neutral lepton exchange, the QCD corrections from
the energy scale $\mu\simeq M_{W_{L}}$ (or $M_{W_{1}}$) to the energy scale
$\mu\simeq 1.0\;$GeV to all dimension-9 operators in the effective Lagrangian
which is responsible for the $0\nu 2\beta$ decays are calculated by the RGE
method, and all the QCD corrections including the contributions from light
neutrnos in the energy scale region $\mu\simeq 1.0\;$GeV to $\mu\simeq
0.10\;$GeV, being of non-perturbative QCD, are taken into account
alternatively by inputting in the experimental measurements for the relevant
current matrix elements of nucleons, which emerge at the effective Lagrangian
at $\mu\simeq 0.10\;$GeV.

With necessary input parameters allowed by experimental data, the theoretical
predictions on $0\nu 2\beta$ decay half life ${}^{76}{\rm Ge}$ and
${}^{136}{\rm Xe}$ is obtained in these two models. In the B-LSSM, the
contributions from heavy neutral leptons are highly suppressed by the tiny
light-heavy neutral lepton mixing parameters for TeV-scale heavy neutral
leptons. Hence $m_{ee}^{\rm BL}$ depends on the light neutrino masses mainly,
and the numerical results show that the $0\nu 2\beta$ decays may be observed
with quite great opportunity in the near future. Whereas if the decays are not
observed in the next generation of $0\nu 2\beta$ experiments, the IH neutrino
masses are excluded completely by the lower bound on $T_{1/2}^{0\nu}(^{76}{\rm
Ge})$ (as shown in Fig. 3).

In the LRSM, the situation is much more complicated than that in the B-LSSM,
so we have learnt a lot of experiences in study the $0\nu 2\beta$ decays for
the model. As for the type I seesaw dominance, the contributions from the
terms proportional to $T_{1i}^{*}$ (see Eq. (46)) are highly suppressed by the
tiny neutrino masses. The numerical results of the contributions from light
neutrinos are similar to the ones for B-LSSM, but the heavy neutral leptons
can make comparatively large contributions through the right-handed current
when $M_{N_{i}}$ ($\lesssim 0.5\;{\rm TeV}$) and $M_{W_{2}}$ ($\lesssim
7\;{\rm TeV}$) both are not too heavy. For the type I+II seesaw dominance, the
terms which do not depend on the neutrino mass $m_{\nu_{i}}$ in Eq. (46) play
the dominant roles. In this case, the contributions from $C_{3L/R}^{RL}(L)$
can be cancelled some by the contributions from $C_{3L}^{RR}(H)$ when the
light-heavy mixing parameter $S_{e}^{2}$ is appropriate, because the signs of
the corresponding NMEs are opposite. Moreover the effects of the cancellation
are affected by $S_{e}^{2}$, $x$ and the right-handed $W$-boson mass
$M_{W_{2}}$ etc in a complicated way. In addition, Fig. 6 shows that the
points either in green (NH) or in red (IH) spread out a lot, and there are
many ‘exotic points’, where cannot be reached for the cases of the B-LSSM and
type I seesaw dominance LRSM. Thus the characteristic feature on the
distribution of $m_{ee}^{\rm LR}$ versus $m_{\nu-{\rm lightest}}$ may help to
realize whether the decays are caused by the LRSM in the type I+II seesaw
dominance or not, particularly when the $0\nu 2\beta$ decays are observed and
the points for $m_{ee}^{\rm LR}$ versus $m_{\nu-{\rm lightest}}$ just fall on
the exotic points.

Finally, according to the numerical results of the present comparative studies
on the $0\nu 2\beta$ decays for the two typical models LRSM and B-LSSM, it may
be concluded that the room of LRSM type models for the foreseeable future
decay experiments is greater than that of B-LSSM type models, and having the
right handed gauge bosons the feature of the LRSM type models is more
complicated than that of the B-LSSM type models.

Acknowledgments: This work was supported in part by the National Natural
Science Foundation of China (NNSFC) under Grants No. 11821505 No. 12047503,
No. 12075301 and No. 11705045. It was also supported in part by the Key
Research Program of Frontier Sciences, CAS, Grant No. QYZDY-SSW-SYS006. The
authors (J.-L. Yang and C.-H. Chang) would like to thank Prof. J.-H. Yu (ITP,
CAS) for helpful discussions and suggestion.

##  Appendix A QCD corrections to the relevant dimension-9 oprators

The essential process $d+d\to u+u+e+e$ for the $0\nu 2\beta$ decays involves
six ‘fermion legs’, so the effective Lagrangian for the decays is composed
generally by a set of independent dimension-9 operators as follows Pas:2000vn
:

$\displaystyle\mathcal{L}_{\rm eff}^{\rm
DBD}=\frac{G_{F}^{2}\cos^{2}\theta_{C}}{2m_{p}}\sum_{X,Y,Z}\Big{[}\sum_{i=1}^{3}C_{iZ}^{XY}(\mu)\cdot\mathcal{O}_{iZ}^{XY}(\mu)+\sum_{j=4}^{5}C_{j}^{XY}(\mu)\cdot\mathcal{O}_{j}^{XY}(\mu)\Big{]},$
(60)

where $\theta_{C}$ is the Cabibbo angle, $\mu\simeq 0.1\;$GeV is the energy
scale where the decays take place, and the independent dimension-9 operators
$\mathcal{O}_{iZ}^{XY}(\mu)$, $\mathcal{O}_{j}^{XY}(\mu)$ are defined as:

$\displaystyle{\mathcal{O}}_{1Z}^{XY}(\mu)=4(\bar{u}P_{X}d)(\bar{u}P_{Y}d)j_{Z},$
$\displaystyle{\mathcal{O}}_{2Z}^{XX}(\mu)=4(\bar{u}\sigma_{\mu\nu}P_{X}d)(\bar{u}\sigma^{\mu\nu}P_{X}d)j_{Z},$
$\displaystyle{\mathcal{O}}_{3Z}^{XY}(\mu)=4(\bar{u}\gamma_{\mu}P_{X}d)(\bar{u}\gamma^{\mu}P_{Y}d)j_{Z},$
$\displaystyle{\mathcal{O}}_{4}^{XY}(\mu)=4(\bar{u}\gamma_{\mu}P_{X}d)(\bar{u}\sigma^{\mu\nu}P_{Y}d)j_{\nu},$
$\displaystyle{\mathcal{O}}_{5}^{XY}(\mu)=4(\bar{u}\gamma_{\mu}P_{X}d)(\bar{u}P_{Y}d)j^{\mu},$
(61)

where $X,Y,Z=L,R$ ; $P_{R/L}=(1\pm\gamma^{5})/2$, the leptonic currents are
defined as

$\displaystyle j_{R/L}=\bar{e}(1\pm\gamma^{5})e^{c},\qquad
j_{\mu}=\bar{e}\gamma_{\mu}\gamma^{5}e^{c}.$ (62)

In the region from $M_{W_{L}}$ to $1.0\;$GeV for the energy scale $\mu$, the
pQCD is applicable, so the QCD corrections to the effective Lagrangian for the
$0\nu 2\beta$ decays can be carried out by renormalization group equation
method Buchalla:1995vs ; Buras:1998raa . The corresponding QCD corrections for
the $0\nu 2\beta$ decays were also calculated in Ref. Gonzalez:2015ady , so
here we describe how the corrections are calculated briefly. In RGE method the
renormalized operator matrix elements $<\mathcal{O}_{i}>^{\rm(R)}$
($\mathcal{O}_{i}$ is defined in Eq. (61)) for pQCD relate to their bare ones
up to one-loop level generally as the following form:

$\displaystyle<\mathcal{O}_{i}>^{\rm(R)}=\Big{[}\delta_{ij}+\frac{\alpha_{s}}{4\pi}b_{ij}(\frac{1}{\varepsilon}+\ln\frac{\mu^{2}}{-p^{2}})\Big{]}<\mathcal{O}_{j}>^{\rm
bare},$ (63)

where $<\mathcal{O}_{j}>^{\rm bare}$ are the ‘bare operator matrix elements’.
The ‘renormalization’ for quark fields $q$ and operator elements
$<\mathcal{O}_{j}>$ are given by

$\displaystyle q^{\rm bare}=Z_{q}^{1/2}q^{\rm
R}\;\;(q=u,d),\;\;\;<\mathcal{O}_{j}>^{\rm
bare}=Z_{q}^{2}Z_{ij}^{-1}<\mathcal{O}_{j}>^{\rm R},$ (64)

where

$\displaystyle
Z_{q}=1-C_{F}\frac{\alpha_{s}}{4\pi}\frac{1}{\varepsilon}+\mathcal{O}(\alpha_{s}^{2}),$
(65)

and $C_{F}=(N^{2}-1)/(2N)$ is the $SU(N)$ color factor ($N=3$). The
singularities in Eq. (63) are required to be cancelled, then we have

$\displaystyle
Z_{ij}Z_{q}^{-2}=\Big{[}\delta_{ij}+\frac{\alpha_{s}}{4\pi}b_{ij}\frac{1}{\varepsilon}\Big{]}.$
(66)

Then $Z_{ij}$ can be read out and written as

$\displaystyle
Z_{ij}=\delta_{ij}+\frac{\alpha_{s}}{4\pi}(b_{ij}-2C_{F}\delta_{ij})\frac{1}{\varepsilon}+\mathcal{O}(\alpha_{s}^{2}).$
(67)

Considering the 4-quark-leg operators $\mathcal{O}(q^{\rm bare})$ in the
effective Lagrangian Eq. (60) are constructed by the bare quark field $q^{\rm
bare}$, and the corresponding coefficients $C^{\rm bare}$ are also bare. Then
$q^{\rm bare}$, $C^{\rm bare}$ relate to the renormalized ones as

$\displaystyle q^{\rm bare}=Z_{q}^{1/2}q^{\rm R},\;\;\;\;C_{i}^{\rm
bare}=Z_{ij}^{C}C_{j}^{\rm R}.$ (68)

Hence we have

$\displaystyle C_{k}^{\rm bare}\mathcal{O}_{k}(q^{\rm
bare})=Z_{q}^{2}Z_{ij}^{C}C_{j}^{\rm R}\mathcal{O}_{i}(q^{\rm R}),$ (69)

and the matrix elements for the QCD corrections are read as

$\displaystyle\quad Z_{q}^{2}Z_{ij}^{C}C_{j}^{\rm R}<\mathcal{O}_{i}(q^{\rm
R})>^{\rm bare}=C_{j}^{\rm R}<\mathcal{O}_{i}(q^{\rm bare})>^{\rm R},$ (70)

Combining Eq. (70) and Eq. (64), we can obtain

$\displaystyle Z_{ij}^{C}=Z_{ij}^{-1}.$ (71)

Due to the fact that the bare quantities $C_{i}^{\rm bare}$ do not depend on
the renormalization energy scale $\mu$, we have

$\displaystyle\quad\;\frac{\rm d}{{\rm d}\ln\mu}C_{i}^{\rm bare}=\frac{\rm
d}{{\rm d}\ln\mu}Z_{ij}^{-1}C_{j}^{\rm R}=0,$ (72)

which can be rewritten as

$\displaystyle\frac{{\rm d}\vec{C}^{\rm R}(\mu)}{{\rm
d}\ln\mu}=\hat{\gamma}^{\rm T}\vec{C}^{\rm R}(\mu).$ (73)

Eq. (73) is the RGE accordingly for the Wilson coefficients, where
$\vec{C}=(C_{1},\;C_{2},\cdot\cdot\cdot)$ is written as a vector form. Then
the anomalous dimension matrix $\hat{\gamma}$ can be written as

$\displaystyle\hat{\gamma}=\frac{1}{\hat{Z}}\frac{\rm d}{{\rm
d}\ln\mu}\hat{Z}.$ (74)

Combining with the one-loop expression in the $\overline{\rm{MS}}$-scheme

$\displaystyle\hat{\gamma}(\alpha_{s})=-2\alpha_{s}\frac{\partial\hat{Z}_{1}(\alpha_{s})}{\partial\alpha_{s}},$
(75)

where $\hat{Z}_{1}$ is the coefficient matrix of $1/\varepsilon$ in Eq. (67),
the ‘anomalous dimension matrix’ to the leading order can be written as

$\displaystyle\gamma_{ij}(\alpha_{s})=\frac{\alpha_{s}}{4\pi}\gamma_{ij},\;\;\;{\rm
with}\;\;\;\gamma_{ij}=-2(b_{ij}-2C_{F}\delta_{ij}).$ (76)

Then solving Eq. (73), the evolution of Wilson coefficients $\vec{C}^{\rm
R}(\mu)$ can be expressed by $\vec{C}^{\rm R}(\Lambda)$ in terms of the
$\mu$-evolution matrix $\hat{U}(\mu,\Lambda)$:

$\displaystyle\vec{C}^{\rm R}(\mu)=\hat{U}(\mu,\Lambda)\cdot\vec{C}^{\rm
R}(\Lambda),$ (77)

where precisely

$\displaystyle\hat{U}(\mu,\Lambda)=\hat{V}{\rm
Diag}\Big{\\{}\Big{[}\frac{\alpha_{s}(\Lambda)}{\alpha_{s}(\mu)}\Big{]}^{\gamma_{i}/(2\beta_{0})}\Big{\\}}\hat{V}^{-1},$
(78)

and

$\displaystyle{\rm Diag}\\{\gamma_{i}\\}=\hat{V}^{-1}\hat{\gamma}\hat{V},$
(79)

with $\hat{\gamma}$ is $\gamma_{ij}$ in matrix form. The running coupling
constant to one-loop level of QCD can be written as

$\displaystyle\alpha_{s}(\mu)=\frac{\alpha_{s}(\Lambda)}{1-\beta_{0}\frac{\alpha_{s}(\Lambda)}{2\pi}\ln(\frac{\Lambda}{\mu})},$
(80)

with $\beta_{0}=(33-2f)/3$, and $f$ is the number of the active quark flavors
which is varied with the energy scale $\mu$, and only the quark $q^{f}$ with
mass $m_{f}$ smaller than the upper bound of the considered energy scale
region are ‘active’. Thus in the region from $\sim 1.0\;$GeV to $m_{c}\sim
1.3\;$GeV we have $f=3$, in the region from $m_{c}$ to $m_{b}\sim 4.6\;$GeV we
have $f=4$ and in the region from $m_{b}$ to $M_{W}$ we have $f=5$. As a
result, the required matrix to describe the QCD RGE evolution from $\mu=M_{W}$
to $\mu\simeq 1.0\;$GeV is

$\displaystyle\hat{U}(\mu,\Lambda=M_{W})=\hat{U}^{f=3}(\mu,\mu_{c})\hat{U}^{f=4}(\mu_{c},\mu_{b})\hat{U}^{f=5}(\mu_{b},M_{W}).$
(81)

To the leading order, the QCD corrections of the operators in Eq. (61)
correspond to Fig. 9 as follows.

Figure 9: One-loop QCD corrections to the dimension-9 operators for the $0\nu
2\beta$ decays in the effective Lagragian.

Calculating the diagrams respectively, the operator matrix elements
corresponding to those in Eq. (63) have the following structures

$\displaystyle{\rm Fig.~{}\ref{QCDplot}}(a1)\Rightarrow\mu^{4-D}\int\frac{{\rm
d}^{D}k}{(2\pi)^{D}}(\bar{u}\gamma_{\alpha}\frac{i}{k_{\beta}\gamma^{\beta}}\Gamma_{i}\frac{i}{k_{\eta}\gamma^{\eta}}\gamma^{\alpha}d)(\bar{u}\Gamma_{j}d)\cdot\frac{-i}{(k-p)^{2}}(ig_{3})^{2}C_{F}$
$\displaystyle\qquad\qquad\;\;=(\bar{u}\gamma_{\alpha}\gamma_{\beta}\Gamma_{i}\gamma^{\beta}\gamma^{\alpha}d)(\bar{u}\Gamma_{j}d)\cdot\frac{1}{4}C_{F}\frac{\alpha_{s}}{4\pi}(\frac{1}{\varepsilon}+\ln\frac{\mu^{2}}{-p^{2}}),$
(82) $\displaystyle{\rm
Fig.~{}\ref{QCDplot}}(a2)\Rightarrow\mu^{4-D}\int\frac{{\rm
d}^{D}k}{(2\pi)^{D}}(\bar{u}\Gamma_{i}d)(\bar{u}\gamma_{\alpha}\frac{i}{k_{\beta}\gamma^{\beta}}\Gamma_{j}\frac{i}{k_{\eta}\gamma^{\eta}}\gamma^{\alpha}d)\cdot\frac{-i}{(k-p)^{2}}(ig_{3})^{2}C_{F}$
$\displaystyle\qquad\qquad\;\;=(\bar{u}\Gamma_{i}d)(\bar{u}\gamma_{\alpha}\gamma_{\beta}\Gamma_{j}\gamma^{\beta}\gamma^{\alpha}d)\cdot\frac{1}{4}C_{F}\frac{\alpha_{s}}{4\pi}(\frac{1}{\varepsilon}+\ln\frac{\mu^{2}}{-p^{2}}),$
(83) $\displaystyle{\rm
Fig.~{}\ref{QCDplot}}(b1)\Rightarrow\mu^{4-D}\int\frac{{\rm
d}^{D}k}{(2\pi)^{D}}(\bar{u}\Gamma_{i}\frac{i}{k_{\eta}\gamma^{\eta}}\gamma_{\alpha}T^{a}d)(\bar{u}\Gamma_{j}\frac{i}{-k_{\sigma}\gamma^{\sigma}}\gamma^{\alpha}T^{a}d)\cdot\frac{-i}{(k-p)^{2}}(ig_{3})^{2}$
$\displaystyle\qquad\qquad\;\;=-(\bar{u}\Gamma_{i}\gamma_{\sigma}\gamma_{\alpha}T^{a}d)(\bar{u}\Gamma_{j}\gamma^{\sigma}\gamma^{\alpha}T^{a}d)\cdot\frac{1}{4}\frac{\alpha_{s}}{4\pi}(\frac{1}{\varepsilon}+\ln\frac{\mu^{2}}{-p^{2}}),$
(84) $\displaystyle{\rm
Fig.~{}\ref{QCDplot}}(b2)\Rightarrow\mu^{4-D}\int\frac{{\rm
d}^{D}k}{(2\pi)^{D}}(\bar{u}\gamma_{\alpha}\frac{i}{k_{\eta}\gamma^{\eta}}\Gamma_{i}T^{a}d)(\bar{u}\gamma^{\alpha}\frac{i}{-k_{\sigma}\gamma^{\sigma}}\Gamma_{j}T^{a}d)\cdot\frac{-i}{(k-p)^{2}}(ig_{3})^{2}$
$\displaystyle\qquad\qquad\;\;=-(\bar{u}\gamma_{\alpha}\gamma_{\sigma}\Gamma_{i}T^{a}d)(\bar{u}\gamma^{\alpha}\gamma^{\sigma}\Gamma_{j}T^{a}d)\cdot\frac{1}{4}\frac{\alpha_{s}}{4\pi}(\frac{1}{\varepsilon}+\ln\frac{\mu^{2}}{-p^{2}}),$
(85) $\displaystyle{\rm
Fig.~{}\ref{QCDplot}}(c1)\Rightarrow\mu^{4-D}\int\frac{{\rm
d}^{D}k}{(2\pi)^{D}}(\bar{u}\Gamma_{i}\frac{i}{k_{\eta}\gamma^{\eta}}\gamma_{\alpha}T^{a}d)(\bar{u}\gamma^{\alpha}\frac{i}{k_{\sigma}\gamma^{\sigma}}\Gamma_{j}T^{a}d)\cdot\frac{-i}{(k-p)^{2}}(ig_{3})^{2}$
$\displaystyle\qquad\qquad\;\;=(\bar{u}\Gamma_{i}\gamma_{\sigma}\gamma_{\alpha}T^{a}d)(\bar{u}\gamma^{\alpha}\gamma^{\sigma}\Gamma_{j}T^{a}d)\cdot\frac{1}{4}\frac{\alpha_{s}}{4\pi}(\frac{1}{\varepsilon}+\ln\frac{\mu^{2}}{-p^{2}}),$
(86) $\displaystyle{\rm
Fig.~{}\ref{QCDplot}}(c2)\Rightarrow\mu^{4-D}\int\frac{{\rm
d}^{D}k}{(2\pi)^{D}}(\bar{u}\gamma_{\alpha}\frac{i}{k_{\eta}\gamma^{\eta}}\Gamma_{i}T^{a}d)(\bar{u}\Gamma_{j}\frac{i}{k_{\sigma}\gamma^{\sigma}}\gamma^{\alpha}T^{a}d)\cdot\frac{-i}{(k-p)^{2}}(ig_{3})^{2}$
$\displaystyle\qquad\qquad\;\;=(\bar{u}\gamma_{\alpha}\gamma_{\sigma}\Gamma_{i}T^{a}d)(\bar{u}\Gamma_{j}\gamma^{\sigma}\gamma^{\alpha}T^{a}d)\cdot\frac{1}{4}\frac{\alpha_{s}}{4\pi}(\frac{1}{\varepsilon}+\ln\frac{\mu^{2}}{-p^{2}}),$
(87)

where $\Gamma_{i}$ are the Lorentz structures of the operators in Eq. (61),
and $T^{a}$ are the generators of $SU(N)$. Since the lepton sector is
irrelevant with the QCD corrections, so in the calculation the leptonic factor
in the operators is irrelevant. According to Eqs. (63$-$76), the anomalous
dimension matrix elements can be extracted from Eqs. (82$-$87).

Summarizing the obtained anomalous dimension matrix elements for
$\mathcal{O}_{1Z}^{XY},\;\mathcal{O}_{2Z}^{XX},\;\mathcal{O}_{3Z}^{XY}$,
$\mathcal{O}_{4}^{XY}$ and $\mathcal{O}_{5}^{XY}$, we have

$\displaystyle\hat{\gamma}_{(12)}^{XX}=-2\left(\begin{array}[]{c}6C_{F}-3,\;\;-\frac{1}{2N}+\frac{1}{4}\\\
-12-\frac{24}{N},\;\;-3-2C_{F}\end{array}\right),\hat{\gamma}_{(31)}^{XY}=-2\left(\begin{array}[]{c}-\frac{3}{N},\;\;-6\\\
\;\;\;\;0,\;\;\;\;\;6C_{F}\end{array}\right),$ (92)
$\displaystyle\gamma_{(3)}^{XX}=-2\Big{(}\frac{3}{N}-3\Big{)},$
$\displaystyle\hat{\gamma}_{(45)}^{XX}=-2\left(\begin{array}[]{c}-\frac{3}{2}-C_{F},\;\;-\frac{3}{2}i-\frac{3i}{N}\\\
-\frac{i}{2}+\frac{i}{N},\;\;\;\;3C_{F}-\frac{3}{2}\end{array}\right),\hat{\gamma}_{(45)}^{XY}=-2\left(\begin{array}[]{c}-\frac{3}{2}-C_{F},\;\;\frac{3}{2}i+\frac{3i}{N}\\\
\frac{i}{2}-\frac{i}{N},\;\;\;\;3C_{F}-\frac{3}{2}\end{array}\right).$ (97)

It is easy to realize that our results on the anomalous dimensions for the
matrix elements $\mathcal{O}_{4,5}^{XY}$ coincide with those in Ref.
Liao:2020roy , but do not coincide with those in Ref. Gonzalez:2015ady . Hence
the calculational details about the anomalous dimensions of
$\mathcal{O}_{4}^{XY},\;\mathcal{O}_{5}^{XY}$ are given below:

$\bullet\mathcal{O}_{4}^{XX}:$

$\displaystyle{\rm(a)-(1)}:\frac{1}{4}C_{F}(\bar{u}\gamma_{\alpha}\gamma_{\beta}\gamma^{\mu}\gamma^{\beta}\gamma^{\alpha}P_{X}d)(\bar{u}\sigma_{\mu\nu}P_{X}d)=C_{F}(\bar{u}\gamma^{\mu}P_{X}d)(\bar{u}\sigma_{\mu\nu}P_{X}d),$
(98)
$\displaystyle{\rm(a)-(2)}:\frac{1}{4}C_{F}(\bar{u}\gamma^{\mu}P_{X}d)(\bar{u}\gamma_{\alpha}\gamma_{\beta}\sigma_{\mu\nu}\gamma^{\beta}\gamma^{\alpha}P_{X}d)=0,$
(99)
$\displaystyle{\rm(b)-(1)}:-\frac{1}{4}(\bar{u}\gamma^{\mu}\gamma^{\sigma}\gamma^{\alpha}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}\gamma_{\sigma}\gamma_{\alpha}T^{a}P_{X}d)$
$\displaystyle\qquad\quad\;\;=-3(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}T^{a}P_{X}d)+3i(\bar{u}\gamma_{\nu}T^{a}P_{X}d)(\bar{u}T^{a}P_{X}d),$
(100)
$\displaystyle(b)-(2):-\frac{1}{4}(\bar{u}\gamma^{\sigma}\gamma^{\alpha}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\gamma_{\sigma}\gamma_{\alpha}\sigma_{\mu\nu}T^{a}P_{X}d)$
$\displaystyle\qquad\quad\;\;=-(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}T^{a}P_{X}d),$
(101)
$\displaystyle(c)-(1):\frac{1}{4}(\bar{u}\gamma^{\mu}\gamma^{\sigma}\gamma^{\alpha}T^{a}P_{X}d)(\bar{u}\gamma_{\alpha}\gamma_{\sigma}\sigma_{\mu\nu}T^{a}P_{X}d)$
$\displaystyle\qquad\quad\;\;=3(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}T^{a}P_{X}d)+3i(\bar{u}\gamma_{\nu}T^{a}P_{X}d)(\bar{u}T^{a}P_{X}d),$
(102)
$\displaystyle(c)-(2):\frac{1}{4}(\bar{u}\gamma^{\alpha}\gamma^{\sigma}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}\gamma_{\sigma}\gamma_{\alpha}T^{a}P_{X}d)=(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}T^{a}P_{X}d),$
(103)
$\displaystyle\Rightarrow(a)+(b)+(c)=C_{F}(\bar{u}\gamma^{\mu}P_{X}d)(\bar{u}\sigma_{\mu\nu}P_{X}d)-3(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}T^{a}P_{X}d)$
$\displaystyle\qquad\qquad\qquad\qquad\;\;\;+3i(\bar{u}\gamma_{\nu}T^{a}P_{X}d)(\bar{u}T^{a}P_{X}d)-(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}T^{a}P_{X}d)$
$\displaystyle\qquad\qquad\qquad\qquad\;\;\;+3(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}T^{a}P_{X}d)+3i(\bar{u}\gamma_{\nu}T^{a}P_{X}d)(\bar{u}T^{a}P_{X}d)$
$\displaystyle\qquad\qquad\qquad\qquad\;\;\;+(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}T^{a}P_{X}d)$
$\displaystyle\qquad\qquad\qquad\qquad=(C_{F}-\frac{3}{2})\mathcal{O}_{4}^{XX}+(-\frac{3}{2}i-\frac{3i}{N})\mathcal{O}_{5}^{XX}.$
(104)

$\bullet\mathcal{O}_{4}^{XY}(X\neq Y):$

$\displaystyle{\rm(a)-(1)}:\frac{1}{4}C_{F}(\bar{u}\gamma_{\alpha}\gamma_{\beta}\gamma^{\mu}\gamma^{\beta}\gamma^{\alpha}P_{X}d)(\bar{u}\sigma_{\mu\nu}P_{Y}d)=C_{F}(\bar{u}\gamma^{\mu}P_{X}d)(\bar{u}\sigma_{\mu\nu}P_{Y}d),$
(105)
$\displaystyle{\rm(a)-(2)}:\frac{1}{4}C_{F}(\bar{u}\gamma^{\mu}P_{X}d)(\bar{u}\gamma_{\alpha}\gamma_{\beta}\sigma_{\mu\nu}\gamma^{\beta}\gamma^{\alpha}P_{Y}d)=0,$
(106)
$\displaystyle{\rm(b)-(1)}:-\frac{1}{4}(\bar{u}\gamma^{\mu}\gamma^{\sigma}\gamma^{\alpha}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}\gamma_{\sigma}\gamma_{\alpha}T^{a}P_{Y}d)$
$\displaystyle\qquad\quad\;\;=-(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}T^{a}P_{Y}d),$
(107)
$\displaystyle{\rm(b)-(2)}:-\frac{1}{4}(\bar{u}\gamma^{\sigma}\gamma^{\alpha}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\gamma_{\sigma}\gamma_{\alpha}\sigma_{\mu\nu}T^{a}P_{Y}d)$
$\displaystyle\qquad\quad\;\;=-3(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}T^{a}P_{Y}d)-3i(\bar{u}\gamma_{\nu}T^{a}P_{X}d)(\bar{u}T^{a}P_{Y}d),$
(108)
$\displaystyle{\rm(c)-(1)}:\frac{1}{4}(\bar{u}\gamma^{\mu}\gamma^{\sigma}\gamma^{\alpha}T^{a}P_{X}d)(\bar{u}\gamma_{\alpha}\gamma_{\sigma}\sigma_{\mu\nu}T^{a}P_{Y}d)=(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}T^{a}P_{Y}d),$
(109)
$\displaystyle{\rm(c)-(2)}:\frac{1}{4}(\bar{u}\gamma^{\alpha}\gamma^{\sigma}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}\gamma_{\sigma}\gamma_{\alpha}T^{a}P_{Y}d)$
$\displaystyle\qquad\quad\;\;=3(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}\sigma_{\mu\nu}T^{a}P_{Y}d)-3i(\bar{u}\gamma_{\nu}T^{a}P_{X}d)(\bar{u}T^{a}P_{Y}d),$
(110)
$\displaystyle\Rightarrow(a)+(b)+(c)=C_{F}\mathcal{O}_{4}^{XY}-4i(\bar{u}\gamma_{\nu}T^{a}P_{X}d)(\bar{u}T^{a}P_{Y}d)$
$\displaystyle\qquad\qquad\qquad\quad\;\;\;=(C_{F}-\frac{3}{2})\mathcal{O}_{4}^{XY}+(\frac{3}{2}i+\frac{3i}{N})\mathcal{O}_{5}^{XY}.$
(111)

$\bullet\mathcal{O}_{5}^{XX}:$

$\displaystyle{\rm(a)-(1)}:\frac{1}{4}C_{F}(\bar{u}\gamma_{\alpha}\gamma_{\beta}\gamma_{\mu}\gamma^{\beta}\gamma^{\alpha}P_{X}d)(\bar{u}P_{X}d)=C_{F}(\bar{u}\gamma^{\mu}P_{X}d)(\bar{u}P_{X}d),$
(112)
$\displaystyle{\rm(a)-(2)}:\frac{1}{4}C_{F}(\bar{u}\gamma_{\mu}P_{X}d)(\bar{u}\gamma_{\alpha}\gamma_{\beta}\gamma^{\beta}\gamma^{\alpha}P_{X}d)=4C_{F}(\bar{u}\gamma^{\mu}P_{X}d)(\bar{u}P_{X}d),$
(113)
$\displaystyle{\rm(b)-(1)}:-\frac{1}{4}(\bar{u}\gamma_{\mu}\gamma_{\sigma}\gamma_{\alpha}T^{a}P_{X}d)(\bar{u}\gamma^{\sigma}\gamma^{\alpha}T^{a}P_{X}d)$
$\displaystyle\qquad\quad\;\;=-i(\bar{u}\gamma^{\nu}T^{a}P_{X}d)(\bar{u}\sigma_{\nu\mu}T^{a}P_{X}d)-(\bar{u}\gamma_{\mu}T^{a}P_{X}d)(\bar{u}T^{a}P_{X}d),$
(114)
$\displaystyle{\rm(b)-(2)}:-\frac{1}{4}(\bar{u}\gamma_{\sigma}\gamma_{\alpha}\gamma_{\mu}T^{a}P_{X}d)(\bar{u}\gamma^{\sigma}\gamma^{\alpha}T^{a}P_{X}d)=-(\bar{u}\gamma_{\mu}T^{a}P_{X}d)(\bar{u}T^{a}P_{X}d),$
(115)
$\displaystyle{\rm(c)-(1)}:\frac{1}{4}(\bar{u}\gamma_{\mu}\gamma_{\sigma}\gamma_{\alpha}T^{a}P_{X}d)(\bar{u}\gamma^{\alpha}\gamma^{\sigma}T^{a}P_{X}d)$
$\displaystyle\qquad\quad\;\;=(\bar{u}\gamma_{\mu}T^{a}P_{X}d)(\bar{u}T^{a}P_{X}d)-i(\bar{u}\gamma^{\nu}T^{a}P_{X}d)(\bar{u}\sigma_{\nu\mu}T^{a}P_{X}d),$
(116)
$\displaystyle{\rm(c)-(2)}:\frac{1}{4}(\bar{u}\gamma_{\alpha}\gamma_{\sigma}\gamma_{\mu}T^{a}P_{X}d)(\bar{u}\gamma^{\sigma}\gamma^{\alpha}T^{a}P_{X}d)=(\bar{u}\gamma_{\mu}T^{a}P_{X}d)(\bar{u}T^{a}P_{X}d),$
(117)
$\displaystyle\Rightarrow(a)+(b)+(c)=5C_{F}\mathcal{O}_{5}^{XX}-2i(\bar{u}\gamma^{\nu}T^{a}P_{X}d)(\bar{u}\sigma_{\nu\mu}T^{a}P_{X}d)$
$\displaystyle\qquad\qquad\qquad\quad\;\;\;=(5C_{F}-\frac{3}{2})\mathcal{O}_{5}^{XX}+(-\frac{i}{2}+\frac{i}{N})\mathcal{O}_{4}^{XX}.$
(118)

$\bullet\mathcal{O}_{5}^{XY}(X\neq Y):$

$\displaystyle{\rm(a)-(1)}:\frac{1}{4}C_{F}(\bar{u}\gamma_{\alpha}\gamma_{\beta}\gamma_{\mu}\gamma^{\beta}\gamma^{\alpha}P_{X}d)(\bar{u}P_{Y}d)=C_{F}(\bar{u}\gamma^{\mu}P_{X}d)(\bar{u}P_{X}d),$
(119)
$\displaystyle{\rm(a)-(2)}:\frac{1}{4}C_{F}(\bar{u}\gamma_{\mu}P_{X}d)(\bar{u}\gamma_{\alpha}\gamma_{\beta}\gamma^{\beta}\gamma^{\alpha}P_{Y}d)=4C_{F}(\bar{u}\gamma^{\mu}P_{X}d)(\bar{u}P_{Y}d),$
(120)
$\displaystyle{\rm(b)-(1)}:-\frac{1}{4}(\bar{u}\gamma_{\mu}\gamma_{\sigma}\gamma_{\alpha}T^{a}P_{X}d)(\bar{u}\gamma^{\sigma}\gamma^{\alpha}T^{a}P_{Y}d)=-(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}T^{a}P_{Y}d),$
(121)
$\displaystyle{\rm(b)-(2)}:-\frac{1}{4}(\bar{u}\gamma_{\sigma}\gamma_{\alpha}\gamma_{\mu}T^{a}P_{X}d)(\bar{u}\gamma^{\sigma}\gamma^{\alpha}T^{a}P_{Y}d)$
$\displaystyle\qquad\quad\;\;=i(\bar{u}\gamma^{\nu}T^{a}\gamma^{5}P_{X}d)(\bar{u}\sigma_{\nu\mu}T^{a}\gamma^{5}P_{Y}d)-(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}T^{a}P_{Y}d),$
(122)
$\displaystyle{\rm(c)-(1)}:\frac{1}{4}(\bar{u}\gamma_{\mu}\gamma_{\sigma}\gamma_{\alpha}T^{a}P_{X}d)(\bar{u}\gamma^{\alpha}\gamma^{\sigma}T^{a}P_{Y}d)=(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}T^{a}P_{Y}d),$
(123)
$\displaystyle{\rm(c)-(2)}:\frac{1}{4}(\bar{u}\gamma_{\alpha}\gamma_{\sigma}\gamma_{\mu}T^{a}P_{X}d)(\bar{u}\gamma^{\sigma}\gamma^{\alpha}T^{a}P_{Y}d)$
$\displaystyle\qquad\quad\;\;=i(\bar{u}\gamma^{\nu}T^{a}\gamma^{5}P_{X}d)(\bar{u}\sigma_{\nu\mu}T^{a}\gamma^{5}P_{Y}d)+(\bar{u}\gamma^{\mu}T^{a}P_{X}d)(\bar{u}T^{a}P_{Y}d),$
(124)
$\displaystyle\Rightarrow(a)+(b)+(c)=5C_{F}\mathcal{O}_{5}^{XX}+2i(\bar{u}\gamma^{\nu}T^{a}P_{X}d)(\bar{u}\sigma_{\nu\mu}T^{a}P_{X}d)$
$\displaystyle\qquad\qquad\qquad\quad\;\;\;=(5C_{F}-\frac{3}{2})\mathcal{O}_{5}^{XY}+(\frac{i}{2}-\frac{i}{N})\mathcal{O}_{4}^{XY}.$
(125)

The Fierz transformation formalisms used in above calculation read

$\displaystyle(P_{X})_{12}(P_{X})_{34}=\frac{1}{2}(P_{X})_{14}(P_{X})_{32}+\frac{1}{8}(\sigma_{\mu\nu}P_{X})_{14}(\sigma^{\mu\nu}P_{X})_{32},$
$\displaystyle(P_{X})_{12}(P_{Y})_{34}=\frac{1}{2}(\gamma^{\mu}P_{X})_{14}(\gamma_{\mu}P_{Y})_{32},$
$\displaystyle(\sigma_{\mu\nu}P_{X})_{12}(\sigma^{\mu\nu}P_{X})_{34}=6(P_{X})_{14}(P_{X})_{32}-\frac{1}{2}(\sigma_{\mu\nu}P_{X})_{14}(\sigma^{\mu\nu}P_{X})_{32},$
$\displaystyle(\gamma^{\mu}P_{X})_{12}(\gamma_{\mu}P_{X})_{34}=-(\gamma^{\mu}P_{X})_{14}(\gamma_{\mu}P_{X})_{32},$
$\displaystyle(\gamma^{\mu}P_{X})_{12}(\gamma_{\mu}P_{Y})_{34}=2(P_{X})_{14}(P_{Y})_{32},$
$\displaystyle(\gamma^{\mu}P_{X})_{12}(P_{X})_{34}=\frac{1}{2}(\gamma^{\mu}P_{X})_{14}(P_{X})_{32}-\frac{i}{2}(\gamma_{\nu}P_{X})_{14}(\sigma^{\nu\mu}P_{X})_{32},$
$\displaystyle(\gamma^{\mu}P_{X})_{12}(P_{Y})_{34}=\frac{1}{2}(\gamma^{\mu}P_{X})_{14}(P_{Y})_{32}+\frac{i}{2}(\gamma_{\nu}P_{X})_{14}(\sigma^{\nu\mu}P_{Y})_{32},$
$\displaystyle(\gamma_{\mu}P_{X})_{12}(\sigma^{\mu\nu}P_{X})_{34}=\frac{3i}{2}(\gamma^{\nu}P_{X})_{14}(P_{X})_{32}-\frac{1}{2}(\gamma_{\mu}P_{X})_{14}(\sigma^{\mu\nu}P_{X})_{32},$
$\displaystyle(\gamma_{\mu}P_{X})_{12}(\sigma^{\mu\nu}P_{Y})_{34}=-\frac{3i}{2}(\gamma^{\nu}P_{X})_{14}(P_{X})_{32}-\frac{1}{2}(\gamma_{\mu}P_{X})_{14}(\sigma^{\mu\nu}P_{X})_{32},$
(126)

and for the generators $T^{a}$ of $SU(N)$ we have

$\displaystyle
T^{a}_{ij}T^{a}_{jk}=\frac{N^{2}-1}{2N}\delta_{ik}=C_{F}\delta_{ik},\;\;\;T^{a}_{ij}T^{a}_{kl}=\frac{1}{2}(\delta_{il}\delta_{kj}-\frac{1}{N}\delta_{ij}\delta_{kl}).$
(127)

Combining Eqs. (79-81) with Eq. (97), we can compute out the numerical QCD RGE
evolution matrices from $\Lambda=M_{W_{L}}$ to $\mu\simeq 1.0\;$GeV as

$\displaystyle\hat{U}_{(12)}^{XX}(\mu,\Lambda)=\left(\begin{array}[]{c}1.96\;\;\quad
0.01\\\ -2.82\quad
0.45\end{array}\right),\;\;\hat{U}_{(31)}^{XY}(\mu,\Lambda)=\left(\begin{array}[]{c}0.87\quad-1.4\\\
0\qquad 2.97\end{array}\right),$ (132) $\displaystyle
U_{3}^{XX}(\mu,\Lambda)=0.76,$
$\displaystyle\hat{U}_{(45)}^{XX}(\mu,\Lambda)=\left(\begin{array}[]{c}0.68\qquad-0.24i\\\
-0.016i\quad
0.71\end{array}\right),\;\;\hat{U}_{(45)}^{XY}(\mu,\Lambda)=\left(\begin{array}[]{c}0.68\qquad
0.34i\\\ 0.023i\quad 1.4\end{array}\right).$ (137)

## References

  * (1) M. Tanabashi et al. [Particle Data Group], “Review of Particle Physics”, Phys. Rev. D 98, 030001 (2018).
  * (2) P. A. Dirac, “The fundamental equations of quantum mechanics”, Proc. Roy. Soc. Lond. A A109, 642-653 (1925).
  * (3) E. Majorana, “Teoria simmetrica dell’elettrone e del positrone”, Nuovo Cim. 14, 171-184 (1937).
  * (4) M. Agostini et al., “Background-free search for neutrinoless double-$\beta$ decay of 76Ge with GERDA”, Nature 544, 47 (2017).
  * (5) V. Guiseppe et al., “The Status and Initial Results of the MAJORANA DEMONSTRATOR Experiment”, AIP Conf. Proc. 1894, 020010 (2017).
  * (6) A. Gando, Y. Gando, T. Hachiya, A. Hayashi, S. Hayashida, H. Ikeda et al. [KamLAND-Zen], “”, Phys. Rev. Lett. 117, 082503 (2016).
  * (7) J. Albert et al. [EXO], “Search for Neutrinoless Double-Beta Decay with the Upgraded EXO-200 Detector”, Phys. Rev. Lett. 120, 072701 (2018).
  * (8) M. Agostini et al. [GERDA], “Final Results of GERDA on the Search for Neutrinoless Double-$\beta$ Decay,” Phys. Rev. Lett. 125, 252502 (2020).
  * (9) Zsigmond A J 2019 LEGEND: The future of neutrinoless double-beta decay search with germanium detectors, Talk presented at TAUP-2019 Toyama, Japan.
  * (10) F. Agostini et al. [DARWIN], “Sensitivity of the DARWIN observatory to the neutrinoless double beta decay of 136Xe”, Eur. Phys. J. C 80, 808 (2020).
  * (11) H. Pas, M. Hirsch, H. Klapdor-Kleingrothaus and S. Kovalenko, “Towards a superformula for neutrinoless double beta decay”, Phys. Lett. B 453, 194-198 (1999).
  * (12) M. González, S. G. Kovalenko and M. Hirsch, “QCD running in neutrinoless double beta decay: Short-range mechanisms”, Phys. Rev. D 93, 013017 (2016); M. González, S. G. Kovalenko and M. Hirsch, Phys. Rev. D 97, 099907 (2018)(E).
  * (13) Y. Liao, X. D. Ma and H. L. Wang, “Effective field theory approach to lepton number violating decays $K^{\pm}\rightarrow\pi^{\mp}l^{\pm}_{\alpha}l^{\pm}_{\beta}$: long-distance contribution”, JHEP 03, 120 (2020).
  * (14) H. Pas, M. Hirsch, H. Klapdor-Kleingrothaus and S. Kovalenko, “A Superformula for neutrinoless double beta decay. 2. The Short range part”, Phys. Lett. B 498, 35-39 (2001).
  * (15) V. Cirigliano, W. Dekens, J. de Vries, M. L. Graesser and E. Mereghetti, “Neutrinoless double beta decay in chiral effective field theory: lepton number violation at dimension seven”, JHEP 12, 082 (2017).
  * (16) V. Cirigliano, W. Dekens, J. de Vries, M. L. Graesser, E. Mereghetti, S. Pastore and U. Van Kolck, “New Leading Contribution to Neutrinoless Double-$\beta$ Decay”, Phys. Rev. Lett. 120, 202001 (2018).
  * (17) V. Cirigliano, W. Dekens, J. de Vries, M. L. Graesser and E. Mereghetti, “A neutrinoless double beta decay master formula from effective field theory”, JHEP 12, 097 (2018).
  * (18) W. Dekens, J. de Vries, K. Fuyuto, E. Mereghetti and G. Zhou, “Sterile neutrinos and neutrinoless double beta decay in effective field theory”, JHEP 06, 097 (2020).
  * (19) J. Barea, J. Kotila and F. Iachello, “$0\nu\beta\beta$ and $2\nu\beta\beta$ nuclear matrix elements in the interacting boson model with isospin restoration”, Phys. Rev. C 91, 034304 (2015).
  * (20) J. Suhonen, “Calculation of the Beta and Beta Beta decay observables of Ca-48 using QRPA with and without particle number projection”, J. Phys. G 19, 139-160 (1993).
  * (21) E. Caurier, F. Nowacki and A. Poves, “Beta beta decay and nuclear structure”, Int. J. Mod. Phys. E 16, 552-560 (2007).
  * (22) F. Simkovic, A. Faessler, V. Rodin, P. Vogel and J. Engel, “Anatomy of nuclear matrix elements for neutrinoless double-beta decay”, Phys. Rev. C 77, 045503 (2008).
  * (23) J. Menendez, A. Poves, E. Caurier and F. Nowacki, “Disassembling the Nuclear Matrix Elements of the Neutrinoless beta beta Decay”, Nucl. Phys. A 818, 139-151 (2009).
  * (24) J. Barea and F. Iachello, “Neutrinoless double-beta decay in the microscopic interacting boson model”, Phys. Rev. C 79, 044301 (2009).
  * (25) T. R. Rodriguez and G. Martinez-Pinedo, “Energy density functional study of nuclear matrix elements for neutrinoless $\beta\beta$ decay”, Phys. Rev. Lett. 105, 252503 (2010).
  * (26) J. Suhonen, “Exotic weak decays of atomic nuclei”, AIP Conf. Proc. 1488, 326-333 (2012).
  * (27) J. Barea, J. Kotila and F. Iachello, “Nuclear matrix elements for double-$\beta$ decay”, Phys. Rev. C 87, 014315 (2013).
  * (28) 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, 095016 (2020).
  * (29) L. Graf, F. F. Deppisch, F. Iachello and J. Kotila, “Short-Range Neutrinoless Double Beta Decay Mechanisms”, Phys. Rev. D 98, 095023 (2018).
  * (30) F. F. Deppisch, M. Hirsch and H. Pas, “Neutrinoless Double Beta Decay and Physics Beyond the Standard Model”, J. Phys. G 39, 124007 (2012).
  * (31) J. Lopez-Pavon, S. Pascoli and C. f. Wong, “Can heavy neutrinos dominate neutrinoless double beta decay?”, Phys. Rev. D 87, 093007 (2013).
  * (32) R. L. Awasthi, M. Parida and S. Patra, “Neutrino masses, dominant neutrinoless double beta decay, and observable lepton flavor violation in left-right models and SO(10) grand unification with low mass $W_{R},Z_{R}$ bosons”, JHEP 08, 122 (2013).
  * (33) T. Asaka, S. Eijima and H. Ishida, “On neutrinoless double beta decay in the $\nu$MSM”, Phys. Lett. B 762, 371-375 (2016).
  * (34) M. Mitra, S. Pascoli and S. Wong, “Effect of cancellation in neutrinoless double beta decay”, Phys. Rev. D 90, 093005 (2014).
  * (35) S. Dell’Oro, S. Marcocci, M. Viel and F. Vissani, “Neutrinoless double beta decay: 2015 review”, Adv. High Energy Phys. 2016, 2162659 (2016).
  * (36) J. Engel and J. Menendez, “Status and Future of Nuclear Matrix Elements for Neutrinoless Double-Beta Decay: A Review”, Rept. Prog. Phys. 80, 046301 (2017).
  * (37) M. J. Dolinski, A. W. Poon and W. Rodejohann, “Neutrinoless Double-Beta Decay: Status and Prospects”, Ann. Rev. Nucl. Part. Sci. 69, 219-251 (2019).
  * (38) P. S. B. Dev, R. N. Mohapatra and Y. Zhang, “Probing the Higgs Sector of the Minimal Left-Right Symmetric Model at Future Hadron Colliders”, JHEP 05, 174 (2016).
  * (39) S. Patra, F. S. Queiroz and W. Rodejohann, “Stringent Dilepton Bounds on Left-Right Models using LHC data”, Phys. Lett. B 752, 186-190 (2016).
  * (40) M. Mitra, R. Ruiz, D. J. Scott and M. Spannowsky, “Neutrino Jets from High-Mass $W_{R}$ Gauge Bosons in TeV-Scale Left-Right Symmetric Models”, Phys. Rev. D 94, 095016 (2016).
  * (41) A. Maiezza, G. Senjanovic and J. C. Vasquez, “Higgs sector of the minimal left-right symmetric theory”, Phys. Rev. D 95, 095004 (2017).
  * (42) S. Khalil and H. Okada, “Dark Matter in B-L Extended MSSM Models”, Phys. Rev. D 79, 083510 (2009).
  * (43) A. Elsayed, S. Khalil and S. Moretti, “Higgs Mass Corrections in the SUSY B-L Model with Inverse Seesaw”, Phys. Lett. B 715, 208-213 (2012).
  * (44) S. Khalil and S. Moretti, “The $B-L$ Supersymmetric Standard Model with Inverse Seesaw at the Large Hadron Collider”, Rept. Prog. Phys. 80, 036201 (2017).
  * (45) L. Delle Rose, S. Khalil, S. J. D. King, C. Marzo, S. Moretti and C. S. Un, “Naturalness and dark matter in the supersymmetric B-L extension of the standard model”, Phys. Rev. D 96, 055004 (2017).
  * (46) J. L. Yang, T. F. Feng, S. M. Zhao, R. F. Zhu, X. Y. Yang and H. B. Zhang, “Two loop electroweak corrections to $\bar{B}\rightarrow X_{s}\gamma$ and $B_{s}^{0}\rightarrow\mu^{+}\mu^{-}$ in the B-LSSM”, Eur. Phys. J. C 78, 714 (2018).
  * (47) J. L. Yang, T. F. Feng, H. B. Zhang, G. Z. Ning and X. Y. Yang, “Top quark decays with flavor violation in the B-LSSM”, Eur. Phys. J. C 78, 438 (2018).
  * (48) J. L. Yang, T. F. Feng, Y. L. Yan, W. Li, S. M. Zhao and H. B. Zhang, “Lepton-flavor violation and two loop electroweak corrections to $(g-2)_{\mu}$ in the B-L symmetric SSM”, Phys. Rev. D 99, 015002 (2019).
  * (49) J. L. Yang, T. F. Feng, S. K. Cui, C. X. Liu, W. Li and H. B. Zhang, “Electric dipole moments of neutron and heavy quarks in the B-LSSM”, JHEP 04, 013 (2020).
  * (50) J. L. Yang, T. F. Feng and H. B. Zhang, “Electroweak baryogenesis and electron EDM in the B-LSSM”, Eur. Phys. J. C 80, 210 (2020).
  * (51) J. L. Yang, T. F. Feng and H. B. Zhang, “Electron and muon $(g-2)$ in the B-LSSM”, J. Phys. G 47, 055004 (2020).
  * (52) S. Weinberg, “Baryon and Lepton Nonconserving Processes”, Phys. Rev. Lett. 43, 1566-1570 (1979).
  * (53) R. Mohapatra and A. Smirnov, “Neutrino Mass and New Physics”, Ann. Rev. Nucl. Part. Sci. 56, 569-628 (2006).
  * (54) R. Mohapatra et al., “Theory of neutrinos: A White paper”, Rept. Prog. Phys. 70, 1757-1867 (2007).
  * (55) R. N. Mohapatra and G. Senjanovic, “Neutrino Mass and Spontaneous Parity Nonconservation”, Phys. Rev. Lett. 44, 912 (1980).
  * (56) R. N. Mohapatra and J. D. Vergados, “A New Contribution to Neutrinoless Double Beta Decay in Gauge Models”, Phys. Rev. Lett. 47, 1713-1716 (1981).
  * (57) C. E. Picciotto and M. S. Zahir, “Neutrinoless Double Beta Decay in Left-right Symmetric Models”, Phys. Rev. D 26, 2320 (1982).
  * (58) M. Hirsch, H. Klapdor-Kleingrothaus and O. Panella, “Double beta decay in left-right symmetric models”, Phys. Lett. B 374, 7-12 (1996).
  * (59) V. Tello, M. Nemevsek, F. Nesti, G. Senjanovic and F. Vissani, “Left-Right Symmetry: from LHC to Neutrinoless Double Beta Decay”, Phys. Rev. Lett. 106, 151801 (2011).
  * (60) J. Chakrabortty, H. Devi, S. Goswami and S. Patra, “Neutrinoless double-$\beta$ decay in TeV scale Left-Right symmetric models”, JHEP 08, 008 (2012).
  * (61) J. Barry and W. Rodejohann, “Lepton number and flavour violation in TeV-scale left-right symmetric theories with large left-right mixing”, JHEP 09, 153 (2013).
  * (62) P. S. Bhupal Dev, S. Goswami, M. Mitra and W. Rodejohann, “Constraining Neutrino Mass from Neutrinoless Double Beta Decay”, Phys. Rev. D 88, 091301 (2013).
  * (63) C. H. Lee, P. S. Bhupal Dev and R. Mohapatra, “Natural TeV-scale left-right seesaw mechanism for neutrinos and experimental tests”, Phys. Rev. D 88, 093010 (2013).
  * (64) W. Huang and J. Lopez-Pavon, “On neutrinoless double beta decay in the minimal left-right symmetric model”, Eur. Phys. J. C 74, 2853 (2014).
  * (65) P. S. Bhupal Dev, S. Goswami and M. Mitra, “TeV Scale Left-Right Symmetry and Large Mixing Effects in Neutrinoless Double Beta Decay”, Phys. Rev. D 91, 113004 (2015).
  * (66) F. Ahmed and M. Horoi, “Interference effects for $0\nu\beta\beta$ decay in the left-right symmetric model”, Phys. Rev. C 101, 035504 (2020).
  * (67) D. Borah and A. Dasgupta, “Neutrinoless Double Beta Decay in Type I+II Seesaw Models”, JHEP 11, 208 (2015).
  * (68) S. Ge, M. Lindner and S. Patra, “New physics effects on neutrinoless double beta decay from right-handed current”, JHEP 10, 077 (2015).
  * (69) R. L. Awasthi, P. S. Bhupal Dev and M. Mitra, “Implications of the Diboson Excess for Neutrinoless Double Beta Decay and Lepton Flavor Violation in TeV Scale Left Right Symmetric Model”, Phys. Rev. D 93, 011701 (2016).
  * (70) G. Bambhaniya, P. S. B. Dev, S. Goswami and M. Mitra, “The Scalar Triplet Contribution to Lepton Flavour Violation and Neutrinoless Double Beta Decay in Left-Right Symmetric Model”, JHEP 04, 046 (2016).
  * (71) P. Pritimita, N. Dash and S. Patra, “Neutrinoless Double Beta Decay in LRSM with Natural Type-II seesaw Dominance”, JHEP 10, 147 (2016).
  * (72) F. F. Deppisch, C. Hati, S. Patra, P. Pritimita and U. Sarkar, “Neutrinoless double beta decay in left-right symmetric models with a universal seesaw mechanism”, Phys. Rev. D 97, 035005 (2018).
  * (73) J. Kotila, J. Ferretti and F. Iachello, “Long-range neutrinoless double beta decay mechanisms”, [arXiv:2110.09141 [hep-ph]].
  * (74) J. Kotila, “Comparison of Microscopic Interacting Boson Model and Quasiparticle Random Phase Approximation $0\nu\beta\beta$ Decay Nuclear Matrix Elements”, Front. Astron. Space Sci. 8, 42 (2021).
  * (75) S. Bertolini, A. Maiezza and F. Nesti, “Present and Future K and B Meson Mixing Constraints on TeV Scale Left-Right Symmetry”, Phys. Rev. D 89, 095028 (2014).
  * (76) T. Tomoda, “Double beta decay”, Rept. Prog. Phys. 54, 53-126 (1991).
  * (77) I. Esteban, M. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni and T. Schwetz, “Dynamics of Phase Separation from Holography”, JHEP 01, 106 (2019).
  * (78) A. M. Sirunyan et al. [CMS], “Search for third-generation scalar leptoquarks and heavy right-handed neutrinos in final states with two tau leptons and two jets in proton-proton collisions at $\sqrt{s}=13$ TeV”, JHEP 07, 121 (2017).
  * (79) A. M. Sirunyan et al. [CMS], “Search for a heavy right-handed W boson and a heavy neutrino in events with two same-flavor leptons and two jets at $\sqrt{s}=$ 13 TeV”, JHEP 05, 148 (2018).
  * (80) M. Aaboud et al. [ATLAS], “Search for heavy Majorana or Dirac neutrinos and right-handed $W$ gauge bosons in final states with two charged leptons and two jets at $\sqrt{s}=13$ TeV with the ATLAS detector”, JHEP 01, 016 (2019).
  * (81) M. Aaboud et al. [ATLAS], “Search for a right-handed gauge boson decaying into a high-momentum heavy neutrino and a charged lepton in $pp$ collisions with the ATLAS detector at $\sqrt{s}=13$ TeV”, Phys. Lett. B 798, 134942 (2019).
  * (82) G. Buchalla, A. J. Buras and M. E. Lautenbacher, “Weak decays beyond leading logarithms”, Rev. Mod. Phys. 68, 1125-1144 (1996).
  * (83) A. J. Buras, “Weak Hamiltonian, CP violation and rare decays”, [arXiv:hep-ph/9806471 [hep-ph]].