\section{Introduction}
\label{sec:intro}
Neutrino mass and mixing, which was confirmed by oscillation experiments \cite{SNO:2002tuh,Super-Kamiokande:2016yck,T2K:2019efw,DayaBay:2012fng,DoubleChooz:2011ymz} can not be understood within the Standard Model (SM) of particle physics since it predicts massless neutrinos. So, there has to be a mechanism beyond the SM which generates nonzero mass for these tiny particles. The seesaw mechanism has become quite famous for explaining the same by extending the SM in the minimal possible way. Some of the variants of this mechanism are type-I \cite{Minkowski:1977sc, Mohapatra:1979ia, Yanagida:1979as, GellMann:1980vs}, type-II\cite{Magg:1980ut, Schechter:1980gr, Cheng:1980qt, Lazarides:1980nt, Mohapatra:1980yp} and type-III  \cite{Foot:1988aq,He:2012ub} seesaw which can be achieved by adding a right-handed neutrino, a scalar triplet and a fermion triplet to the SM respectively. However, a heavy right-handed scale associated with these seesaw mechanisms renders them unverifiable at the collider experiments. Thus arises the necessity of bringing down the seesaw scale to a verifiable TeV range.
The seesaw mechanisms assume neutrinos are Majorana particles, which can be probed via the lepton number violating process of neutrinoless double beta decay \cite{PhysRevD.25.2951}. Such a rare transition occurs when two neutrons simultaneously decay into two protons and two electrons without any neutrinos. It can be induced either by light left-handed neutrinos, called the standard mechanism or by exotic particles like heavy right-handed neutrinos or sterile neutrinos, called new physics contribution. In the standard mechanism case, the experimental limits on the half-life of the decay can only be saturated by quasi-degenerate \cite{Bilenky:2001rz} light neutrinos, which are disfavored by cosmological data sets \cite{Planck:2018vyg,Abazajian:2022ofy,RoyChoudhury:2019hls,RoyChoudhury:2018gay}. On the other hand, identifying the correct neutrino mass hierarchy, considering the sum of light neutrino masses, would require a multi-ton scale detector that is beyond feasible in the near future. Any comparative future experimental observation of $0\nu\beta\beta$ decay would only be attributed to new physics contribution. 
The current lower limit on the decay half life of $\text{Ge}^{76}$ is $T^{0\nu}_{1/2} > 1.8 \times 10^{26}$ yrs at $90\%$ C.L. from GERDA~\cite{GERDA:2020xhi}. Experiments using the isotope $\text{Xe}^{136}$ like EXO-200~\cite{,EXO-200:2019rkq} and KamLAND-Zen~\cite{KamLAND-Zen:2016pfg,KamLAND-Zen:2022tow} have derived the lower bounds on half-life as $T^{0\nu}_{1/2} > 3.5 \times 10^{25}$ yrs and $T^{0\nu}_{1/2} > 1.07 \times 10^{26}$ yrs respectively. With this motivation, we consider a Left-Right symmetric model with a double seesaw mechanism \cite{PhysRevLett.56.561,PhysRevD.34.1642} as new physics and study the new contributions to $0\nu\beta\beta$ decay process.

Left-Right Symmetric Model (LRSM) \cite{Mohapatra:1974gc,Pati:1974yy,Senjanovic:1975rk,Senjanovic:1978ev} is a well-suited candidate for physics beyond SM for several reasons. To name a few, it can explain the theoretical origin of maximal parity violation in weak interaction, it can incorporate neutrino mass due to the presence of a right-handed neutrino state, it appears as a subgroup of SO(10) Grand Unified Theory, and it can be broken down to SM gauge symmetry at low energies. Moreover, it delivers rich phenomenology if the left-right symmetry breaking occurs at few TeV scale \cite{Keung:1983uu,  Ferrari:2000sp,  Schmaltz:2010xr, Nemevsek:2011hz, Chen:2011hc, Chakrabortty:2012pp, Das:2012ii, AguilarSaavedra:2012gf, Han:2012vk, Chen:2013fna, Rizzo:2014xma,  Deppisch:2014zta, Deppisch:2015qwa, Gluza:2015goa, Ng:2015hba, Patra:2015bga,Dobrescu:2015qna,PhysRevLett.115.211802, Brehmer:2015cia, Dev:2015pga,  Coloma:2015una, Deppisch:2015cua,Dev:2015kca, Mondal:2015zba, Aguilar-Saavedra:2015iew,    Lindner:2016lpp, Lindner:2016lxq, Mitra:2016kov, Anamiati:2016uxp,   Khachatryan:2014dka, Aad:2015xaa, Khachatryan:2016jqo}. The spontaneous symmetry breaking of LRSM to SM plays a vital role in generating neutrino mass through the seesaw mechanism. The seesaw scheme varies with the choice of scalars considered in the left-right model and regulates the associated phenomenology. In general, symmetry breaking can be done with the help of Higgs doublets or Higgs triplets or with the combination of both doublets and triplets. In the case of Higgs doublets, neutrinos don't get Majorana mass, and thus the model forbids any signatures of lepton number violation or lepton flavour violation. In the case of Higgs triplets, neutrino mass is generated via the type-I plus type-II seesaw mechanism. Even though Majorana mass is generated for light and heavy neutrinos in this case, the seesaw can't be probed by experiments considering the high scale associated with it.

The seesaw scale can be brought down to the TeV range in the case of a linear seesaw and inverse seesaw, some of which are discussed in ref.\cite{Hirsch:2009mx,Gu:2010xc,Dev:2009aw,Deppisch:2015cua,Humbert:2015yva,Parida:2012sq,Brdar:2018sbk,ThomasArun:2021rwf,Sahu:2022xkq,Ezzat:2021bzs}. However in case of a linear seesaw and inverse seesaw the light neutrinos are Majorana. In contrast, the heavy neutrinos are pseudo-Dirac, due to which the heavy neutrinos do not play a dominant role in lepton number violation. To study dominant new contributions to LNV and LFV decays, the Higgs and fermion sectors of LRSM have been extended in various refs \cite{Tello:2010am,Barry:2013xxa,BhupalDev:2013ntw,Hernandez:2021uxx,Brdar:2019fur}.

We explore here a double seesaw mechanism \cite{PhysRevLett.56.561,PhysRevD.34.1642} within a left-right symmetric model without Higgs triplets allowing significant lepton number violation and new physics contribution to neutrinoless double beta decay. We keep the scalar sector of the model minimal while adding only one sterile neutrino per generation in the fermion sector.
Even though the Higgs and fermion sectors are the same as in the case of a linear and inverse seesaw, the choice of bare Majorana masses for sterile neutrinos can induce large Majorana masses for heavy RH neutrino as well. The non-zero masses for RH neutrinos are generated through the double seesaw mechanism by implementing seesaw approximation twice. In the first step, the Majorana mass matrix and masses of the RH neutrinos are generated via the type-I seesaw mechanism. In this case, light neutrino mass becomes linearly dependent on a heavy sterile neutrino mass scale. This is how the double seesaw mainly differs from the canonical seesaw mechanism, where the light neutrino masses are inversely proportional to heavy RH neutrino masses. Another essential feature of our model is that we express mass relations between light and heavy Majorana neutrinos in terms of oscillation parameters and the lightest neutrino mass. Thus, it enables us to derive meaningful information about the absolute scale of the lightest neutrino mass and mass hierarchy from the new contributions to the neutrinoless double beta decay process by saturating the current experimental limits.

The plan of the paper can be summarized as follows.
In Section \ref{sec:lrsm}, we give a brief description of the left-right 
symmetric model with the double seesaw mechanism. In Section \ref{numassmixing} 
we explain the implementation of the double seesaw mechanism  
and the origin of Majorana masses for light and heavy right-handed (RH) 
and sterile neutrinos. The 
generation of the masses of gauge bosons associated with the 
$SU(2)_{\rm R}$ gauge group as well the constraints on their masses and 
on their mixing with the Standard Model gauge bosons are also considered 
in this Section. The general expression for the neutrinoless double beta 
decay half-life, including the new physics (i.e., the ``non-standard'') 
contributions is given in Section \ref{sec:0nubb}, in which we 
also, present and discuss briefly the nuclear matrix elements 
of the process and their current uncertainties.  
Detailed phenomenological analysis of the non-standard contributions 
together with numerical estimates of their magnitude are presented in 
Section \ref{sec:PA}. We also give predictions of the considered model 
for the neutrinoless double beta decay ``generalised'' effective Majorana mass 
and half-life accounting, in particular, for the uncertainties of the 
relevant nuclear matrix elements.  
Section \ref{sec:comments} contains brief comments on the 
potential lepton flavour violation and collider 
phenomenology of the considered model. 
Section \ref{sec:concl} contains a summary of the 
results obtained in the present study.
In Appendix \ref{app:lrdssm}, we give a detailed derivation of the 
the masses and mixing of the light and heavy Majorana neutrinos 
in the considered left-right symmetric model with a double seesaw 
mechanism of neutrino mass generation.

\section{Left-Right Symmetric Model with Double Seesaw}
\label{sec:lrsm}
Left-Right symmetric models were proposed with the motivation of restoring 
parity (or left-right) symmetry at a high scale \cite{Pati:1974yy,Mohapatra:1974gc,Senjanovic:1978ev,Senjanovic:1975rk}. Therefore, in the model the 
left- and right-handed fermion fields are  assigned to 
$SU(2)_L$ and $SU(2)_R$ doublets, respectively, which are related 
by a discrete symmetry. 
The complete gauge group, which is an extension of the SM 
gauge group can be written as:
\begin{align}
\mathcal{G}_{LR}\equiv SU(2)_L \times SU(2)_R \times U(1)_{B-L}\,,
\label{eq:LB-LRSM}
\end{align}
where $SU(3)_C$ is omitted for simplicity. The electric charge for any particle in this model is defined as 
\begin{align}
Q=T_{3L}+T_{3R}+\frac{B-L}{2}\,.
\end{align}
where $T_{3L}$ ($T_{3R}$) is the third component of the isospin 
associated with the $SU(2)_L$ ($SU(2)_R$) gauge group.
The model's fermion sector comprises all the Standard Model fermions plus a right-handed neutrino $N_R$. The fermion fields  with their respective quantum numbers can be written as follows:
\begin{eqnarray}
&&q_{L}=\begin{pmatrix}u_{L}\\
d_{L}\end{pmatrix}\equiv[2,1,1/3]\,, ~ q_{R}=\begin{pmatrix}u_{R}\\
d_{R}\end{pmatrix}\equiv[1,2,1/3]\,,\nonumber \\
&&\ell_{L}=\begin{pmatrix}\nu_{L}\\
e_{L}\end{pmatrix}\equiv[2,1,-1] \,, ~ \quad \ell_{R}=\begin{pmatrix} N_{R}\\
e_{R}\end{pmatrix}\equiv[1,2,-1] \,. \nonumber
\end{eqnarray}

The scalar sector is responsible for the spontaneous 
symmetry breaking of LRSM to SM 
and plays a crucial role in deciding the type of seesaw mechanism 
through which neutrino masses can be generated. The left-right symmetry 
breaking can be done either with the help of doublets $H_L$, $H_R$ or 
triplets $\Delta_L$, $\Delta_R$, or with the combination of both doublets and triplets. The next step of symmetry breaking, i.e., the breaking of SM symmetry 
to $U(1)_{em}$, is done with the help of the doublet  $\phi$ contained 
in the bidoublet $\Phi$. The doublets $H_L$, $H_R$ and the bidoublet $\Phi$ 
have the form,
\begin{eqnarray}
&&
 H_L=
 \begin{pmatrix} 
 h_L^+\\
 h_L^0
 \end{pmatrix}\equiv[2,1,1]  \,,\nonumber \\
&&
H_R=
\begin{pmatrix} 
h_R^+  \\
h_R^0
\end{pmatrix}\equiv[1,2,1] \,, \nonumber \\
&&
\Phi=
\begin{pmatrix} 
\phi_{1}^0     &  \phi_{2}^+ \\
\phi_{1}^-     &  \phi_{2}^0
\end{pmatrix} \equiv[2,2,0]\,,\nonumber \\
 &&
\end{eqnarray}
The symmetry breaking steps can be sketched as follows:
\begin{eqnarray*}
&\mbox{\bf \underline{\blue{Spontaneous symmetry breaking of LRSM:}}}& \\
&\hspace*{-2cm}  \begin{array}[t]{c} \pmb{SU(2)_L}\\ \{T_L, T_{3L} \} \\ g_L\end{array} \pmb{\times} 
\underbrace{\begin{array}[t]{c} \pmb{SU(2)_R}\\ \{T_R, T_{3R} \} \\ g_R\end{array} \pmb{\times} 
\begin{array}[t]{c} \pmb{U(1)_{B-L}}\\ {\rm B-L}    \\ g_{BL}\end{array}} \\
\\ 
&\hspace*{2.5cm} \downarrow  \langle H_R(1,2,1)\rangle \hspace*{0.2cm}  \hspace*{0.2cm} 
           & \\
&\hspace*{-2.0cm} \underbrace{\begin{array}[t]{c} \pmb{SU(2)_L}\\ \{T_L, T_{3L} \} \\ g\equiv g_L\end{array} \hspace*{0.2cm} \pmb{\times} \hspace*{0.2cm}
\begin{array}[t]{c} \pmb{U(1)_Y}\\ Y       \\ g^\prime \end{array}} & \\ 
&\hspace*{3.5cm} \downarrow  \langle \phi(1_{L},1/2_{Y})\rangle \subset \Phi(2_L,2_R,0_{B-L})\quad 
            & \\
&\hspace*{1.0cm}\begin{array}[t]{c} \pmb{U(1)_{\rm em}}\\ \mbox{(Q,~e) }\quad \end{array} 
\begin{array}[t]{c}  \\  \hspace*{0.0cm} \red{Q=T_{3L} +Y} \end{array} 
&
\end{eqnarray*}

The first step of symmetry breaking, i.e. 
$SU(2)_R \times U(1)_{B-L} \to U(1)_Y$ is achieved by assigning a non-zero 
 vacuum expectation value (VEV) $\langle H^0_R\rangle$ to the neutral 
component $h^0_R$ of $H_R$  as $v_R$.
The scale of this symmetry breaking 
determines the mass of the heavily charged and neutral gauge
bosons associated with the $SU(2)_R$ symmetry, ${W_R}$ and $Z_R$. 
The scalar doublet $H_L$ doesn't play any role but is present because of 
the left-right invariance. The electroweak symmetry breaking 
i.e., $SU(2)_L \times U(1)_Y \to U(1)_{em}$, is done by assigning 
 non-zero VEVs 
$\langle \phi^0_1 \rangle \equiv v_1$ and $\langle \phi^0_2 \rangle \equiv v_2$ 
to the neutral components of $\Phi$,
with $v_{SM} = \sqrt{v^2_1+v^2_2}\simeq$~246 GeV. 
The neutral components of the scalar bidoublet generate masses for the quarks 
and charged leptons via the following Yukawa Lagrangian:
\begin{eqnarray}
-\mathcal{L}_{Yuk} &\supset& \overline{q_{L}} \left[Y_1 \Phi + Y_2 \widetilde{\Phi} \right] q_R 
+\,\overline{\ell_{L}} \left[Y_3 \Phi + Y_4 \widetilde{\Phi} \right] \ell_R+ \mbox{h.c.}\,,
\label{eqn:LR-Yuk}
\end{eqnarray}
where $\widetilde{\Phi} = \sigma_2 \Phi^* \sigma_2$ and $\sigma_2$ 
is the second Pauli matrix. 
When the scalar bidoublet $\Phi$ acquires  non-zero VEVs,
\begin{equation}
\langle \Phi \rangle
=
\begin{pmatrix}
 v_1 & 0 \\
0 & v_2 
\end{pmatrix}\,,
\end{equation}
it gives masses to quarks and charged leptons in the following manner:
\begin{eqnarray}
&&  M_u =  Y_1 v_1 + Y_2 v_2\,, \quad \quad \nonumber \\
&&    M_d =  Y_1 v_2 + Y_2 v_1\,, \quad \quad \nonumber \\
&&  M_e =  Y_3 v_2 + Y_4 v_1 \,.\quad \quad
\end{eqnarray}
Here $M_u$ ($M_d$) and $M_e$ are the up-type (down-type) 
quark and charged lepton mass matrices.
The Lagrangian in Eq. (\ref{eqn:LR-Yuk}) also yields Dirac mass for 
the light neutrinos as
\begin{equation}
M^\nu_D\equiv M_D = Y_3 v_1 + Y_4 v_2 \,.
\end{equation}

In contrast to Yukawa couplings, which are complex, $v_1$ and $v_2$ are here 
assumed to be real. Typically, in the context of left-right symmetric theories, 
one investigates mainly the generation of Majorana neutrino masses. 
However, we note  that when $v_2 \ll v_1$ and $|Y_3| \ll |Y_4|$ one can 
have small Dirac neutrino masses.
In this scenario, the charged lepton and neutrino masses can be written as:
\begin{eqnarray}
M_e &\approx & Y_4 v_1\,, \\
M^\nu_D &=& v_1 \left( Y_3 + M_e  \frac{v_2}{|v_1|^2} \right)\,. 
\end{eqnarray}

\section{Neutrino Masses and Mixing}
\label{numassmixing}

In order to implement the double (or cascade) seesaw mechanism \cite{PhysRevLett.56.561,PhysRevD.34.1642} of neutrino mass generation within the manifest left-right symmetric model, we extend the fermion sector with the addition of one sterile neutrino $S_L \equiv[1,1,0] $ per generation. The relevant interaction Lagrangian $\mathcal{L}_{LRDSM}$ is given by:
\begin{eqnarray}
-\mathcal{L}_{LRDSM} &=& \mathcal{L}_{M_D}+\mathcal{L}_{M_{RS}}  + 
\mathcal{L}_{M_S}\,, 
\end{eqnarray}
where the individual terms can be expanded as follows.
\begin{itemize}
 \item $\mathcal{L}_{M_D}$ is the Dirac mass term connecting 
left-handed and right-handed neutrino fields $\nu_L-N_R$:
   \begin{eqnarray}
     \mathcal{L}_{M_D} &=& \sum_{\alpha, \beta} \overline{\nu_{\alpha L}} [M_D]_{\alpha \beta} N_{\beta R} \mbox{+ h.c.} \nonumber \\
     &&\subset \sum_{\alpha, \beta} \overline{\ell_{\alpha L}} \left( (Y_\ell)_{\alpha \beta} \Phi + (\widetilde{Y}_\ell)_{\alpha \beta} \tilde{\Phi} \right) 
     \ell_{\beta R} \mbox{+ h.c.}
     \end{eqnarray}
 \item  $\mathcal{L}_{M_{RS}}$  is another Dirac mass term 
connecting $N_R$ and $S_L$ and 
in the considered left-right symmetric theory
it has the form:
     \begin{eqnarray}
     \mathcal{L}_{M_{RS}} &=& \sum_{\alpha, \beta} \overline{S_{\alpha L}} [M_{RS}]_{\alpha \beta} N_{\beta R} \mbox{+ h.c.} \nonumber \\
    &&\subset \sum_{\alpha, \beta} \overline{S_{\alpha L}} (Y_{RS})_{\alpha \beta} H_R 
    \ell_{\beta R} \mbox{+ h.c.}
     \end{eqnarray}
 \item The bare Majorana mass term $\mathcal{L}_{M_S}$
for sterile neutrinos $S_L$ is given by:
      \begin{eqnarray}
      \mathcal{L}_{M_{S}} &=& \frac{1}{2} \sum_{\alpha, \beta} \overline{S^c_{\alpha R}} [M_{S}]_{\alpha \beta} S_{\beta L} \mbox{+ h.c.}
      \nonumber \\
      &&\subset \sum_{\alpha, \beta}  \frac{1}{2} (M_S)_{\alpha \beta} 
\overline{S^c_{\alpha R}} S_{\beta L} \mbox{+ h.c.}\,,
      \end{eqnarray}
\end{itemize}
where $S^c_{\alpha R} \equiv C(\overline{S_{\alpha L}})^T$, 
$C$ being the charge conjugation matrix  
($C^{-1} \gamma_\mu C = -\, \gamma^T_\mu$).


\subsection{The Double Seesaw Approximation}

After the spontaneous symmetry breaking, the complete $9 \times 9$ 
neutral fermion mass matrix in the flavor basis of $\left(\nu_L, N^c_R, S_L \right)$ can be written as:
\begin{equation}
 \mathcal{M}_{LRDSM}=
 \left[ 
\begin{array}{c  c} 
  \begin{array}{c c} 
     {\bf 0} & M_D \\ 
     M^T_D & {\bf 0}
  \end{array} & 
  \begin{array}{c} 
     {\bf 0} \\ M_{RS}
  \end{array} \\ 
  \begin{array}{c c} 
     \,~{\bf 0}  &\quad  M^T_{RS}
  \end{array} & 
  \begin{array}{c} 
     {M_S}
  \end{array} \\
 \end{array} 
\right]
\label{eqn:dss1}       
\end{equation}
 We assume in what follows that  $|M_D| \ll |M_{RS}| \ll |M_s|$. 
This allows us to apply to the mass matrix  $\mathcal{M}_{LRDSM}$  
twice the seesaw approximate block diagonalization 
procedure for getting the mass matrices of light and heavy neutrinos, 
as discussed below.

\begin{itemize}
\item {\bf First Seesaw Approximation:}
 We implement the first seesaw block diagonalization procedure 
on the lower right $6\times 6$ sub-matrix of  $\mathcal{M}_{LRDSM}$ 
as indicated below.
\begin{eqnarray}
&& \mathcal{M}^{}_{LRDSM} 
=
 \left[ 
\begin{array}{c | c} 
  \begin{array}{c} 
     {\bf 0}
  \end{array} & 
  \begin{array}{c c} 
     M_D & {\bf 0} 
  \end{array} \\ 
  \hline
  \begin{array}{c} 
     M^T_D  \\  {\bf 0}
  \end{array} & 
  \begin{array}{cc} 
     {\blue {0}} & \blue{M_{RS}} \\
     \blue{M^T_{RS}}  &  \blue{M_S}
  \end{array} \\
 \end{array} 
\right]  \nonumber \\
&&\hspace*{-0.2cm}
 \mathop{\xrightarrow{\hspace*{1.5cm}}}^{M_S > M_{RS} \gg M_D}_{\mbox{1st seesaw}} 
  \left[ 
\begin{array}{c | c} 
  \begin{array}{c c} 
     {\bf 0}
  \end{array} & 
  \begin{array}{c c} 
     M_D & {\bf 0} 
  \end{array} \\ 
  \hline
  \begin{array}{c} 
     M^T_D  \\  {\bf 0}
  \end{array} & 
  \begin{array}{cc} 
     \blue{- M_{RS} M^{-1}_S M^T_{RS}} & {\bf 0} \\
     {\bf 0}  &  \blue{M_S}
  \end{array} \\
 \end{array} 
\right]
\label{eqn:dss-b}  
\end{eqnarray}
\item {\bf Second Seesaw Approximation:-}
Denoting $- M_{RS} M^{-1}_S M^T_{RS}= M_R$, 
 which is the expression for the 
mass matrix for right-handed neutrinos, we repeat the diagonalization 
procedure with seesaw condition, $|{M_R}| \gg |M_D|$. 
We get the resultant matrix structure as
\begin{eqnarray}
 \left[ 
\begin{array}{c | c}  
  \begin{array}{c c} 
     {\bf 0} & \blue{M_D} \\ 
     \blue{M^T_D} & \blue{M_R}
  \end{array} &          \begin{array}{c} 
                       {\bf 0}  \\  {\bf 0}
                         \end{array}                \\
  \hline
  \begin{array}{c} 
   \hspace*{-0.0cm} {\bf 0}  \hspace*{0.6cm} {\bf 0}
  \end{array}          & {M_S} 
 \end{array} 
\right]
 \mathop{\xrightarrow{\hspace*{1.5cm}}}^{M_R \gg  M_D}_{\mbox{2nd seesaw}} 
 \left[ 
\begin{array}{c | c | c}  
     \blue{\tiny -M_D M^{-1}_R M^T_D} & 0  & 0\\
  \hline
  0  & \blue{M_{R}} & 0 \\
  \hline 
  0 & 0 & \blue{M_S}
 \end{array} 
\right] 
\label{eqn:dss-c}       
\end{eqnarray}
\end{itemize}
 Using the above results, the light neutrino, the heavy neutrino 
and sterile fermions mass matrices  
$m_\nu$,  $m_N$ and $m_S$  can be expressed as: 
\begin{eqnarray}
&&m_\nu \cong 
- M_D \left( -M_{RS} M^{-1}_S M^T_{RS} \right)^{-1} M^T_D\, \nonumber \\
  &&\hspace*{0.5cm}
=\frac{M_D}{M^T_{RS}}  M_S \frac{M^T_D}{M_{RS}}  ,\nonumber \\
&& m_N \equiv M_R \cong -M_{RS} M^{-1}_S M^T_{RS} ,\nonumber \\
&& m_S \cong M_S\, .
\label{massformulae}
\end{eqnarray}
\begin{table}[h!]
 \centering
	\begin{tabular}{|c|c|c||c|c|c||c|c|}
	\hline \hline
    $M_D$ & $M_{RS}$ & $M_S$ & $m_\nu$(eV) & $m_{N}$ & $m_{S}$ & $V^{\nu N}$ & $V^{N S}$ \\
	\hline
	$10^{-4}$  & $10^{3}$  & $10^{4}$ & $0.1$  & $10^{2}$ & $10^{4}$ & $10^{-6}$ & $0.1$ \\
    $10^{-5}$  & $10^{2}$  & $10^{3}$ & $0.01$  & $10$ & $10^{3}$ & $10^{-6}$ & $0.1$ \\    $10^{-5}$  & $10^{1}$  & $10^{2}$ & $0.1$  & $1$ & $10^{2}$ & $10^{-5}$ & $0.1$ \\
    \hline \hline
	\end{tabular}
	\label{tab:tab1}
\caption{A representative set of model parameters in left-right symmetric models and the order of magnitude estimation of various neutrino masses within the double seesaw mechanism. All the masses are expressed in units of GeV except the light neutrino masses, which are in the eV scale.}
\label{tab:dss}
\end{table}

In the double seesaw 
expression  for the light neutrino Majorana mass matrix 
as given in Eq. (\ref{massformulae}), different choices of $M_D$ and $M_{RS}$ 
 are possible. Following \cite{Smirnov:1993af,Altarelli:2004za,Lindner:2005pk,Ludl:2015tha,Bajc:2016eiw,Smirnov:2018luj}, 
we have considered in the present article the case of 
$M_D$ and $M_{RS}$ being proportional to identity such that  $M_D = k_d I$ 
and $M_{RS}= k_{rs}I$ , where $k_d$ and $k_{rs}$ are real 
constants with $|k_d| < |k_{rs}|$. 
This means, $ M_D M^{-1}_{RS} = \frac{k_d}{k_{rs}}I$.
As discussed in \cite{Lindner:2005pk,Brdar:2018sbk}, the equality and 
simultaneous diagonal structures of $M_D$ and $M_{RS}$  may arise as a 
consequence of $Z_2 \times Z_2$ symmetry \cite{Ludl:2015tha}. With the 
introduction of additional permutation symmetry in the diagonal elements of 
$M_D$ and $M_{RS}$, one can get equal diagonal elements. 
As we have indicated, these kinds of considerations have been reasoned 
for the double seesaw mechanism, e.g.,  
in the references \cite{Smirnov:1993af,Altarelli:2004za,Lindner:2005pk,Ludl:2015tha,Bajc:2016eiw,Smirnov:2018luj}.

 With the choices for the forms of $M_D$ 
and $M_{RS}$ made above, the relation between 
light neutrino  and sterile neutrino mass matrices $m_\nu$ and $m_S$
can be written as $ m_{\nu}= \frac{k^2_d}{k^2_{rs}} m_S $. 
The mass matrix $m_N $ can also be determined from Eq. (\ref{massformulae})
and the relationship between 
light neutrino  and heavy right-handed neutrino mass matrices 
$m_{\nu}$ and $m_N$ has the form $ m_{N}= - k^2_d\frac{1}{m_{\nu}}$ .

In the basis in which the charged lepton mass matrix is diagonal 
we will work with in what follows,
the light neutrino Majorana mass matrix is diagonalized with the help of a 
unitary mixing matrix -- the Pontecorvo, Maki, Nakagawa, Sakata (PMNS) 
mixing matrix 
$U_{\rm PMNS} \equiv U_{\nu}$ 
\cite{Pontecorvo:1957qd,Maki:1962mu,Pontecorvo:1967fh}:
$$m^{\rm diag}_\nu= U^\dagger_{\rm PMNS} m_\nu U^*_{\rm PMNS} 
= \mbox{diag}\left(m_1,m_2,m_3 \right)\,,~~m_i > 0\,,$$
so the physical masses $m_i$ are related to the mass matrix $m_\nu$ in 
the flavour basis as 
$$m_\nu = U_{\rm PMNS} m^{\rm diag}_\nu U^T_{\rm PMNS}\,.$$
 The right-handed neutrino Majorana mass matrix $m_N$ is diagonalized as 
$\widehat{m_N} = {U_N}^\dagger m_N {U_N}^*$, 
$\widehat{m_N} =  \mbox{diag}(m_{N_1},m_{N_2},m_{N_3})$, 
$m_{N_j}$ being the mass of the heavy RH Majorana neutrino $N_j$, $j=1,2,3$.
It proves convenient to work with positive masses 
of $N_j$, $m_{N_j} > 0$. 
Given the relation $ m_{N}= - k^2_d\frac{1}{m_{\nu}}$ and the 
positivity of the eigenvalues of $m_\nu$, 
the requirement that the eigenvalues of $m_N$ are also positive
implies that the unitary transformation matrices 
diagonalizing the light neutrino and the heavy right-handed neutrino 
mass matrices $m_\nu$ and $m_N$ are related in the following way:
\begin{eqnarray}
 U_N = i\,U^*_{\nu} \equiv i\,U^*_{PMNS}\,.
\label{eq:UN}
 \end{eqnarray} 
 Since $m_S = (k^2_{rs}/k^2_d)m_\nu$,              
the diagonalization of sterile neutrino Majorana mass matrix $m_S$, 
$\widehat{m_S} = {U_S}^\dagger m_S {U_S}^*$,
where $\widehat{m_S} =  \mbox{diag}(m_{S_1},m_{S_2},m_{S_3})$, $m_{S_k} > 0$, 
$k=1,2,3$, can be performed with the help of the same mixing matrix $U_{PMNS}$:
\begin{eqnarray}
 U_S = U_{\nu}\equiv U_{PMNS}\,. 
\label{eq:US}
 \end{eqnarray}
 Thus, in the considered 
scenario, 
the light neutrino masses $m_i$, the heavy RH neutrino masses 
$m_{N_j}$ and the sterile neutrino masses ($m_{S_k}$) are related  as follows:
\begin{equation}
m_{i} = \frac{k^2_d}{m_{N_i}} = \frac{k^2_d}{k^2_{rs}}\, m_{S_i}\,,~i=1,2,3\,.
\label{eq:massrel}
\end{equation}

 In what follows, we will use the standard parametrization of 
the PMNS matrix (see, e.g., \cite{ParticleDataGroup:2018ovx}):
\begin{eqnarray}
&&U_{\rm {PMNS}}= \nonumber \\
&& \hspace*{-0.3cm} 
\begin{pmatrix} c_{13}c_{12}&c_{13}s_{12}&s_{13}e^{-i\delta}\\
-c_{23}s_{12}-c_{12}s_{13}s_{23}e^{i\delta}&c_{12}c_{23}-s_{12}s_{13}s_{23}e^{i\delta}&s_{23}c_{13}\\
s_{12}s_{23}-c_{12}c_{23}s_{13}e^{i\delta}&-c_{12}s_{23}-s_{12}s_{13}c_{23}e^{i\delta}&c_{13}c_{23}
\end{pmatrix} \mbox{P}
\label{PMNS} 
\end{eqnarray}
where the mixing angles are denoted by $s_{ij}=\sin \theta_{ij}$, 
$c_{ij}=\cos \theta_{ij}$,  
 $0 \leq \theta_{ij}\leq \pi/2$,
$\delta$ is the Dirac CP violation phase, $0 \leq \delta \leq 2\pi$,
$\mbox{P} $ is the diagonal phase matrix 
containing the two Majorana CP violation phases
$\alpha$ and $\beta$ \cite{Bilenky:1980cx},
$\mbox{P}=\mbox{diag}\left(1, e^{i\alpha/2}, e^{i \beta/2} \right)$. 
The Majorana phases take values in the interval 
$[0, \pi]$. The experimental values of 
different oscillation parameters for 
the light neutrino mass spectrum 
with normal ordering (NO) and inverted ordering (IO) 
(see, e.g., \cite{ParticleDataGroup:2018ovx}) 
are taken from Ref.~\cite{PhysRevD.101.116013} and are 
presented in Table \ref{tab2_oscipara}. 
\begin{table}[]
\centering
\begin{tabular}{|l|c|c|}
\hline
Parameter & Best fit values & 3$\sigma$ range
\\
\hline\hline
$\Delta m^2_{21}\: [10^{-5}\,\mbox{eV}]$ & 7.34 & 6.92--7.90 \\
\hline
$|\Delta m^2_{31}|\: [10^{-3}\,\mbox{eV}]$ (NO) &  2.485 &  2.389--2.578\\
$|\Delta m^2_{32}|\: [10^{-3}\,\mbox{eV}]$ (IO)&  2.465  &  2.374--2.556 \\
\hline
$\sin^2\theta_{12} / 10^{-1}$ (NO) & 3.05 & 2.65--3.47\\
$\sin^2\theta_{12} / 10^{-1}$(IO) & 3.03 & 2.64--3.45\\
\hline
  $\sin^2\theta_{23} / 10^{-1}$ (NO)
	  &	5.45 & 4.36--5.95 \\
  $\sin^2\theta_{23} / 10^{-1}$ (I))
	  & 5.51 & 4.39--5.96 \\
\hline
$\sin^2\theta_{13} / 10^{-2}$ (NO) & 2.22 & 2.01--2.41 \\
$\sin^2\theta_{13} / 10^{-2}$ (IO) & 2.23 & 2.03--2.43 \\
    \hline
  \end{tabular}
  %
  \caption{The current updated estimates of experimental values of 
Neutrino oscillation parameters for global best fits and $3\sigma$ range taken from~\cite{PhysRevD.101.116013}.}
  \label{tab2_oscipara}  
\end{table}

\begin{itemize}
 \item \underline{Masses of light neutrinos}

 It proves convenient to express the masses of the two heavier 
neutrinos in terms of the lightest neutrino mass and the neutrino 
mass squared differences measured in neutrino oscillation experiments.
In the case of NO light neutrino mass spectrum, 
$m_1< m_2 < m_3$, we have: 
 \begin{eqnarray}
 &&m_1 = \mbox{lightest neutrino mass} \,,\qquad  \nonumber \\
 &&m_2 = \sqrt{m_1^2 +\Delta m_{\rm sol}^2}\,,\qquad \nonumber \\
 &&m_3 = \sqrt{m_1^2 +\Delta m_{\rm atm}^2}\,,
\label{eq1:NH_m2_m3}
\end{eqnarray}
where $\Delta m_{\rm sol}^2=\Delta m_{21}^2$ and 
$\Delta m_{\rm atm}^2=\Delta m_{31}^2$.
Similarly, for inverted mass ordering, $m_3 < m_1 < m_2$, we get: 
\begin{eqnarray}
 &&m_3 = \mbox{lightest neutrino mass} \,,\qquad  \nonumber \\
 &&m_1 = \sqrt{m_3^2 -\Delta m_{\rm sol}^2-\Delta m_{\rm atm}^2}\,,
\qquad \nonumber \\
 &&m_2 = \sqrt{m_3^2  -\Delta m_{\rm atm}^2 }\;, 
\label{eq1:IH_m1_m2}
\end{eqnarray}
with $\Delta m_{\rm sol}^2=\Delta m_{21}^2$ and 
$\Delta m_{\rm atm}^2=\Delta m_{32}^2$.

 Depending on the value of the lightest neutrino mass 
the neutrino mass spectrum can also be normal hierarchical (NH) 
when $m_1 \ll m_2 < m_3$, inverted hierarchical (IH) 
if $m_3 \ll m_1 < m_2$, or else quasi-degenerate (QD)
when $m_1 \cong m_2 \cong m_3$, $m^2_{1,2,3} >> \Delta m^2_{31(23)}$, 
i.e., $m_{1,2,3}\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  0.1$ eV.


 All considered types of neutrino mass spectrum are compatible with the existing data. The best upper limit on the lightest 
neutrino mass $m_{1(3)}$ has been obtained in the KATRIN experiment. It is in the range of the QD spectrum and effectively 
reads: $m_{1,2,3} < 0.80$ eV (90\% C.L.).


\item \underline {Masses of heavy RH neutrinos}
 
 Since $m_i$ and $m_{N_j}$ are inversely proportional to each other, 
for NO light neutrino mass spectrum, $m_1 < m_2 < m_3$, 
$m_{N_1}$ has to be the largest RH neutrino mass. 
We can express $m_{N_2}$ and $m_{N_3}$ in terms of $m_{N_1}$ 
and the light neutrino masses:
 \begin{eqnarray}
 m_{N_2} = \frac{m_{1}}{m_{2}} m_{N_1}\,, \quad
m_{N_3} = \frac{m_{1}}{m_{3}} m_{N_1}\,, \quad m_{N_3} <  m_{N_2} <  m_{N_1}\,.
\label{eqn:Mne1} 
\end{eqnarray}
In the case of IO spectrum, $m_3 < m_1 < m_2$, 
$m_{N_3}$ is the largest mass. The mass relations in this case become:
  \begin{eqnarray}
m_{N_1} = \frac{m_{3}}{m_{1}} m_{N_3}\,, \quad 
 m_{N_2} = \frac{m_{3}}{m_{2}} m_{N_3}\,, \quad  m_{N_1} <  m_{N_2} <  m_{N_3}\,. 
\label{eqn:Mne2} 
\end{eqnarray}


 \item \underline {Masses of sterile neutrinos}

 
 Since $m_i$ and $m_{S_k}$ are directly proportional to each other, 
in the NO case  $m_{S_3}$ is the heaviest sterile neutrino mass 
and the analogous mass relations read: 
\begin{eqnarray}
 m_{S_1} = \frac{m_{1}}{m_{3}} m_{S_3}\,, \quad 
 m_{S_2} = \frac{m_{2}}{m_{3}} m_{S_3}\,, \quad  m_{S_1} <  m_{S_2} < m_{S_3}\,.
\label{eqn:MSe1} 
\end{eqnarray}
For the IO spectrum we have $m_{S_3} <  m_{S_1} <  m_{S_2}$ and
\begin{eqnarray}
m_{S_1} = \frac{m_{1}}{m_{2}} m_{S_2}\,, \quad
m_{S_3} = \frac{m_{3}}{m_{2}} m_{S_2}\,, \quad  m_{S_3} <  m_{S_1} < m_{S_2}\,. 
\label{eqn:MSe2} 
\end{eqnarray}

\item \underline {Neutrino Mixing}
\end{itemize}

The diagonalisation of $\mathcal{M_{\rm LRDSM}}$ 
 leads to the following relation between 
the fields of the neutral fermions written in the flavour 
(weak interaction eigenstate) basis and in the mass eigenstate basis:
\begin{eqnarray}
&& \begin{pmatrix}
\nu_{\alpha L}\\ N^c_{\beta L}\\ S_{\gamma L}
\end{pmatrix}
= 
\begin{pmatrix}
V_{\alpha i}^{\nu \nu} & V_{\alpha j}^{\nu N} & V_{\alpha k}^{\nu S}\\V_{\beta i}^{ N \nu} & V_{\beta j}^{NN} & V_{\beta k}^{N S}\\V_{\gamma i}^{S \nu} & V_{\gamma j}^{S N} & V_{\gamma k}^{SS}
\end{pmatrix}               
\begin{pmatrix}
\nu_{iL}\\ N^c_{jL}\\ S_{kL}
\end{pmatrix} \,.
\label{nuNS}
\end{eqnarray}
 Here  $N^c_{\beta L} \equiv C(\overline{N_{\beta R}})^T$, 
 $N^c_{j L} \equiv C(\overline{N_{j R}})^T = N_{jL}$ 
($N_j$ are Majorana fields),
$C$ being the charge conjugation matrix, 
the indices $\alpha, \beta, \gamma$ run over three generations of 
light left-handed neutrinos, heavy right-handed neutrinos and 
sterile neutrinos in flavor basis respectively, 
whereas the indices $i,j,k$ run over corresponding mass eigenstates.
 The mixing matrix elements in Eq. (\ref{nuNS}) 
are given in Eq. (\ref{mixingmatrix1}) in 
 Appendix \ref{app:lrdssm}.
The mixing between the right-handed neutrinos and sterile neutrinos 
$(N^c_L-S_L)$ is given by the term,
\begin{equation}
V^{N S} \propto M_{RS} M^{-1}_S 
\end{equation}
while the mixing between the fields of the 
left-handed flavour neutrinos and the  
heavy right-handed neutrinos $(\nu_L-N^c_L)$ is determined by
\begin{equation}
V^{\nu N} \propto M_D M^{-1}_R = -M_D {M^T_{RS}}^{-1}  M_{S} M^{-1}_{RS} 
\end{equation}
The mixing between sterile and light neutrinos ($\nu_L - S_L$) is vanishing, 
$V_{\alpha k}^{\nu S}\cong 0$ and $V_{\gamma i}^{S \nu} \cong 0$.

The possible sets of numerical values for different mixing matrices, masses 
and mixing that can give rise to dominant contributions to LNV decays are 
listed in Table \ref{tab:dss}. By choosing one representative 
set of model parameters from 
the table, we get the mixing as given below. 
\begin{eqnarray}
\begin{pmatrix}
V_{\alpha i}^{\nu \nu} & V_{\alpha j}^{\nu N} & V_{\alpha k}^{\nu S}\\V_{\beta i}^{ N \nu} & V_{\beta j}^{NN} & V_{\beta k}^{N S}\\V_{\gamma i}^{S \nu} & V_{\gamma j}^{S N} & V_{\gamma k}^{SS}
\end{pmatrix}  \simeq \begin{pmatrix}
{\cal O}(1.0)       & {\cal O}(10^{-6})     & 0  \\
{\cal O}(10^{-6})      & {\cal O}(1.0)    & {\cal O}(0.1)    \\
0       & {\cal O}(0.1)    & {\cal O}(1.0)  
\end{pmatrix}
\label{mixing}
\end{eqnarray} 
In the above matrix, the non-zero elements come from 
$V_{\alpha i}^{\nu \nu}$, $V_{\beta j}^{N N}$,  $V_{\beta k}^{N S}$, 
$V_{\gamma j}^{S N}$ and $V_{\gamma k}^{S S}$ while all other terms are 
negligibly small. These non-zero mixings would contribute sizeably to 
the predicted neutrinoless double beta decay rate.
Thus, Eq. (\ref{nuNS}) can be rewritten for the 
fields of flavour neutrinos
$\nu_{\alpha L}$ and the heavy RH neutrinos $N^c_{\beta L}$  as:
\begin{eqnarray}
\nu_{\alpha L} &\cong &V_{\alpha i}^{\nu \nu} \nu_{iL} + V_{\alpha j}^{\nu N} N_{jL}\,,
\nonumber \\
 N_{\beta L}^c &=& V_{\beta i}^{ N \nu} \nu_{iL}+ V_{\beta j}^{N N} N_{jL} 
+ V_{\beta k}^{N S} S_{kL}\,. 
 \label{contibution}
\end{eqnarray}

 As we have indicated, in the considered model we have
$|V_{\alpha j}^{\nu N}|\sim 10^{-6}$. Correspondingly,
the contribution to the $0\nu\beta\beta$ decay amplitude 
arising from the coupling of $N_{jL}$ 
to the electron in the LH (i.e., V-A) charged lepton current 
involves the factor $(V_{e j}^{\nu N})^2$ and  
is negligible. 


\subsection{Gauge Boson Masses}
\label{sec:gbmass}

We briefly summarize here the gauge bosons masses and mixing in our model 
which will be used in estimating half-life of neutrinoless double beta decay 
process. Besides the SM gauge bosons $W_L^{\pm}$ and $Z$, there are 
right-handed gauge bosons $W_R^{\pm}$ and $Z^{\prime}$ which get their masses 
from left-right symmetry breaking. 
Following ref~\cite{Senjanovic:1978ev} and choosing VEVs of the Higgs fields as,
\begin{equation}
  \langle H^0_{R} \rangle = v_{R} \,,
\quad \langle \phi^0_{1,2} \rangle = v_{1,2} \,,
\end{equation}
the mass matrix for charged gauge bosons, 
in the basis $(W_L^+  \ W_R^+)$ can be written as,
\begin{equation*}
  \mathbb M_{CGB}
=
\left\lgroup
\begin{matrix} \frac{g_L^2 v^2}{2} & & - g_L g_R v_1 v_2 \\ - g_L g_R v_1 v_2 & & \frac{g_R^2}{2}(\frac{1}{2}v_R^2 + v^2) \end{matrix}
\right\rgroup \,,
\end{equation*}
where $v^2 = v_1^2 + v_2^2$ 
 and $g_R = g_L$. 
The physical mass for extra charged 
gauge boson is given by 
\begin{eqnarray}
M_{W_R} &\simeq & \frac{1}{2} g_R v_R \,. 
\end{eqnarray}
The mixing angle between $W_R$ and $W_L$ is defined as 
$$\tan 2 \theta_{LR} \approx 8 \frac{g_L}{g_R} \frac{v_1 v_2}{v_R^2}$$.
 
The neutral gauge boson mass matrix is given by 
\begin{equation}
\mathbb M_{NGB} =
\left\lgroup
 \begin{matrix} 
\frac{g_{L}^2 v^2}{2} && - \frac{g_L  g_{R}}{2} v^2 && 0 \\ 
- \frac{g_L g_{R}}{2} v^2 && \frac{g_R^2}{2} (\frac{1}{2}v_R^2 + v^2) && - \frac{g_R g_{BL}}{4} v_R^2  \\ 
0 && - \frac{g_R g_{BL}}{4} v_R^2 && \frac{g_{BL}^2 v^2_R}{4} \end{matrix}
\right\rgroup \,. \nonumber
\end{equation}
 As can be easily checked, this mass matrix has one zero eigenvalue 
corresponding to the photon $A_\mu$.
After few simplification, the mass eigenstates  $Z_\mu$, $Z^{'}_\mu$ 
and $A_\mu$ are related to the weak eigenstates 
$(W^0_{L\mu}, W^0_{R\mu}, Z_{BL\mu})$ in the following way,
\begin{eqnarray}
W_{L\mu}^0 &=& \cos \theta_W Z_{L\mu} +\sin \theta_W A_{\mu}\,, \nonumber \\
W_{R\mu}^0 &=& \cos \theta_R Z_{R\mu} -\sin \theta_W \sin \theta_R Z_{L\mu} 
+ \cos \theta_W \sin \theta_R A_{\mu}\,, \nonumber \\
Z_{BL\mu}^0 &=& -\sin \theta_R Z_{R\mu}- \sin \theta_W \cos \theta_R Z_{L\mu} + 
\cos \theta_W \cos \theta_R A_\mu\,, \nonumber 
\end{eqnarray}
where 
\begin{eqnarray}
Z_{L\mu} &\equiv& Z_\mu \cos\xi + Z^{'}_\mu \sin\xi\,,
\nonumber \\
 Z_{R\mu} &\equiv& -\, Z_\mu \sin\xi + Z^{'}_\mu \cos\xi\,.
\end{eqnarray}
Here, the mixing angles are defined as $\tan \theta_R = g_{BL}/g_R$, 
$\tan \theta_W=g_Y/g_L$ with $g_Y = g_{BL}g_R/\sqrt{g^2_{BL}+g_R^2}$,  
while the mixing angle between the $Z$ and the heavy $Z^{'}$ reads: 
\begin{equation}
\tan 2 \xi \approx \frac{v_1 v_2}{v^2_R} \ \frac{-4 g_R^2 \sqrt{g_L^2g_R^2 + g^2_{BL}(g_L^2+g_R^2)}}{(g_{BL}^2+g_R^2)^2}\,.
\end{equation}
The physical mass for extra neutral gauge boson $Z^\prime$ is given by:
\begin{eqnarray}
M^2_{Z^{'}} &\simeq & \frac{1}{2} \left(g^2_{BL} + g^2_R \right) 
   \bigg[v^2_R +\frac{g^2_R v^2}{g^2_R+g^2_{BL}} \bigg]
\end{eqnarray}
The value of $\tan 2 \xi$ has to be smaller than $10^{-3}$ in order to 
satisfy the electroweak precision constraints in the limit 
$v_R^2 \gg v_1^2 +v_2^2$. With $v^2=v^2_1+v^2_2\simeq (246$~GeV)$^2$
and $g_R\simeq g_L = 0.653$, we have 
$\tan\theta_W=\frac{g_{BL}}{\sqrt{g^2_{BL}+g^2_{R}}}$, 
which implies 
$\frac{g^2_{BL}}{g^2_{L/R}}=\frac{\sin^2\theta_W}{1-2\sin^2\theta_W}\approx0.43$, 
where $sin^2\theta_W=0.231$. 
 Using this result and $g_L = g_R$, we get for 
the angle describing the  $Z-Z^{'}$ mixing: 
$|\tan 2\xi| \cong 2.67 v_1v_2/v^2_R$. The upper limit 
$|\tan 2\xi| < 10^{-3}$ implies:
\begin{equation}
\dfrac{v_1\,v_2}{v^2_R} < 3.75\times 10^{-4}\,.
\label{eq:v1v2ov2R}
\end{equation}
This in turn leads to the following upper limit on the 
$W_L - W_R$ mixing angle $\theta_{LR}$:
\begin{equation}
\theta_{LR} \cong 4\,\dfrac{v_1\,v_2}{v^2_R} < 1.50\times 10^{-3}\,.
\label{eq:thLR}
\end{equation}

  The left-handed gauge boson masses are similar 
to those of the SM gauge bosons with 
$g_Y = \left(g_R g_{BL} \right)/\left( g^2_R+g^2_{BL}\right)^{1/2} $, 
while the masses of the extra heavy gauge bosons are related as follows,
\begin{eqnarray}
M_{W_R} &\simeq & \frac{1}{2} g_R v_R \,, \\
M_{Z^{'}} &\simeq & \frac{\sqrt{g^2_{BL} + g^2_R}}{g_R} M_{W_R} 
{\simeq} 1.2 \, M_{W_R}\,.
\end{eqnarray}
The current experimental bound on $M_{W_R} > 5$ TeV is 
obtained in high energy 
collider experiements at LHC~\cite{ATLAS:2018dcj,ATLAS:2019isd,CMS:2018agk},
while the  low energy precision measurements \cite{Li:2020wxi,Dekens:2021bro} 
imply a lower bound on the $Z^{'}$ mass, i.e. $M_{Z^{'}}>6$ TeV.

\section{Neutrinoless double beta decay}
\label{sec:0nubb}

 Neutrinoless double beta decay process can be induced by the exchange of 
light active Majorana neutrinos, which is usually referred to as  
``the standard mechanism'', or by some other lepton number violating 
``non-standard mechanism'' associated with BSM physics.
In this section, we discuss the
standard and the new physics contributions to $0\nu\beta\beta$ decay 
amplitude and rate that arise in our model 
due to the exchange of the light Majorana neutrinos $\nu_i$, 
heavy Majorana neutrinos $N_{1,2,3}$ 
and sterile Majorana neutrinos $S_{1,2,3}$.

The charged current (CC) interaction Lagrangians for leptons and quarks, 
relevant for our further discussion,  
are given by:
\begin{eqnarray}
\mathcal{L}^{\rm \ell}_{\rm CC}
 &=& \sum_{\alpha=e, \mu, \tau}
\bigg[\frac{g_L}{\sqrt{2}} \overline{\ell}_{\alpha L} \gamma_\mu {\nu}_{\alpha L} W^{\mu}_L 
      +\frac{g_R}{\sqrt{2}} \overline{\ell}_{\alpha R} \gamma_\mu {N}_{\alpha R} W^{\mu}_R \bigg] + \text{h.c.} 
      \nonumber \\
&=& \frac{g_L}{\sqrt{2}}
\overline{e}_{L} \gamma_\mu {\nu}_{e L} W^{\mu}_L  
   + \frac{g_R}{\sqrt{2}}\overline{e}_{R} \gamma_\mu {N}_{e R} W^{\mu}_R + \text{h.c.} + \cdots 
\label{eqn:ccint-flavor-lepton}
\\
\mathcal{L}^{\rm q}_{\rm CC} &=& \bigg[\frac{g_L}{\sqrt{2}}\, \overline{u}_{\,L}\, \gamma_\mu {d}_{\,L}\, W^{\mu}_L 
      +\frac{g_R}{\sqrt{2}}\, \overline{u}_{R}\, \gamma_\mu d_{R}\, W^{\mu}_R \bigg] + \text{h.c.} 
\label{eqn:ccint-flavor-quark}
\end{eqnarray}
Using Eq. (\ref{contibution}), $\mathcal{L}^{\rm \ell}_{\rm CC}$
be rewritten as:
\begin{eqnarray}
\mathcal{L}^{\rm \ell}_{\rm CC}&=& \frac{g_L}{\sqrt{2}}
\bigg[ \overline{e}_{L} \gamma_\mu 
        \{\mbox{V}^{\nu\nu}_{e i} \nu_i + 
                      \mbox{V}^{\nu N}_{e i} N_i \} 
              W^{\mu}_L \bigg] +\mbox{h.c.} \nonumber \\
    & &  + \frac{g_R}{\sqrt{2}}
\bigg[ \overline{e}_{R} \gamma_\mu 
\{\mbox{V}^{N \nu}_{e i}\nu_i +
      \mbox{V}^{N N}_{e i}N_i +
                      \mbox{V}^{N S}_{e i} S_i \} 
              W^{\mu}_R \bigg] + \mbox{h.c.}
\label{eqn:ccint-mass}
\end{eqnarray}

In the present model, the heavy neutrino masses are around $1-1000$~GeV, 
$\nu_L-N^c_L$ mixing  $|V_{\nu N}| \leq 10^{-6}$ and $\nu_L-S_L$ mixing 
is vanishing. Other contributions involving the light-heavy neutrino mixings 
and the $W_L-W_R$ mixing are negligible. 
The  lepton Lagrangian that is relevant for the dominant contributions 
to $0\nu\beta\beta$ decay rate is:
\begin{eqnarray}
\mathcal{L}^{\rm \ell}_{\rm CC}&=& \frac{g_L}{\sqrt{2}}\,
\bigg[ \overline{e}_{\,L}\, \gamma_\mu 
        \{\mbox{V}^{\nu\nu}_{e\, i}\, \nu_i  
                      \}\, 
              W^{\mu}_L \bigg] +\mbox{h.c.} \nonumber \\
    & &  + \frac{g_R}{\sqrt{2}}
\bigg[ \overline{e}_{R} \gamma_\mu 
\{ \mbox{V}^{N\, N}_{e j}N_j +
                      \mbox{V}^{N S}_{e k} S_k \} 
              W^{\mu}_R \bigg] + \mbox{h.c.}
\label{eqn:ccint-mass}
\end{eqnarray}
Thus, in the considered model, the dominant contributions to 
the $0\nu\beta\beta$ decay amplitude are given by:
\begin{itemize}
\item the standard mechanism due to the exchange of light neutrino $\nu_i$, 
mediated by left-handed gauge boson $W_L$, i.e. due to purely left-handed 
(LH) CC interaction;
\item new contributions due to the exchange of heavy neutrinos $N_{1,2,3}$ 
and sterile neutrinos $S_{1,2,3}$, mediated by right-handed gauge boson $W_R$, 
i.e., due to purely right-handed (RH) CC interaction. 
 The contribution due to exchange of virtual $S_{1,2,3}$ 
is possible due to the mixing between $N^c_L$ and $S_L$.
\end{itemize}

 The so-called $<\lambda>$- and $<\eta>$- mechanism contributions 
$0\nu\beta\beta$ decay amplitude arising  
from the product of LH and RH lepton currents 
\cite{Doi:1985dx}
are sub-dominant being strongly suppressed.
The $<\lambda>$-mechanism contribution involves the factor 
$|V^{N \nu}_{ei}|(M_{W_L}/M_{W_R})^2 < 2.6\times 10^{-10}$, 
where we have used $|V^{N \nu}_{ei}| = 10^{-6}$, 
$M_{W_L} = 80.38$ GeV and $M_{W_R} > 5$ TeV, 
while the $<\eta>$-mechanism contribution is suppressed by the factor 
$|V^{N \nu}_{ei} \sin\theta_{LR}| < 10^{-9}$. As a consequence, 
we neglect these contributions in the analysis which follows.


The Feynman diagrams for the dominant contributions of interest  
to the $0\nu\beta\beta$ decay amplitude are shown in Fig. \ref{FD}, 
where the first diagram from the left corresponds to the standard mechanism, while 
the second and third diagrams correspond to the new contributions mediated 
by $N_{1,2,3}$ and $S_{1,2,3}$, respectively.
\begin{figure*}[!ht]
       \includegraphics[width=1.0\textwidth]{Feyn_Diad_DSS.png}
\caption{Feynman diagrams for the process of neutrinoless double beta decay 
mediated by the exchane of virtual (a) light Majorana neutrinos $\nu_i$
(the standard mechanism), (b) heavy  neutrinos $N_R$ 
(heavy Majorana neutrinos $N_{1,2,3}$) and 
(c) heavy sterile neutrinos $S_L$ (heavy Majorana neutrinos $S_{1,2,3}$). 
}
\label{FD}
\end{figure*}

When $0\nu\beta \beta$ decay is mediated by only light Majorana neutrinos 
$\nu_i$, the inverse half-life for this process can be expressed as,
\begin{eqnarray}
 \left[T_{1/2}^{0\nu}\right]^{-1} &=& g^4_{\rm A}\,G^{0\nu}_{01}~|{\cal M}^{0\nu}_\nu\big|^2 ~ |\eta_\nu|^2   
\label{std_halflife}
\end{eqnarray}
where $g_{\rm A}$ is the axial coupling constant, 
 $G^{0\nu}_{01}$ is the phase space factor, 
${\cal M}^{0\nu}_\nu$ is the Nuclear Matrix Elements (NME) for 
light neutrino exchange and $\eta_\nu$ is a dimensionless particle physics 
parameter that is a measure of lepton number violation. 
Considering both the standard mechanism and the new contributions 
to this decay process in our model, the inverse half life can be written as:
\begin{eqnarray}
 \left[T_{1/2}^{0\nu}\right]^{-1} &=& g^4_{\rm A}\,G^{0\nu}_{01}\bigg[ |{\cal M}^{0\nu}_\nu \cdot \eta_\nu|^2 + 
 |{\cal M}^{0\nu}_N \cdot \left(\eta_{N}+ \eta_{S} \right) \big|^2 \bigg]
\label{eq:Hlife-a} ,
\end{eqnarray}
where ${\cal M}^{0\nu}_N$ is the Nuclear Matrix Elements (NME) for the 
heavy neutrino exchange and 
 $\eta_{N}$ and $\eta_{S}$ are lepton number violating parameters associated 
with the exchange of the heavy neutrinos $N_{1,2,3}$ and $S_{1,2,3}$.

In the analysis and the numerical estimates which follow, we will use 
a mildly quenched value of the axial coupling constant $g_{\rm A} = 1.00$, 
the unquenched value being, as is well known, $g_{\rm A} = 1.27$.
If it turns out that  $g_{\rm A}$ is actually not quenched, that will reduce 
the estimates of the $0\nu\beta\beta$ decay half-lives made in the present study by a factor of 2.60. 

The interference term of the light neutrino $\nu_{1,2,3}$ and 
the heavy neutrinos $N_{1,2,3}$ and $S_{1,2,3}$ exchange contributions
to the $0\nu\beta\beta$ decay amplitude,
which are generated respectively by LH ($V-A$) CC and RH ($V+A$) 
CC interactions, is strongly suppressed, being proportional to the 
electron mass~\cite{Halprin:1983ez} (see also~\cite{Faessler:2011qw}) 
and we have neglected it in Eq. (\ref{eq:Hlife-a}). 

 The values of $G^{0\nu}_{01}$ and the NMEs for both light and heavy 
neutrino exchange mechanism are distinct for different isotopes and 
can be found, e.g., in \cite{Ejiri:2019ezh}. 
 We present in Table \ref{tab:nucl-matrix} the values obtained 
by six different groups of authors using different methods of NME 
calculation. Of particular importance 
for the estimates of the relative magnitude of the new non-standard 
contributions in the $0\nu\beta\beta$ decay amplitude
with respect to the contribution of the standard mechanism 
is the ratio ${\cal M}^{0\nu}_N/{\cal M}^{0\nu}_\nu$.
As it follows from Table \ref{tab:nucl-matrix}, 
the ratio  ${\cal M}^{0\nu}_N/{\cal M}^{0\nu}_\nu$ 
predicted by each of the six cited groups using different 
methods of NME calculation   
is essentially the same for the four isotopes 
$^{76}$Ge,  $^{82}$Se,  $^{130}$Te, $^{136}$Xe -- 
it varies with the isotope by not more than 
$\sim 15\%$. At the same time, for a given isotope 
the ratio of interest obtained by the six different 
methods of NME calculation quoted in 
Table \ref{tab:nucl-matrix} varies by a factor 
of up to $\sim 3.5$. In view of this we will use 
the values of the NMEs for $^{76}$Ge as reference values 
in our numerical analysis. For the minimal and maximal values 
of the ratio  ${\cal M}^{0\nu}_N/{\cal M}^{0\nu}_\nu$ 
for $^{76}$Ge we get from Table \ref{tab:nucl-matrix}:
\begin{equation}
22.2 \lesssim \dfrac{{\cal M}^{0\nu}_N}{{\cal M}^{0\nu}_\nu} 
\lesssim 76.3\,,~~~{\rm  ^{76}Ge}\,. 
\label{eq:MnuMN}
\end{equation}
They correspond respectively to ${\cal M}^{0\nu}_\nu = 4.68$ and 5.26.

\begin{table}[h]
 \centering
\vspace{10pt}
 \begin{tabular}{|l|cc|cc|cc|cc|c|}
 \hline \hline
  &  $^{76}$Ge 
& & $^{82}$Se  & & $^{130}$Te & &$^{136}$Xe &  \\[2mm] 
{methods}  & {${\cal M}^{0\nu}_\nu$}  
 & {${\cal M}^{0\nu}_N$} & {${\cal M}^{0\nu}_\nu$} & {${\cal M}^{0\nu}_N$}&{${\cal M}^{0\nu}_\nu$} &{${\cal M}^{0\nu}_N$} &{${\cal M}^{0\nu}_\nu$}  & {${\cal M}^{0\nu}_N$} \\ 
 \hline 
 dQRPA \cite{Fang:2018tui} & $3.12$  & $187.3$ & $2.86$  & $175.9$& $2.90$& $191.4$ &$1.11$ & $66.9$\\ 
  \hline
  QRPA-Tu \cite{PhysRevC.87.045501,PhysRevD.90.096010}  & $5.16$  & $287.0$ &$4.64$  & $262.0$&$ 3.89$ & $264.0$ & $2.18 $ & $152.0$ \\ 
  \hline
  QRPA-Jy \cite{PhysRevC.91.024613}  & $5.26$  & $401.3$ & $3.73$ &$287.1$ &$4.00$ &$338.3$ & $2.91 $ & $186.3$ \\ 
  \hline
   IBM-2 \cite{PhysRevC.87.014315} & $4.68$  & $104$ & $3.73$ & $82.9$ & $3.70 $& $91.8 $ & $3.05 $ & $72.6 $ \\ 
  \hline
  CDFT \cite{PhysRevC.90.054309,PhysRevC.91.024316,PhysRevC.95.024305} & $ 6.04$ & $209.1 $  & $5.30 $ & $189.3 $ & $4.89 $ & $193.8 $ & $4.24 $ & $166.3 $ \\ 
\hline 
 ISM \cite{Menendez:2017fdf} & $2.89 $ & $130 $  & $2.73 $ & $121 $ & $2.76 $ & $146 $ & $ 2.28$ & $116 $ \\ 
\hline
\hline
$G^{0\nu}_{01}$ $[{10^{-14} \rm yrs}^{-1}] $ \cite{Horoi:2017gmj} & $0.22$ & & 1 & & $1.4$ & &$1.5$ &\\
\hline \hline
 \end{tabular}
 \caption{Values of Nuclear Matrix Elements for various isotopes calculated 
by different methods for light and heavy neutrino exchange. 
Here QRPA-Jy uses CD-Bonn short range correlations (SRC) and the rest 
use Argonne SRC, with minimally quenched $g_A = 1$. The last row shows the 
phase space factor $G^{0\nu}_{01}$ for various 
isotopes \cite{Horoi:2017gmj,Ejiri:2019ezh}.}
 \label{tab:nucl-matrix}
\end{table} 

The dimensionless particle physics parameters $\eta_\nu$, 
$\eta_N$ and $\eta_s$ in Eq. (\ref{eq:Hlife-a}) are functions of 
neutrino masses, 
mixing parameters and CPV phases and can be expressed as:
\begin{eqnarray}
\label{eq:eta_LL} 
|\mathcal{\eta}_{\nu}| = \sum_{i=1,2,3} \frac{{\mbox{V}^{\nu \nu}_{ei}}^2\, m_{\nu_i}}{m_e} \,, \\
 \label{eq:etaN}
 |\mathcal{\eta}_{N}| = m_p \left(\frac{M_{W_L}}{M_{W_R}}\right )^4 \sum_{j=1,2,3} \frac{{\mbox{V}^{N N}_{ej}}^2}{m_{N_j}} \,, \\                   
\label{eq:etaS}  
|\mathcal{\eta}_{S}| = m_p \left(\frac{M_{W_L}}{M_{W_R}}\right )^4 \sum_{k=1,2,3} \frac{{\mbox{V}^{N S}_{ek}}^2}{m_{S_k}} \,.
\end{eqnarray}
where $m_e$ and $m_p$ are the electron and proton masses. 
The quantity  
$m_e|\mathcal{\eta}_{\nu}| \equiv | m^\nu_{\beta\beta, L}|$ is the effective 
Majorana mass (EMM) associated with the standard mechanism of  
$0\nu\beta\beta$ decay (see, e.g., \cite{Bilenky:1987ty,Bilenky:2001rz}). 
The mixing parameters in Eq. (\ref{eq:eta_LL}) are given in 
Appendix \ref{app:lrdssm}. In the framework of our model we have: 
${V}^{\nu \nu} \approx U_\nu$, 
${V}^{N N}\approx U_N$ and 
${V}^{N S}\equiv M_{RS} M^{-1}_S U_S$.
We recall that 
 $U_N = i\,U^*_\nu$, 
$U_S = U_\nu$ and 
$U_\nu \equiv U_{PMNS}$ (Eqs. (\ref{eq:UN}) and (\ref{eq:US})).

 The expressions for  $|\mathcal{\eta}_{N}|$ and  $|\mathcal{\eta}_{S}|$
in Eq. (\ref{eq:eta_LL}) are obtained under the condition
$\langle p^2 \rangle \ll  M^2_i$, where $\sqrt{\langle p^2 \rangle}$ 
is the average momentum  exchanged in the process of $0\nu\beta\beta$ decay 
and $M_i$ here is a  generic notation for the masses of 
$N_{1,2,3}$ and $S_{1,2,3}$. The chiral structure of the matrix elements 
involving virtual  $N_{1,2,3}$ and $S_{1,2,3}$ propagators in the case of the 
heavy neutrino exchange contribution is given by:
\begin{eqnarray}
& &P_{R}\frac{\slashed{p}+M_i}{p^2-M_i^2}P_{R} = \frac{M_i}{p^2-M_i^2}P_{R}\,,  
\end{eqnarray}
where $P_{R} = (1 + \gamma_5)/2$ is the RH projection operator.
A typical value of the neutrino virtuality is 
$\langle p^2 \rangle \cong (190 \,\mbox{MeV})^2$ 
(see, e.g., \cite{Babic:2018ikc}). Thus, in the case of interest, we have 
$ p^2 \ll  M^2_i$ and the heavy state propagators reduce 
to a good approximation to $1/M_i$.

 It proves convenient for our further analysis 
to rewrite the inverse half-life in terms of one 
particle physics parameter -- generalised effective Majorana mass (GEMM) -- 
that contains the lepton number violating information in it:
\begin{eqnarray}
\left[T_{1/2}^{0\nu}\right]^{-1} &=& G^{0\nu}_{01}\bigg[ |{\cal M}^{0\nu}_\nu  \eta_\nu|^2 + 
{\cal M}^{0\nu}_N \big| \eta_{N} + \eta_{S}  \big|^2  \bigg] \nonumber \\
&=&  G^{0\nu}_{01} \bigg| \frac{{\cal M}^{0\nu}_\nu}{m_e} \bigg|^2   \bigg[\big|m^{\nu}_{\beta \beta, L} \big|^2 
+  \big|m^{N}_{\beta \beta, R} +  m^{S}_{\beta \beta, R} \big|^2 \bigg] \nonumber \\
&=&G^{0\nu}_{01} \bigg |\frac{{\cal M}^{0\nu}_\nu}{m_e}\bigg |^2  
|m^{\rm eff}_{\beta \beta, L,R}|^2\,,
\label{halflife:no-int}
\end{eqnarray}

where \cite{Babic:2018ikc}
\begin{eqnarray}
m^{N}_{\beta \beta ,R} &=& \sum_{j}  
m_p m_e \frac{\mathcal{M}_N^{0\nu}}{\mathcal{M}_\nu^{0\nu}}
\frac{M^4_{W_L}}{M^4_{W_R}}  \frac{{\mbox{V}^{N N}_{ej}}^2}{m_{N_{j}}}\,,
\label{eqn:mee_NR}   
\nonumber \\
 m^{S}_{\beta \beta ,R} &=& \sum_{k} 
m_p m_e \frac{\mathcal{M}_N^{0\nu}}{\mathcal{M}_\nu^{0\nu}}
\frac{M^4_{W_L}}{M^4_{W_R}}  \frac{{\mbox{V}^{N S}_{ek}}^2}{m_{S_{k}}}\,.
\label{eqn:mee_SR} 
\end{eqnarray}

 It follows from Eqs. (\ref{eq:etaN}), (\ref{eq:etaS}) 
and (\ref{eqn:mee_NR}) that 
the new physics contributions to the $0\nu\beta\beta$ decay amplitude
are suppressed, in particular, by the factor  
$M^4_{W_L}/M^4_{W_R} < (1.6\times 10^{-2})^4$,
where  $M_{W_L} = 80.38$~GeV is the SM $W$-bosons mass and  
we have used the lower bound $M_{W_R} > 5$ TeV  
\cite{ATLAS:2018dcj,ATLAS:2019isd,CMS:2018agk,Li:2020wxi,Dekens:2021bro}.
Fixing  $M_{W_R}$ at, e.g., = 5.5 TeV, we have for the ratio 
$\left( M^4_{W_L}/M^4_{W_R} \right)  \sim {\cal O}(10^{-8})$. 
Taking further the masses of $S_k$ and $N_j$ in the ranges 
respectively of $(10^2 - 10^4)$ GeV and $(1 - 10^2)$ GeV,
the mixing ${V}^{N N}_{ej} \approx U_N$ and 
${V}^{N S}_{ek}\equiv M_{RS} M^{-1}_S U_S$ from Appendix \ref{app:lrdssm}, 
one finds that the new physics contributions can be  
in the $0.01-0.1$~eV range (see Table \ref{tab:half-life}), 
i.e., within the experimental search sensitivity.


We note that, since the dominant contributions to 
$0 \nu\beta\beta$ decay arises from more than one contribution, 
it is also possible that there might be interference between them 
in the decay rate of the process. The interference of light 
neutrino ($\nu_i$) contribution due to purely $V-A$ interaction 
involving LH  currents 
with either of the heavy neutrino $N_j$ and $S_k$ contributions, 
which are generated by  $(V+A)$ interaction with RH currents,  
is suppressed, as we have indicated earlier. However, the interference 
between the contributions of the heavy neutrinos $N_j$ and $S_k$   
both involving RH currents, in general, can't be neglected.  
 In the case when this interference is not taken into consideration, 
the generalised effective Majorana mass 
is determined by the sum of individual 
contributions of the three types of neutrinos $\nu_i$, $N_j$ and $S_k$:
 \begin{equation}
|m^{\rm eff}_{\beta \beta, L,R}| \equiv m^{\nu+N+S}_{ee} = 
\left (\big|m^{\nu}_{\beta \beta, L} \big|^2 
 +  \big|m^{N}_{\beta \beta, R} \big|^2 + \big| m^{S}_{\beta \beta, R} \big|^2 
\right)^{\frac{1}{2}}\,.
\label{eq:no-int}
 \end{equation}
 Accounting for the interference, the 
generalised effective Majorana mass can be written as:
 \begin{eqnarray}
|m^{\rm eff}_{\beta \beta, L,R}| \equiv m^{\nu+|N+S|}_{ee} &=& 
\left( \big|m^{\nu}_{\beta \beta, L} \big|^2 
 +  \big|m^{N}_{\beta \beta, R} + m^{S}_{\beta \beta, R} \big|^2 \right )^{\frac{1}{2}}
  \nonumber \\
 &=& \left( (m^{\nu+N+S}_{ee})^2 +  
2 \text{Re} (m^{N}_{\beta \beta, R} \cdot m^{S^{*}}_{\beta \beta, R})
\right)^{\frac{1}{2}}\,.
 \label{eq:int}
 \end{eqnarray}   
 In order to assess the relevance of the interference term 
$2 \text{Re} (m^{N}_{\beta \beta, R} \cdot m^{S^{*}}_{\beta \beta, R})$
in our study, we consider both the cases  of neglecting it and of taking 
it into account.

 We express next the three terms 
in the generalised effective Majorana mass in terms of 
the PMNS mixing angles, Dirac and Majorana CPV phases present 
in the PMNS matrix, the three light neutrino masses and, 
in the case of the non-standard contributions, the masses of 
$N_{1,2,3}$ and of $S_{1,2,3}$.


The  effective Majorana mass term for standard mechanism can be written as 
(see, e.g., \cite{Bilenky:1987ty,Bilenky:2001rz}):
\begin{subequations}
\begin{eqnarray}
\left| m^\nu_{\beta\beta, L} \right| &=& \left|\sum_{i=1}^{3} U_{ei}^2m_i \right|
=  \left|m_1c_{12}^2c_{13}^2 + m_2s_{12}^2c_{13}^2e^{i\alpha} + m_3s_{13}^2e^{i(\beta-2\delta)} \right|
 \end{eqnarray}
  \label{effmassnu}
\end{subequations}
 where $m_1,m_2,m_3$ are masses of the light Majorana neutrinos 
$\nu_{1,2,3}$ and we have used the standard parametrization of the PMNS matrix.
 Defining 
\begin{equation}
C_N =  m_e\,m_p\,\frac{M^{0\nu}_N}{M^{0\nu}_\nu} \frac{ M^4_{W_L}}{M^4_{W_R}}\,, 
\label{eq:CN}
\end{equation}
the expression for $m_{\beta \beta, R}^N$ can be cast in the form:
\begin{eqnarray}
\left| m^N_{\beta\beta,R} \right| &=&\frac{C_N}{{m_{N_1}}} \left| \bigg[U_{e1}^2 +
\frac{U_{e2}^2 e^{i\, \alpha}{m_{N_1}}}{m_{N_2}} 
+\frac{U_{e3}^2 e^{i\, \beta}{m_{N_1}}}{m_{N_3}} 
\bigg] \right| \nonumber \\
&=& \frac{C_N}{{m_{N_1}}} \left| \bigg[U_{e1}^2 +
\frac{U_{e2}^2 e^{i\, \alpha}{m_2}}{m_1} 
+\frac{U_{e3}^2 e^{i\, \beta}{m_3}}{m_1} 
\bigg] \right|  \nonumber \\
&=& \frac{C_N}{{m_{N_1}m_1}} \left| m^\nu_{\beta\beta,L} \right|\,, 
\hspace{1cm}\text{NO case}\,,
\label{eq:M1} 
\end{eqnarray}
\begin{eqnarray}
\left| m^N_{\beta\beta,R} \right|&=&\frac{C_N}{{m_{N_3}}}\left| \bigg[\frac{U_{e1}^2 {m_{N_3}}}{m_{N_1}} +
\frac{U_{e2}^2 e^{i\, \alpha}{m_{N_3}}}{m_{N_2}} 
+U_{e3}^2 e^{i\, \beta} \bigg]\right|  \nonumber \\
&=& \frac{C_N}{{m_{N_3}}}\left| \bigg[\frac{U_{e1}^2 m_{1}}{m_{3}} +
  \frac{U_{e2}^2 e^{i\, \alpha} m_{2}}{m_{3}}
 +U_{e3}^2 e^{i\, \beta} \bigg]\right|  \nonumber \\
&=& \frac{C_N}{{m_{N_3}m_3}} \left| m^\nu_{\beta\beta,L} \right|\,, 
\hspace{1cm}\text{IO case}\,, 
\label{eq:M2} 
\end{eqnarray}
 where we have used Eqs. (\ref{eqn:Mne1}) and (\ref{eqn:Mne2}).
We see that in the considered setting 
the contribution due to exchange of the heavy Majorana neutrinos 
$N_{1,2,3}$ is proportional to the standard contribution due to the 
light Majorana neutrino exchange, 
$|m^N_{\beta\beta,R}| \propto |m^\nu_{\beta\beta,L}|$.


 We consider next $m^S_{\beta\beta, R}$. 
It follows from 
Eq. (\ref{eq:massrel}) that $m_{N_i}= \frac{k^2_{rs}}{m_{S_i}}$.
As it is described in Appendix (\ref{app:lrdssm}), 
the mixing ${V}^{N S}_{ek}\equiv M_{RS} M^{-1}_S U_S$. 
We can diagonalize $M_S$ as 
$M_S = {U_S}M^{D}_S {U_S}^T$. 
Since $M_{RS}=k_{rs}I$, the mixing 
\begin{eqnarray}
&&{V}^{N S}_{ek} = k_{rs} ({U_S}M^{D}_S {U_S}^T)^{-1} U_S\, \nonumber \\
&&\hspace*{0.95cm}=k_{rs}{U^{*}_S} \,\text {diag}(1/m_{S_1},1/m_{S_2},1/m_{S_3})
\label{massmixing}
\end{eqnarray}
 Using the relation 
$m_{N_i}= \frac{k^2_{rs}}{m_{S_i}}$
and Eqs. (\ref{eqn:Mne1}) - (\ref{eqn:MSe2}), 
the expression for $m^S_{\beta \beta, R}$ can be written as:
\begin{eqnarray}
\left| m^S_{\beta\beta, R} \right|&=& C_N\, k^2_{rs}\left| \bigg[\frac{U_{e1}^2}{m^3_{S_1}} +
\frac{U_{e2}^2 e^{i\, \alpha}}{m^3_{S_2}} 
+\frac{U_{e3}^2 e^{i\, \beta}}{m^3_{S_3}} 
\bigg]\right|   \nonumber \\
&=& \left| C_N \bigg[\frac{U_{e1}^2 \,  m_{N_1}}{m^2_{S_1}} +
\frac{U_{e2}^2 e^{i\, \alpha} \,  m_{N_2}}{m^2_{S_2}} 
+\frac{U_{e3}^2 e^{i\, \beta} \,  m_{N_3}}{m^2_{S_3}} 
\bigg]\right|   \nonumber \\
&=& \left| \frac{C_N \,  m_{N_1} \,m_1 \,m^2_3}{m^2_{S_3}\,m^3_1} 
\bigg[U_{e1}^2 +
U_{e2}^2 e^{i\, \alpha}\,\frac{m^3_1}{m^3_{2}}
+ U_{e3}^2 e^{i\, \beta}\, \frac{m^3_1}{m^3_3}
\bigg]\right|\,~~~\text{ NO case}\,, 
\label{eq:M3} 
\end{eqnarray}
\begin{eqnarray}
\left| m^S_{\beta\beta,R} \right| &=& \left|C_N \, m_{N_3} \bigg[\frac{U_{e1}^2}{m^2_{S_1}} \frac{m_{3}}{m_{1}} +
\frac{U_{e2}^2 e^{i\, \alpha}}{m^2_{S_2}}\frac{m_{3}}{m_{2}} 
+\frac{U_{e3}^2 e^{i\, \beta}}{m^2_{S_3}} 
\bigg]\right|   
\nonumber \\
&=& \left|\frac{C_N \, m_{N_3}\,m_3\,m^2_2}{m^2_{S_2}\,m^3_3} 
 \bigg[U_{e1}^2\,\frac{m^3_3}{m^3_1} +
 U_{e2}^2 e^{i\, \alpha}\, \frac{m^3_3}{m^3_{2}}
 + U_{e3}^2 e^{i\, \beta}
 \bigg]\right|\,,~~~\text{ IO case}\,. 
\label{eq:M4} 
\end{eqnarray}

 It follows from Eqs. (\ref{eq:M1}) - (\ref{eq:M4})
that $|m^N_{\beta\beta,R}|$ and $|m^S_{\beta\beta,R}|$
exhibit very unusual dependence on the lightest neutrino mass 
$m_{1(3)}$:  $|m^N_{\beta\beta,R}|\propto 1/m_{1(3)}$ 
and  $|m^S_{\beta\beta,R}|\propto 1/m^2_{1(3)}$.
Correspondingly, the new physics contributions to the 
$0 \nu\beta\beta$ decay amplitude are strongly enhanced at relatively 
small values of $m_{1(3)}$. We will discuss this 
dependence in greater detail in the next Section. 
We will show, in particular, that in the considered scenario with  
$m_{S_k} \sim (10^2 - 10^4)$ GeV and $m_{N_j}\sim (1-100)$ GeV of interest,  
the lightest neutrino mass $m_{1(3)}$  cannot be smaller than 
$\sim 10^{-4}$ eV. We will also show that due to 
the indicated enhancement the new contributions   
dominate over the standard mechanism contribution for 
$m_{1(3)}\sim (10^{-4} - 10^{-2})$ eV.


\begin{table}[h!]
\begin{center}
	\begin{tabular}{cc|c|c}
	\hline
	Isotope & $T_{1/2}^{0 \nu}~\text{ yrs}$ & $m_{\beta \beta}^{0 \nu}~[\text{eV}]$ & Collaboration \\
	\hline
	$^{76}$Ge	& $> 1.8\times10^{26}$	& $< (0.08 - 0.18)$ 	& GERDA~\cite{GERDA:2020xhi} 			\\
	$^{76}$Ge	& $> 2.7\times10^{25}$	& $< (0.2 - 0.433)$ 	& MAJORANA DEMONSTRATOR~\cite{Majorana:2019nbd} 			\\
		& $> 8.3\times10^{25}$	& $< (0.113 - 0.269)$ 	& ~\cite{Majorana:2022udl} 			\\
	$^{82}$Se	& $> 3.5\times10^{24}$	& $< (0.311 - 0.638)$ 	& CUPID-0~\cite{PhysRevLett.123.032501} 			\\
	$^{130}$Te	& $> 2.2\times10^{25}$	& $< (0.09 - 0.305)$ 	& CUORE~\cite{CUORE:2021mvw} 			\\
	$^{136}$Xe	& $> 3.5 \times10^{25}$	& $< (0.093 - 0.286)$ 	& EXO~\cite{EXO-200:2019rkq} 				\\
	$^{136}$Xe	& $> 1.07\times10^{26}$	& $< (0.061 - 0.165)$ 				& KamLAND-Zen~\cite{KamLAND-Zen:2016pfg} 		\\
		& $> 2.3\times10^{26}$	& $< (0.036 - 0.156)$ 				& ~\cite{KamLAND-Zen:2022tow} 		\\
	\\
	\hline
	\label{tab:half-life}
	\end{tabular}	
\caption{The current lower limits on the half life $T_{1/2}^{0 \nu}$ and upper 
limits on the effective mass parameter $m_{\beta\beta}^{0 \nu}$ of neutrinoless 
double beta decay for different isotopes. The range for the 
effective Majorana mass parameter comes from uncertainties in the 
nuclear matrix element.}
\end{center}
\end{table}

\section{Phenomenological analysis}
\label{sec:PA}

 In the present Section, we will discuss the effects of the new physics 
contributions to the $0 \nu\beta\beta$ decay amplitude 
on the predictions for the effective Majorana mass 
and the $0 \nu\beta\beta$ decay half-life. 
We recall that if $0 \nu\beta\beta$ decay will be observed, 
the data on the  half-life of $0 \nu\beta\beta$ 
decay generated by the standard mechanism can provide
important  information 
on the absolute scale of light neutrino masses and 
on the neutrino mass ordering \cite{Pascoli:2002xq,Pascoli:2001by}.
With additional input data about the values of the lightest neutrino mass 
$m_{1(3)}$ (or the sum of the neutrino masses), it might be possible to get 
information about the values of the Majorana phases in the PMNS matrix 
as well \cite{Bilenky:2001rz,Pascoli:2005zb}.
In what follows, we will investigate, in particular, how the quoted  
results are possibly modified by the new contributions to the 
 $0 \nu\beta\beta$ decay amplitude.


\subsection{Mass parameter ranges}
\noindent

The new non-standard contributions to the 
 $0 \nu\beta\beta$ decay amplitude,
$|m^N_{\beta\beta,R}|$ and $|m^S_{\beta\beta,R}|$, 
as it follows from Eqs. (\ref{eq:M1}) - (\ref{eq:M4}), 
have very peculiar dependence on the lightest neutrino 
mass $m_{1(3)}$. They are strongly enhanced and, 
as we are going to show below, are considerably larger 
that the standard mechanism contribution $|m^\nu_{\beta\beta,L}|$
at $m_{1(3)} \ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$<$}\  10^{-3}$ eV, where 
$|m^N_{\beta\beta,R}| > |m^\nu_{\beta\beta,L}|$,
$|m^S_{\beta\beta,R}| >> |m^\nu_{\beta\beta,L}|$ 
and $|m^S_{\beta\beta,R}| >> |m^N_{\beta\beta,R}|$.
At $m_{1(3)} \ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  5\times 10^{-2}$ eV, however, we have 
$|m^N_{\beta\beta,R}|,|m^S_{\beta\beta,R}| \ll |m^\nu_{\beta\beta,L}|$.
This implies that the most stringent conservative 
experimental upper limit on  $|m^{0\nu}_{\beta\beta}| < 0.156$ eV 
reported by the KamLAND-Zen collaboration [26] 
(see Table 4) applies to $|m^\nu_{\beta\beta,L}|$
since it corresponds to light neutrino masses 
$m_{1,2,3}\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  0.1$ eV. Actually, it follows from 
the quoted upper limit that 
\cite{ParticleDataGroup:2018ovx,Penedo:2018kpc}
$m_{1,2,3} \ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$<$}\  0.156/(\cos 2\theta_{12} - \sin^2\theta_{13}) \cong 0.55$ eV,
where we have used the $3\sigma$ allowed ranges of 
$\sin^2\theta_{12}$ and $\sin^2\theta_{13}$
given in Table \ref{tab2_oscipara} neglecting the minor 
differences in the ranges corresponding to NO and IO neutrino 
mass spectra. Thus, the largest light neutrino mass 
$m_{3(2)}$ is allowed to vary approximately between 
$\sqrt{\Delta m^2_{31(23)}}\cong 5\times 10^{-2}$ eV and 
0.55 eV.

The Cosmic Microwave Background (CMB)
data of WMAP and PLANCK experiments, combined with
supernovae and other cosmological and astrophysical data 
can be used to obtain information in the form of an upper 
limit on the sum of neutrinos masses and thus on  $m_{3(2)}$ 
(see e.g., ref.~\cite{Vagnozzi:2017ovm}). 
Depending on the model complexity and the input data used 
one typically finds
\cite{ParticleDataGroup:2020ssz} (see also \cite{Capozzi:2017ipn}): 
$\sum_j m_j <  (0.12 - 0.54)$ eV (95\% CL). 
The quoted conservative upper limit on $\sum_j m_j$ 
implies  $m_{3(2)}\lesssim 0.18$ eV. In our phenomenological and numerical 
analysis, we will use somewhat larger values of  $m_{3(2)}$, 
keeping in mind the existence of more stringent limits.
We recall further that in the model considered by us 
$m_{N_i} = k^2_d/m_i$, $m_{S_i} = (k^2_{rs}/k^2_d)m_i$, 
where   $k_d$ and $k_{rs}$ are real constant parameters.
Correspondingly, in the case of NO light neutrino mass spectrum, 
$m_1 < m_2 < m_3$, we have $m_{N_3}< m_{N_2} < m_{N_1}$ and 
$m_{S_1} < m_{S_2} < m_{S_3}$. For IO spectrum, $m_3 < m_1 < m_2$,  
we have instead: $m_{N_2}< m_{N_1} < m_{N_3}$ and 
$m_{S_3} < m_{S_1} < m_{S_2}$. In the double seesaw model 
under discussion, we should always have in the NO (IO) case   
$m_{N_{1(2)}} \ll m_{S_{3(2)}}$, i.e., $m_{S_{3(2)}} \ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  10\,m_{N_{1(2)}}$.

In what follows, we will consider the values of 
$m_{S_{3(2)}}$ and $m_{N_{1(3)}}$
in the intervals $(1 - 10)$ TeV  and $(10^2 - 10^3)$ GeV, respectively, 
while the mass of the lightest RH Majorana neutrino $N_{3(2)}$ 
will be assumed to satisfy $m_{N_{3(2)}} \geq 1 $ GeV.
The minimal value of $m_{N_{3(2)}}$ of 1 GeV should correspond to the 
maximal allowed value of $m_{3(2)} \cong 0.55$ eV considered by us.
As a consequence, we have: 
$k^2_d = {\rm min}(m_{N_{3(2)}})\,{\rm max}(m_{3(2)}) = 0.55$ eV GeV.  
We get similar value of $k^2_d$ if we use $m_{N_{3(2)}} = 10$ GeV 
and $m_{3(2)} \cong 0.05$ eV.
Given the value of $k^2_d$, the requirement that the mass 
of the heaviest RH Majorana neutrino 
$m_{N_{1(3)}}$ should not exceed $10^3$ GeV implies a lower limit on the mass 
of the lightest Majorana neutrino $m_{1(3)}$: 
$m_{1(3)} = k^2_d/m_{N_{1(3)}} \ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  0.55\times 10^{-3}$ eV.
Thus, for consistency with the chosen ranges of value of 
the heavy Majorana fermions in the model, the value of the lightest 
neutrino mass should not be smaller than about $5.5\times 10^{-4}$ eV.
In the numerical analysis, we will perform 
we will exploit the range $m_{1(3)} = (10^{-4} - 1.0)$ eV.
 
In the analysis which follows, we will use the values of the 
neutrino oscillation parameters given in Table \ref{tab2_oscipara}.
We set the Dirac phase $\delta = 0$. 
The Majorana phases $\alpha$ and $\beta$ are varied in the interval 
$[0,\pi]$. For the parameters  $M_{W_R}$, $m_{N_{1(3)}}$ $m_{S_{3(2)}}$ 
and the ratio $M^{0\nu}_N/M^{0\nu}_\nu$ the following reference values 
will be utilised: $M_{W_R} = 5.5$ TeV, $m_{N_{1(3)}} = 300$ GeV, 
$m_{S_{3(2)}} = 3$ TeV and $M^{0\nu}_N/M^{0\nu}_\nu \cong 22.2 - 76.3$ 
(concerning $M^{0\nu}_N/M^{0\nu}_\nu$, see Eq. (\ref{eq:MnuMN}) and 
the discussion related to it).


\subsection{Light neutrino contribution}
\label{sec:nuexch}

The phenomenology of the light neutrino contribution to 
the $0 \nu\beta\beta$ decay half-life, including the properties 
of the corresponding effective Majorana mass $| m^\nu_{\beta\beta,L}|$
have been extensively studied and are well known 
(see, e.g., \cite{ParticleDataGroup:2018ovx}). 
In this subsection, we summarise the main features of $| m^\nu_{\beta\beta,L}|$.

{\bf Normal Ordering}

In this case $| m^\nu_{\beta\beta,L}|$ (see 
Eq. \ref{effmassnu}) can be rewritten in terms of neutrino 
mass square differences as,
\begin{equation}
\left| m^\nu_{\beta\beta,L} \right|=
 \left| m_1 c_{12}^2c_{13}^2 +
 \sqrt{m^2_1 + \Delta m^2_{21}} s_{12}^2c_{13}^2e^{i\alpha} +
\sqrt{m^2_1 + \Delta m^2_{31}}s_{13}^2e^{i\beta} \right|\,.
\label{eq:mnuNO}
\end{equation}
The best fit values and the 3$\sigma$ allowed ranges of  
$s_{12}^2 \equiv \sin^2\theta_{12}$, $s_{13}^2 \equiv \sin^2\theta_{13}$ 
and of  $\Delta m^2_{21}$ and $\Delta m^2_{31}$ are given 
in Table \ref{tab2_oscipara}.

The case of hierarchical 
light neutrino mass spectrum corresponds to $m_1\ll m_2 < m_3$. 
In this case 
$m_2 \approx \sqrt{\Delta m^2_{21}} \approx 8.57\times 10^{-3}$ eV and
$m_3 \approx \sqrt{\Delta m^2_{31}} \approx 4.98\times 10^{-2}$ eV, 
and thus  $m_1 \lesssim 8\times 10^{-4}$ eV. 
Depending on the values of the Majorana phases, 
$| m^\nu_{\beta\beta,L}|$ can take values in the interval
$(0.4 - 4.8)\times 10^{-3}$ eV, where we have used the 
the $3\sigma$ allowed ranges of the relevant oscillation parameters.
At $m_1 = 10^{-4}$ eV we have:
$0.91\times 10^{-3}~{\rm eV} \lesssim | m^\nu_{\beta\beta,L}| 
\lesssim 4.37\times 10^{-3}$ eV. 

The effective Majorana mass $| m^\nu_{\beta\beta,L}|$ exhibits strong dependence 
on the values of the Majorana phases $\alpha$ and $\beta$
in the case of NO neutrino mass spectrum with partial hierarchy 
corresponding to $m_1 = (10^{-3} - 10^{-2})$ eV.
Indeed, for $\alpha = \pi$ and $\beta = 0$,
$| m^\nu_{\beta\beta,L}|$
is strongly suppressed for values of $m_1$ lying 
in the interval $(1.3 - 9.0)\times 10^{-3}$ eV, 
where  $| m^\nu_{\beta\beta,L}| \lesssim 2\times 10^{-4} $ eV 
due to cancellations (partial or complete) 
between the three terms in the expression of 
$| m^\nu_{\beta\beta,L}|$. 
Using the best fit values of the neutrino oscillation parameters, 
we find that a complete cancellation takes place 
and  $| m^\nu_{\beta\beta,L}| = 0$ at $m_1 \cong 2.26\times 10^{-3}$ eV. 
At the same time, at $m_1 = 2.26~(9.0)\times 10^{-3}$ eV, 
for example,  $| m^\nu_{\beta\beta,L}|\approx 5~(10)\times 10^{-3}$ eV
if  $\alpha = 0$ and $\beta = 0$.

As $m_1$ increases beyond $10^{-2}$ eV, 
$| m^\nu_{\beta\beta,L}|$ increases almost linearly with $m_1$ 
and at $m_1\cong 0.1$ eV enters the quasi-degenerate (QD)
neutrino mass spectrum region where 
$| m^\nu_{\beta\beta,L}| \ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  0.05$ eV.


{\bf Inverted Ordering}

In this case we have:
\begin{equation}
\left| m^\nu_{\beta\beta,L} \right|=
 \left|\sqrt{m^2_3 + \Delta m^2_{23} - \Delta m^2_{21}}\, c_{12}^2\,c_{13}^2 +
 \sqrt{m^2_3 + \Delta m^2_{23}}\, s_{12}^2\,c_{13}^2\,e^{i\alpha} +
m_3\,s_{13}^2\,e^{i\beta}) \right|\,.
\label{eq:mnuIO}
\end{equation}
The behavior of $| m^\nu_{\beta\beta,L}|$ as a function of 
the lightest neutrino mass $m_3$ is very different 
from the behavior in the NO case.
Given the fact that 
$m_2 = \sqrt{m^2_3 + \Delta m^2_{23}}\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  5\times 10^{-2}$ eV, 
$m_1 = \sqrt{m^2_3 + \Delta m^2_{23} - \Delta m^2_{21}}\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  5\times 10^{-2}$ eV, 
$s^2_{13} \cong  0.022$ and at $3\sigma$ we have 
$(c^2_{12} - s^2_{12}) \ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  0.31$, complete 
cancellation between the three terms in Eq. (\ref{eq:mnuIO})
is not possible. Actually, at $m^2_3 \ll \Delta m^2_{23}$, or 
equivalently,  
at $m_3 \lesssim 1.6 \times 10^{-2}$ eV, 
$| m^\nu_{\beta\beta,L}|$ practically does not depend on $m_3$. 
At these values of $m_3$ we have 
$\sqrt{\Delta m^2_{23}} \cos2\theta_{12} \lesssim | m^\nu_{\beta\beta,L}| 
\lesssim \sqrt{\Delta m^2_{23}}$. 
Using the $3\sigma$ allowed ranges of $\sqrt{\Delta m^2_{23}}$
and $\cos2\theta_{12}$ from Table \ref{tab2_oscipara} we get:
$1.51\times 10^{-2}~{\rm eV} \lesssim | m^\nu_{\beta\beta,L}| 
\lesssim 5.06 \times 10^{-2}~{\rm eV}$. 

If the  $0 \nu\beta\beta$ decay were generated by the standard mechanism only, 
the fact that in the case of hierarchical light neutrino mass spectrum 
the minimal value of $| m^\nu_{\beta\beta,L}|$ for the IH spectrum 
is approximately by a factor of 3.4 larger than the maximal 
value of  $| m^\nu_{\beta\beta,L}|$ for the NH spectrum 
opens up the possibility of obtaining information about the type of 
neutrino mass spectrum from a measurement of $| m^\nu_{\beta\beta,L}|$ 
\cite{Pascoli:2002xq}.

As $m_3$ increases beyond $1.6 \times 10^{-2}$ eV, 
$|m^\nu_{\beta\beta,L}|$ also increases and at 
$m_3 \cong 0.1$ eV enters the QD region where
$|m^\nu_{\beta\beta,L}|\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  0.03$ eV
growing linearly with $m_3$.

\subsection{The contribution due to the exchange of $N_{1.2.3}$}
\label{sec:Nexch}

The contribution due to the exchange of virtual $N_{1,2,3}$ 
in the NO and IO cases 
are given respectively in Eqs. (\ref{eq:M1}) and (\ref{eq:M2}) 
can be cast in the form:
\begin{eqnarray}
\left| m^N_{\beta\beta,R} \right| 
&=& \frac{C_N}{{m_{N_{1(3)}}m_{1(3)}}} \left| m^\nu_{\beta\beta,L} \right|\,, 
\hspace{1cm}\text{NO (IO)~case}\,,
\label{eqn:mnNOIO} 
\end{eqnarray}
where $C_N$ is defined in Eq. (\ref{eq:CN}) 
and $|m^\nu_{\beta\beta,L}|$ is the effective Majorana mass 
associated with the standard mechanism discussed in the 
preceding sub-section. Using the $M_{W_L} = 80.38$ GeV, 
the reference values of 
$M_{W_R} = 5.5$ TeV 
we get:
\begin{equation} 
\frac{C_N}{m_{N_{1(3)}}\,m_{1(3)}} \cong 
0.729\, (\frac{m_{N_{1(3)}}}{300~{\rm GeV}})^{-1}\,
 (\frac{m_{1(3)}}{10^{-4}~{\rm eV}})^{-1}\, 
\frac{M^{0\nu}_N}{M^{0\nu}_\nu}\,.
\label{eq:CNmnmi}
\end{equation}
Taking into account the minimal and maximal reference values of 
$M^{0\nu}_N/M^{0\nu}$ in Eq. (\ref{eq:MnuMN}) and 
fixing $m_{N_{1(3)}}$ to the reference value of 300 GeV,  
we find the possible range of values of the factor  
$C_N/(m_{N_{1(3)}}m_{1(3)})$ for four values of the lightest 
neutrino mass $m_{1(3)}$:\\
$C_N/(m_{N_{1(3)}}m_{1(3)}) \cong (16.2 - 55.6)$ for $m_{1(3)} = 10^{-4}$ eV,\\ 
$C_N/(m_{N_{1(3)}}m_{1(3)}) \cong (1.62 - 5.56)$ for $m_{1(3)} = 10^{-3}$ eV,\\ 
$C_N/(m_{N_{1(3)}}m_{1(3)}) \cong (0.32 - 1.11)$ for 
$m_{1(3)} = 5\times 10^{-3}$ eV,\\ 
$C_N/(m_{N_{1(3)}}m_{1(3)}) \cong (0.16 - 0.56)$ for $m_{1(3)} = 10^{-2}$ eV.\\ 

It is clear from these estimates that
for  $ 10^{-4}{\rm eV} \leq m_{1(3)} \leq 10^{-3}$ eV, 
the  contribution due to exchange of virtual $N_{1,2,3}$ is larger than 
the standard mechanism contribution: 
$| m^N_{\beta\beta,R}| > | m^\nu_{\beta\beta,L}|$. 
For $m_{1(3)}\sim (10^{-4} - 5\times 10^{-4})$ eV we have actually:
$| m^N_{\beta\beta,R}| \gg | m^\nu_{\beta\beta,L}|$. 
In this interval of values of $m_{1}$ 
in the NO case, $| m^N_{\beta\beta,R}|$ lies in the region 
corresponding to the IO neutrino mass spectrum 
if only the standard mechanism 
(i.e., only light Majorana neutrino exchange)
were operative in $0 \nu\beta\beta$ decay.   
The predicted values of 
$| m^N_{\beta\beta,R}|$ in the IO case 
are larger than the experimental limits on 
effective Majorana mass reported by the GERDA and KamLAND-Zen 
experiments (see Table \ref{tab:half-life}).

In the region  $m_{1(3)}\sim (10^{-3} - 10^{-2})$ eV 
we have roughly $| m^N_{\beta\beta,R}| \sim | m^\nu_{\beta\beta,L}|$ 
(see below), with $| m^N_{\beta\beta,R}|$ decreasing as $1/m_{1(3)}$. 
In the NO case, $| m^\nu_{\beta\beta,L}|$ can be strongly suppressed,
i.e., depending on the values of the Majorana phases 
it can have value $| m^\nu_{\beta\beta,L}| \leq 10^{-4}$ eV,
and in this case $| m^N_{\beta\beta,R}|$ will also be 
suppressed. At $m_{1} = 10^{-3}$ eV though at which 
$| m^\nu_{\beta\beta,L}| \cong 3\times 10^{-4}$ eV, $| m^N_{\beta\beta,R}|$ 
can be somewhat larger than $| m^\nu_{\beta\beta,L}|$  
owing to the relevant NME element ratio 
and can have a value $| m^N_{\beta\beta,R}| \cong 1.5\times 10^{-3}$ eV.
In the IO case, $| m^N_{\beta\beta,R}| \ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  | m^\nu_{\beta\beta,L}|$ 
in the discussed region. It can be larger than $| m^\nu_{\beta\beta,L}|$ 
by a factor of 2.

At  $m_{1(3)} > 10^{-2}$ eV, $| m^\nu_{\beta\beta,L}| \ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  | m^N_{\beta\beta,R}|$
and at  $m_{1(3)} \geq 5\times 10^{-2}$ eV, we have 
 $| m^\nu_{\beta\beta,L}| >> | m^N_{\beta\beta,R}|$ and 
the contribution due to the exchange of $N_{1.2.3}$ is 
subleading and practically negligible in both NO and IO cases.

For values of $m_{N_{1(3)}}$ smaller (larger)
than the considered 300 GeV, 
$|m^N_{\beta\beta,R}|$ will have  values which are larger (smaller) 
than those discussed above by the factor
$300~{\rm GeV}/m_{N_{1(3)}}$. Since in the considered scenario 
the mass of the lightest $N_{j}$ 
is assumed to satisfy $m_{N_{3(2)}} \geq 1$ GeV and 
is given by $m_{N_{3(2)}} = (m_{1(3)}/m_{3(2)}) m_{N_{1(3)}}$, 
$m_{1(3)}\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  5.5\times 10^{-4}$ eV, 
and from the data it follows that $m_{3(2)} \ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  5\times 10^{-2}$ eV,
for consistency one should have also $m_{N_{1(3)}}\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  100$ GeV. 

\subsection{The contribution due to the exchange of $S_{1.2.3}$}
\label{sec:Sexch}

The  important parameter 
for the contribution due to the exchange of $S_{1.2.3}$ 
in the NO (IO) case is the dimensionful factor
\begin{equation} 
C^{\rm NO(IO)}_{S} \equiv \frac{C_N\,m_{N_{1(3)}}\,m_{1(3)}\,m^2_{3(2)}} 
{m^2_{S_{3(2)}}\, m^3_{1(3)}}\,,~\hspace{1cm}\text{NO (IO)}\,.
\label{eq:CNS}
\end{equation}
Taking into account Eq. (\ref{eq:CNmnmi}), 
$C^{\rm NO(IO)}_{S}$ can be cast in the form:
\begin{equation} 
C^{\rm NO(IO)}_{S} = 0.729\times 10^{-6}~{\rm eV}\, 
\frac{m_{N_{1(3)}}}{300~{\rm GeV}}\, 
\left (\frac{m_{S_{3(2)}}}{3~{\rm TeV}}\right)^{-2}\,
\left(1 + \frac{\Delta m^2_{31(23)}}{m^2_{1(3)}}\right)\, 
\frac{{\rm M^{0\nu}_N}}{{\rm M^{0\nu}_\nu}}\,,~~\text{NO (IO)}\,.
\label{eq:CNS2}
\end{equation}
Setting $m_{N_{1(3)}}$, $m_{S_{3(2)}}$ and  $M_{W_R}$ 
to the reference values 
of 300 GeV, 3 TeV, and 5.5 TeV and using the values of 
$\Delta m^2_{31(23)} \cong 2.5\times 10^{-3}~{\rm eV^2}$ 
(see Table \ref{tab2_oscipara}) and of 
${\rm M^{0\nu}_N/M^{0\nu}}$ in Eq. (\ref{eq:MnuMN}) we get:\\ 
$C^{\rm NO(IO)}_{S} \cong (4.05 - 13.90)$ eV for $m_{1(3)} = 10^{-4}$ eV,\\ 
$C^{\rm NO(IO)}_{S} \cong (0.040 - 0.139)$ eV for $m_{1(3)} = 10^{-3}$ eV,\\ 
$C^{\rm NO(IO)}_{S} \cong (4.21\times 10^{-4} - 1.45\times 10^{-3})$ eV for 
$m_{1(3)} = 10^{-2}$ eV,\\ 
$C^{\rm NO(IO)}_{S} \cong (3.24\times 10^{-5} - 1.11\times 10^{-4})$ eV for 
$m_{1(3)} = 5\times 10^{-2}$ eV,\\ 
$C^{\rm NO(IO)}_{S} \cong (2.0 - 6.9)\times 10^{-5})$ eV for 
$m_{1(3)} = 10^{-1}$ eV.

It follows from these numerical estimates that $C^{\rm NO(IO)}_{S}$, and thus
$|m^S_{\beta\beta,R}|$, decreases  rapidly with the increase of 
$m_{1(3)}$ in the interval  $(10^{-4} - 5\times 10^{-2})$ eV.

We recall that the contributions due to the exchange of $S_{1.2.3}$ 
in the NO and IO cases are given by
\begin{eqnarray} 
\left| m^S_{\beta\beta,R} \right|
&=& C^{\rm NO}_{\rm S} \left|U_{e1}^2 +
U_{e2}^2 e^{i\, \alpha}\,\frac{m^3_1}{m^3_{2}}
+ U_{e3}^2 e^{i\, \beta}\, \frac{m^3_1}{m^3_3}
\right|\,~~~\text{ NO case}\,, 
\label{eq:mSNO} 
\end{eqnarray}
\begin{eqnarray}
\left| m^S_{\beta\beta,R} \right|
&=& C^{\rm IO}_{\rm S} 
\left|U_{e1}^2 \frac{m^3_3}{m^3_{1}}+
U_{e2}^2 e^{i\, \alpha}\,\frac{m^3_3}{m^3_{2}}
+ U_{e3}^2 e^{i\, \beta}\,
\right|\,~~~\text{ IO case}\,, 
\label{eq:mSIO} 
\end{eqnarray}
and that $|U_{e1}|^2 \cong 0.7$ and 
$|U_{e3}|^2 \cong 0.022$.

Consider the NO case. 
We note first that with the increasing of $m_1$ beyond $10^{-2}$ eV, 
the contribution $| m^S_{\beta\beta,R}|$ to the 
$0 \nu\beta\beta$ decay amplitude becomes sub-dominant and 
negligible. 
For $m_{1} = (10^{-4} - 10^{-2})$ eV, 
the ratio  $m^3_1/m^3_3 \ll 1$, 
while $m^3_1/m^3_2 \ll 1$ in the interval 
$m_{1} = (10^{-4} - 4.5\times 10^{-3})$ eV. 
This implies that for $m_1 \ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$<$}\  4.5\times 10^{-3}$ eV,
the second and third terms in the expression
(\ref{eq:mSNO}) for $| m^S_{\beta\beta,R}|$
are practically negligible  
and $| m^S_{\beta\beta,R}| \cong  C^{\rm NO}_{\rm S} |U_{e1}|^2$ 
with essentially no dependence on the Majorana phases.
In the interval $m_1 = (4.5\times 10^{-3} - 10^{-2})$ eV 
the ratio  $m^3_1/m^3_2$ increases with $m_1$ and 
at $m_1 =  10^{-2}$ eV we have  $m^3_1/m^3_2 \cong 0.44$. 
Thus,  even at this value the maximal effect of the Majorana phases 
is to change the value of $|U_{e1}|^2\cong 0.7$ to  
$|U_{e1}|^2 \pm  m^3_1/m^3_2 |U_{e2}|^2$, or 
to $0.70 \pm 0.13$, i.e. by at most 18\%, 
where we have used the best fit values of $\sin^2\theta_{12}$ and 
$\sin^2\theta_{13}$. Thus,  in the interval of values 
of $m_1$ of interest, where the contribution of $| m^S_{\beta\beta,R}|$ 
is important, there can not be significant compensation between 
the three terms in the expression for $| m^S_{\beta\beta,R}|$.

From the numerical estimates of 
$| m^\nu_{\beta\beta,L}|$, $| m^N_{\beta\beta,R}|$ and $| m^S_{\beta\beta,R}|$ 
in the preceding and current sub-sections, it follows that 
in the interval of interest $m_{1} = (10^{-4} - 10^{-2})$ eV, 
we have $| m^S_{\beta\beta,R}| > (\gg) | m^\nu_{\beta\beta,L}|,| m^N_{\beta\beta,R}|$.
This is particularly important in the interval 
  $ 10^{-3}~{\rm eV} < m_{1} < 10^{-2})~{\rm eV}$, where 
$| m^\nu_{\beta\beta,L}| <  3\times 10^{-4}$ eV, 
while $| m^S_{\beta\beta,R}|\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  3\times 10^{-4}$ eV and, depending on 
the NME, at  $m_{1} = 10^{-3}$ eV 
can be as large as  $| m^S_{\beta\beta,R}|\cong 9.7\times 10^{-2}~{\rm eV} 
\gg | m^\nu_{\beta\beta,L}|,| m^N_{\beta\beta,R}|$.

The situation is very different in the IO case.
In the interval  $m_{3} = (10^{-4} - 10^{-2})$ eV, where 
the factor  $C^{\rm IO}_{\rm S}$ has a relatively large value,  
we have $(m_3/m_{2(1)})^3 \lesssim 7.5\times 10^{-3}$, 
which implies that actually 
$| m^S_{\beta\beta,R}| \cong C^{\rm IO}_{\rm S}\,|U_{e3}|^2\cong  
 2.2\times 10^{-2}\, C^{\rm IO}_{\rm S}$. 
As a consequence of the suppression due to 
$|U_{e3}|^2$ we have in the interval of values of 
 $m_{3}$ of interest $| m^S_{\beta\beta,R}| \ll |m^N_{\beta\beta,R}|$.

\subsection{The contribution of the 
interference term}
\label{sec:Nexch}
The contribution of the interference term  
$2 \text{Re} (m^{N}_{\beta \beta, R} \cdot m^{S^{*}}_{\beta \beta, R})$
in Eq. (\ref{eq:int}) in the  $0\nu\beta\beta$ decay rate
may be non-negligible only in the  interval of values of 
$m_{1(3)} = (10^{-4} - 10^{-2})$ eV,
where the new non-standard contributions are significant.
Using the analytical expressions for 
$| m^N_{\beta\beta,R}|$ and $| m^S_{\beta\beta,R}|$ 
in Eqs. (\ref{eq:M1}) - (\ref{eq:M4}) and the results 
reported in the preceding subsections, it is not difficult to 
estimate the relative magnitude of the contribution of  
this term. Our results show that 
it varies significantly with the type of neutrino mass spectrum,   
the values of the lightest neutrino mass $m_{1(3)}$
and of the Majorana phases $\alpha$ and $\beta$.

 The relative contribution of the interference term of interest 
is determined by the ratio: 
\begin{equation}
R \equiv \frac{2\text{Re} (m^{N}_{\beta \beta, R} \cdot m^{S^{*}}_{\beta \beta, R})}
{| m^\nu_{\beta\beta,L}|^2 +
|m^{N}_{\beta \beta, R}|^2 + |m^{S}_{\beta \beta, R}|^2}\,.
\label{eq:intover}
\end{equation}
Using the ratio $R$, the generalised effective Majorana mass 
defined in Eq. (\ref{eq:int}) can be written as:
 \begin{equation}
 m^{\nu+|N+S|}_{ee} = m^{\nu+N+S}_{ee}\,\sqrt{1 + R}\,.
\label{eq:intR}
\end{equation}

In the case of NO spectrum, the sign 
of the interference term of interest depends on the
Majorana phases $\alpha$ and $\beta$. 
For $\alpha = \beta = 0$, the ratio $R < 0$  
and thus, the interference terms give a negative contribution to the 
$0\nu\beta\beta$ decay rate. The magnitude of this contribution 
increases quickly when $m_1$ increases from $10^{-4}$ eV to $10^{-3}$ eV
with $R$ changing from (-0.044) to (-0.48). The effect of the 
interference term peaks at $m_1 \cong 2\times 10^{-3}$ eV  
where $R \cong - 0.85$. Thus, at this value of $m_1$ 
we have the maximal suppression of 
$m^{\nu+N+S}_{ee}$ by the factor $\sqrt{1 + R}$:
$m^{\nu+|N+S|}_{ee} \cong  0.39\,m^{\nu+N+S}_{ee}$.
The ratio $R$  decreases rapidly 
when $m_1$ increases beyond $5\times 10^{-3}$ eV 
at which $R\cong -0.48$. The quoted values of $R$ 
at  $m_1 = 10^{-4}$ eV and $m_1 = 10^{-3}$ eV  
are essentially independent of the value of the ratio 
$M^{0\nu}_N/M^{0\nu}_\nu$ lying in the reference interval 
(22.2 - 76.3). The value of $R$ quoted at  $m_1 = 5\times 10^{-3}$ eV  
corresponds to $M^{0\nu}_N/M^{0\nu}_\nu = 76.3$; 
for  $M^{0\nu}_N/M^{0\nu}_\nu = 22.2$ it is significantly  smaller 
in magnitude: $R\cong -0.089$.

The effect of the interference term is quite different 
for $\alpha =\pi$, $\beta = 0$. In this case, 
the interference terms give a positive contribution to the 
$0\nu\beta\beta$ decay rate for $m_1 < 2.26\times 10^{-3}$ eV 
where $R > 0$. At $m_1 \cong 2.26\times 10^{-3}$ eV 
it goes through zero ($R = 0$) since at this value 
$m^\nu_{\beta\beta,L} \cong 0$ and thus  $m^N_{\beta\beta,R} \cong 0$.
Correspondingly, at $m_1 \cong 2.26\times 10^{-3}$ eV 
the generalised effective Majorana mass (Eq. (\ref{eq:int}))
$m^{\nu+|N+S|}_{ee} \cong m^S_{\beta\beta,R} \cong C^{\rm NO}_{\rm S} U_{e1}^2$. 
Taking into account that 
$U_{e1}^2 \cong \cos^2\theta_{12} \cong 0.7$ 
and using Eq. (\ref{eq:CNS2}), for the reference values of 
$m_{N_{1}} = 300$ GeV, $m_{S_{3}} = 3$ TeV,  $M_{W_R} = 5.5$ TeV 
and $M^{0\nu}_N/M^{0\nu}_\nu = 22.2~(76.3)$ 
we find $m^{\nu+|N+S|}_{ee} \cong 5.3\,(18.8) \times 10^{-3}$ eV.
At  $m_1 > 2.26\times 10^{-3}$ eV the interference term is 
negative ($R < 0$). It increases in magnitude as $m_1$ increases 
in the interval $m_1 = 3.5\times 10^{-3} - 10^{-2}$ eV and, e.g.,
at $m_1 = 10^{-2}$ eV we have 
$R \cong -\,0.02~(-\,0.16)$ for 
$M^{0\nu}_N/M^{0\nu}_\nu = 22.2~(76.3)$. At  $m_1 > 10^{-2}$ eV 
we have $|R| \ll 1$ and the interference term has a sub-leading 
(practically negligible) contribution in the  $0\nu\beta\beta$ 
decay rate.

 The results for the ratio $R$ of interest are very different in the 
IO case. It is maximal in magnitude at $m_3 = 10^{-4}$ eV, 
where $R\cong -\,0.54$. However, for $m_3 \sim 10^{-4}$ eV,
the predicted values of the generalised effective Majorana mass 
$|m^{\nu+|N+S|}_{ee}|$ (see Eq. (\ref{eq:int}), 
as we will show in the next Section,  
are strongly disfavored (practically ruled out) by the existing upper limits 
from the KamLAND-Zen and GERDA experiments (see Table  \ref{tab:half-life}
and Fig. 2). In the region of values of  $m_3 \ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  10^{-3}$ eV, 
where the predictions for the generalised effective Majorana mass 
are compatible with the current experimental upper limits 
one has  $|R| < 0.06$, with the value of 
$|R|$ decreasing rapidly with the increasing of $m_3$.
Thus, in the IO case, the interference term under discussion 
has at most, a sub-leading (practically negligible) effect  
on the  $0\nu\beta\beta$ half-life in the 
interval of values of $m_3$ where the predictions of the model considered 
are compatible with the existing lower limits on the half-life.

\subsection{Numerical Results}
\noindent

It follows from the analyses performed in the preceding 
four subsections in particular, that in the NO case 
the contribution due to the  $S_{1,2,3}$
exchange, $|m^S_{\beta\beta,R}|$,  
dominates over the light neutrino $\nu_i$ and 
$N_{1,2,3}$ exchange contributions 
for $10^{-4}~{\rm eV} \leq m_1 \lesssim 1.5\times 10^{-3}$ eV. 
As a consequence, in the indicated interval of values of 
$m_1$ the generalised effective Majorana mass 
$|m^{\nu+|N+S|}_{ee}|$ exhibits weak dependence
on the Majorana phases $\alpha$ and $\beta$
since $|m^S_{\beta\beta,R}|$ practically does not depend 
on these phases. At $m_1 \ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  2\times 10^{-3}$ eV, 
for $\alpha = \beta = 0$,
the $S_{1,2,3}$ contribution is subleading and 
\begin{equation}
 m^{\nu+|N+S|}_{ee} \cong
\sqrt{| m^\nu_{\beta\beta,L}|^2 + |m^N_{\beta\beta,R}|^2}  
\cong | m^\nu_{\beta\beta,L}|\,
\left (1 + \frac{C_N}{m_{N_{1}}m_{1}}\right )^{\frac{1}{2}}\,,
~~{\rm NO}\,,~\alpha = \beta = 0\,,
\label{eq:GEMMIO}
\end{equation}
where we have used Eq. (\ref{eqn:mnNOIO}).  
For $\alpha = \pi$, $\beta = 0$, however, 
$|m^\nu_{\beta\beta,L}|$ is strongly suppressed 
in the interval $m_1 \cong (1.5\times 10^{-3}  - 9\times 10^{-3})$ eV 
and goes through zero at  $m_1 \cong 2.26\times 10^{-3}$ eV,
where the value of $m_1$ is obtained using 
the best fit values of the neutrino oscillations parameters. 
Therefore $|m^S_{\beta\beta,R}|$ gives significant contribution to 
$m^{\nu+|N+S|}_{ee}$ in the indicated interval and determines 
the minimal value of $m^{\nu+|N+S|}_{ee}$. 
At $m_1 \cong 2.26\times 10^{-3}$ eV, e.g., we have 
$|m^{\nu+|N+S|}_{ee}| \cong |m^S_{\beta\beta,R}| \cong C^{\rm NO}_{\rm S}|U_{e1}|^2$.

In contrast, in the IO case 
the contribution due to the  $S_{1,2,3}$ exchange   
$|m^S_{\beta\beta,R}|$ and of the interference term 
$ 2\text{Re} (m^{N}_{\beta \beta, R} \cdot m^{S^{*}}_{\beta \beta, R})$
in the interval of values of $m_3$ of interest 
are practically negligible. Thus, for 
the generalised effective Majorana mass is given by:
\begin{equation}
 m^{\nu+|N+S|}_{ee} \cong
\sqrt{| m^\nu_{\beta\beta,L}|^2 + |m^N_{\beta\beta,R}|^2}  
\cong | m^\nu_{\beta\beta,L}|\,
\left (1 + \frac{C_N}{m_{N_{3}}m_{3}}\right )^{\frac{1}{2}}\,,~~{\rm IO}\,.
\label{eq:GEMMIO}
\end{equation}
In this case $|m^{\nu+|N+S|}_{ee}|$ depends 
significantly on the Majorana phases.

The conclusions regarding the new non-standard contributions 
due to the exchange of virtual heavy Majorana fermions 
$N_{1,2,3}$ and $S_{1,2,3}$ to the $0\nu\beta\beta$ decay  
generalised effective Majorana mass and 
half-life reached in the phenomenological 
analysis are confirmed by our numerical results. 
\begin{figure*}[!ht]
\centering
\includegraphics[width=0.31\textwidth]{EffectiveMass_SM.pdf} 
\hspace*{0.1cm}
\includegraphics[width=0.31\textwidth]{EffectiveMass_SMNRSL_070223.pdf} 
\hspace*{0.1cm}
 \includegraphics[width=0.31\textwidth]{EffectiveMass_SMNRSL_Int_070223.pdf} 
\hspace*{0.1cm}
\includegraphics[width=0.31\textwidth]{SM_Halflife.pdf}
\hspace*{0.1cm}
\includegraphics[width=0.31\textwidth]{SMNRSL_Halflife070223.pdf} 
\hspace*{0.1cm}
\includegraphics[width=0.31\textwidth]{SMNRSL_Int_Halflife070223.pdf}

\caption{Plots showing effective Majorana mass parameter (upper panel) 
and half-life (lower panel) of $0\nu\beta\beta$ decay as functions of 
the lightest neutrino mass $m_{1(3)}$ in the case of NO (IO) 
light neutrino mass spectrum.
The left upper panel shows the dependence on $m_{1(3)}$ of the 
standard mechanism effective Majorana mass, while in the 
middle and right upper panels, the dependence of 
the generalised effective Majorana mass in which the contributions 
due to the exchange of the heavy Majorana fermions $N_{1,2,3}$ and 
$S_{1,2,3}$ are included without accounting for (middle upper panel) 
and accounting for (upper right panel) their interference.
The brown and blue bands correspond respectively to the NO (NH) and IO (IH) 
types of light neutrino mass spectrum. The overlap of the two
bands indicates the region of the QD spectrum.
The lower panels show the dependence on $m_{1(3)}$ of the
$0\nu\beta\beta$ decay half-lives corresponding to the 
respective upper panels. The green horizontal band represents 
the current bound on effective Majorana 
mass from the experiments KamLAND-Zen and GERDA as given in 
Table~\ref{tab:half-life}, whereas the vertical pink bands represent the
bound corresponding to the upper limit 
on the sum of light neutrino masses of 0.12 eV 
reported by the Planck \cite{Planck:2018vyg} and the prospective 
bound of 0.20 eV that can be set by the  KATRIN~\cite{Aker:2019uuj}
experiment. See text for further details.
}
\label{plot:onubb2}
\end{figure*}
\noindent
These are illustrated in Fig. \ref{plot:onubb2}.
In the three upper panels of Fig.  \ref{plot:onubb2} 
we show i) $m^\nu_{ee} \equiv |m^\nu_{\beta\beta,L}|$ (left panel),
ii) $ m^{\nu+N+S}_{ee} \equiv \sqrt{|m^\nu_{\beta\beta,L}|^2 + 
| m^N_{\beta\beta,R}|^2 + | m^S_{\beta\beta,R}|^2}$ (middle panel), and 
iii)  $ m^{\nu+|N+S|}_{ee} \equiv
\sqrt{| m^\nu_{\beta\beta,L}|^2 + |m^N_{\beta\beta,R} + m^S_{\beta\beta,R}|^2}$  
(right panel), as functions of the lightest neutrino mass $m_{1(3)}$ 
in the case of NO (IO) light neutrino mass spectrum.
Thus, the upper left  panel shows the dependence on $m_{1(3)}$ of the 
standard mechanism effective Majorana mass, while in the upper 
middle and right panels, the dependence of 
the generalised effective Majorana mass (GEMM) in which the contributions 
due to the exchange of the heavy Majorana fermions $N_{1,2,3}$ and 
$S_{1,2,3}$ are included without accounting for (middle panel) 
and accounting for (right panel) their interference.
The brown and blue bands correspond respectively to the NO (NH) and IO (IH) 
types of light neutrino mass spectrum. The overlap of the two
bands indicates the region of the QD spectrum.
Following the discussion of NME in Section 4, 
the ratio of nuclear matrix elements 
${\cal M}^{0\nu}_N/{\cal M}^{0\nu}_\nu$ is varied in the interval 22.2 - 76.3,
as given in Eq. (\ref{eq:MnuMN}).
The minimal (maximal) value in this interval, 
${\cal M}^{0\nu}_N/{\cal M}^{0\nu}_\nu = 22.2~(76.3)$, 
corresponds to ${\cal M}^{0\nu}_\nu  = 4.68~(5.26)$. 
The results for ${\cal M}^{0\nu}_N/{\cal M}^{0\nu}_\nu = 22.2~(76.3)$ 
are indicated with solid (dashed) lines.
For the parameters  $M_{W_R}$, $m_{N_{1(3)}}$ and $m_{S_{3(2)}}$ 
the reference values of 5.5 TeV, 300 GeV, and 3TeV, respectively,
are used. All plots are obtained by varying the neutrino oscillation 
parameters in their respective 3$\sigma$ allowed ranges.  
The Majorana phases $\alpha$ and $\beta$ are varied in the interval 
$[0,\pi]$, while the Dirac phase $\delta$ is set to zero.
For both ${\cal M}^{0\nu}_N/{\cal M}^{0\nu}_\nu = 22.2$ and 76.3,
the curves showing the maximal (minimal) values of 
GEMM as functions of $m_{1(3)}$ 
correspond to $\alpha = \beta = 0$ ($\alpha =\pi$, $\beta = 0$).
The lower panels show the dependence on $m_{1(3)}$ of the
$0\nu\beta\beta$ decay half-lives corresponding to the 
respective upper panels. 
The green horizontal band represents the current bound on effective Majorana 
mass from the experiments KamLAND-Zen and GERDA as given in 
Table~\ref{tab:half-life}, whereas the vertical pink bands represent the
bound corresponding to the upper limit 
on  the sum of light neutrino masses of 0.12 eV 
reported by the Planck experiment \cite{Planck:2018vyg} and the prospective 
bound of 0.20 eV that can be set by the  KATRIN~\cite{Aker:2019uuj}
experiment. 

 A comparison between the upper left and right panels in 
Fig. \ref{plot:onubb2} shows that the presence of the new non-standard 
contributions change drastically the dependence  
of the effective Majorana mass of the standard mechanism 
$|m^\nu_{ee}|$ on the lightest neutrino mass $m_{1(3)}$ 
at $m_{1(3)} <  10^{-2}$ eV, where the new 
contributions dominate over the standard contribution.
At  $m_{1(3)} \ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\   10^{-2}$ eV the new contributions 
are strongly suppressed and practically negligible and
we have   $|m^{\nu+|N+S|}_{ee}|\cong |m^\nu_{ee}|$, 
as also is clearly seen in Fig. \ref{plot:onubb2}. 
 
 The non-standard contributions are so large 
at relatively small values of the  $m_{1(3)}$ that 
for $m_{1(3)} \lesssim 2\times 10^{-4}$ eV 
in the NO case they are ruled out even for the minimal value 
of $M^N_{0\nu}/M^\nu_{0\nu} = 22.2$ 
by the existing upper limits from the KamLAND-Zen and GERDA 
experiments. In the IO case the new contributions are 
also ruled our  for $M^N_{0\nu}/M^\nu_{0\nu} = 76.3$; 
for  $M^N_{0\nu}/M^\nu_{0\nu} = 22.2$ they are ruled out for 
$\alpha = \beta = 0$, while for $\alpha = \pi$, $\beta = 0$, 
they are compatible with the  KamLAND-Zen and GERDA 
upper limits. 

 For NO spectrum and $\alpha = \beta = 0$, 
the inequality  $|m^{\nu+|N+S|}_{ee}|> (>>) |m^\nu_{ee}|$
always holds in the interval of values of
 $m_{1} \cong  (3\times 10^{-4} - 8\times 10^{-3})$ eV 
where the non-standard contributions are significant.
In  this interval and for  $M^N_{0\nu}/M^\nu_{0\nu} = 76.3~(22.2)$, 
$|m^{\nu+|N+S|}_{ee}|\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  0.009~(0.007)$ eV, 
and with the exception of a very narrow interval 
around $m_{1} \cong 4.0~(2.5)\times 10^{-3}$ eV at which the 
quoted minimum of $|m^{\nu+|N+S|}_{ee}|$ takes place,  
we have $|m^{\nu+|N+S|}_{ee}|\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  0.010$ eV.  
In most of the considered intervals of values of 
$m_1$ the half-life $T^{\nu +|N_R + S_L|}_{1/2} \lesssim 10^{28}$ yrs. 

 In the case of NO spectrum, $\alpha =\pi$, $\beta = 0$ and 
$M^N_{0\nu}/M^\nu_{0\nu} = 76.3~(22.2)$,  
the value of $|m^{\nu+|N+S|}_{ee}|\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  0.010$ eV, 
and  $T^{\nu +|N_R + S_L|}_{1/2} \lesssim 10^{28}$ yrs,
at $m_{1} \lesssim 2.0~(1.5)\times 10^{-3}$ eV.
For the minimal value of  $|m^{\nu+|N+S|}_{ee}|$ 
we find ${\rm min}(|m^{\nu+|N+S|}_{ee}|) \cong 9~(3)\times 10^{-4}$ eV.
It takes place at $m_{1} \cong 9.0~(8.0)\times 10^{-3}$ eV.
We recall that  $|m^\nu_{ee}|$ goes through zero at 
$m_{1} \cong 2.26\times 10^{-3}$ eV, while 
$|m^{\nu+|N+S|}_{ee}|\cong 18.3~(5.3)\times 10^{-3}~$ eV 
at this value of $m_1$.
 
For IO spectrum, we have 
$|m^{\nu+|N+S|}_{ee}|\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  0.015$ eV 
for $m_3 < 9\times 10^{-3}$ eV, where 
$|m^{\nu+|N+S|}_{ee}| >  |m^\nu_{ee}|$ 
for any of the considered values of $M^N_{0\nu}/M^\nu_{0\nu}$
and of $\alpha$ and $\beta$. The approximate 
equality  $|m^{\nu+|N+S|}_{ee}| \cong  |m^\nu_{ee}|$ 
holds at $m_3\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  10^{-2}$ eV.
For all considered values of $m_3$,  $M^N_{0\nu}/M^\nu_{0\nu}$ and 
the Majorana phases the predicted half-life  
$T^{\nu +|N_R + S_L|}_{1/2} \lesssim 10^{28}$ yrs, while  
in the case of $\alpha = \beta = 0$, 
we have  $T^{\nu +|N_R + S_L|}_{1/2} \lesssim 2\times 10^{27}$ yrs. 

It follows from our numerical analysis  that most of the parameter 
space of the considered model, the predictions for the 
$0\nu\beta\beta$ decay generalised effective Majorana 
mass and half-life are within the sensitivity range of the 
planned next generation of neutrinoless double beta decay 
LEGEND-200 (LEGEND-1000), nEXO, KamlAND-Zen-II, CUPID, NEXT-HD 
(see \cite{Giuliani:2019uno,Agostini:2022zub} 
and references quoted therein).


\section{Comments on LFV and LHC signatures}
\label{sec:comments}

The considered model has rich lepton flavour violating (LFV) 
and collider phenomenology. A detailed investigation of the 
model's phenomenology is beyond the scope of the present study.
We limit ourselves here with a few brief comments.

  The LFV processes as like 
$\mu \rightarrow e + \gamma$, $\mu \rightarrow 3e$ decays and 
$\mu - e$ conversion in nuclei  can be mediated by heavy RH and 
sterile neutrinos $N_{1,2,3}$ and $S_{1,2,3}$. 
Although we expect the contributions due to  $N_{1,2,3}$ and especially 
due to $S_{1,2,3}$ to be rather suppressed, there might be 
a relatively large region of the 
model's parameter space where they might still be in the range of 
sensitivity of the next generation of experiments 
MEG II, Mu3e, Mu2e, COMET and PRISM/PRIME 
(see, e.g., \cite{Calibbi:2017uvl} and the references therein).

At LHC, the main channel for the production of the heavy RH neutrinos $N_{1,2,3}$ is
via on-shell $Z_R$ production and $W_R$ fusion and can be expressed as 
$p+p\rightarrow W_R^{\pm}\rightarrow l^{\pm} + N_j$, $l=e,\mu,\tau$.
This $N_j$ further decays as 
$N_j\rightarrow W_R^*\rightarrow l^{\prime \pm} + 2j$, $l^{\prime}=e,\mu,\tau$,
which is considered as the ``smoking gun'' signature of lepton number 
and lepton flavour violation at LHC. 
This rapid decay of $N_j$ happens in the case its mass is 
sufficiently large. Our model satisfies this requirement
as we have taken  ${\rm max}(M_{N_j})\sim 100$ GeV.
We recall that the mass of $W_R$ is constrained by experiments CMS,ATLAS and 
low energy precision measurements as $M_{W_R}\gtrsim 5$ TeV ~\cite{ATLAS:2018dcj,ATLAS:2019isd,CMS:2018agk,Li:2020wxi,Dekens:2021bro} and considering 
the relation $M_{Z_R}\simeq 1.2 M_{W_R}$, the mass of $Z_R$ can be constrained 
as $M_{Z_R}\gtrsim 6$ TeV. 
If the mass of $N_j$ lies in the range $5-20$ GeV, 
then it takes some time to decay and travels some distance resulting in a 
displaced vertex of leptons \cite{Helo:2013esa,Izaguirre:2015pga}. 
So, the observable in this case would be a prompt charged lepton and 
a displaced leptonic vertex. The current status of displaced vertex searches 
at LHC can be found in ref. \cite{ATLAS:2012av,CMS:2013czn,ATLAS:2012cdk}.  
Another distinguishing feature in the signatures of small mass ($<100$ GeV) 
and large mass ($\sim 800$ GeV) RH neutrinos $N_j$ is the angle between 
the produced charged leptons. In the former case, parallel tracks of charged 
leptons are expected, whereas in the later case, back-to-back emissions 
are expected \cite{Almeida:2000pz}.




\section{Summary}
\label{sec:concl}

In the present article, we have derived predictions for 
the neutrinoless double beta ($0\nu\beta\beta$)
decay generalised effective Majorana mass
and half-live in a left-right (L-R) symmetric model with 
the double seesaw mechanism at the TeV scale.
The gauge group of the model is the standard L-R symmetric 
extension of the Standard Model (SM) 
gauge group: $\mathcal{G}_{LR}\equiv SU(2)_L \times SU(2)_R \times U(1)_{B-L}$.
The fermion sector has the usual 
for the L-R symmetric models, three families of 
left-handed (LH) and right-handed (RH) 
quark and lepton fields, including right-handed neutrino fields $N_{\beta R}$,
$\beta = e,\mu,\tau$, assigned respectively to  
$SU(2)_L$  and $SU(2)_R$ doublets. It included also   
three $SU(2)_{L,R}$ singlet LH fermion fields $S_{\gamma L}$.
The Higgs sector is composed of two $SU(2)_L$  and $SU(2)_R$
Higgs doublets $H_L$ and $H_R$, and of a bi-doublet  $\Phi$.
The vacuum expectation value (VEV) of  the  $SU(2)_R$ Higgs doublet 
$H_R$ breaks the  $\mathcal{G}_{LR}$ gauge symmetry  
to the SM gauge symmetry  $SU(2)_L \times U(1)_{Y}$, 
while the VEVs of the two neutral components of the 
bi-doublet $\Phi$ break  $SU(2)_L \times U(1)_{Y}$ to $U(1)_{\rm em}$.
The Yukawa couplings of the LH and RH fermion doublets 
to the bi-doublet $\Phi$ generate (via the VEVs of 
the neutral components of $\Phi$) Dirac mass terms for the 
quarks and charged leptons, as well as 
a $\nu_{\alpha L} - N_{\beta R}$ Dirac mass term 
$M^\nu_{\rm D}$ involving the 
LH active flavor neutrino fields 
$\nu_{\alpha L}$ and the RH fields $N_{\beta R}$, 
$\alpha,\beta=e,\mu,\tau$. 
The singlet LH fermion fields $S_{\gamma L}$ are assumed to have a 
Majorana mass term
$M_{\rm S}$ and Yukawa coupling 
with the RH doublets containing $N_{\beta R}$ which involves $H_R$. 
This Yukawa coupling produces  a $S_{\gamma L} - N_{\beta R}$ Dirac mass term 
$M_{\rm RS}$ when  $H_R$ develops a non-zero VEV.
Under the condition $|M_{\rm RS}| \ll |M_{\rm S}|$,
the RH neutrinos  $N_{\beta R}$ get a Majorana mass term 
$M_{\rm R} \cong -\,M_{\rm RS}M^{-1}_{\rm S} M^T_{\rm RS}$ 
via a seesaw-like mechanism.
This in turn generates a Majorana mass term 
for the LH flavour neutrinos 
$m_\nu \cong -\, M^\nu_{\rm D}M^{-1}_{\rm R}(M^\nu_{\rm D})^T$
via a second seesaw mechanism 
\footnote{Hence the term ``double or cascade seesaw mechanism 
of neutrino mass generation''.}. 
In such a way, the model contains in addition to the three 
light Majorana neutrinos $\nu_i$ having masses $m_i$, 
two sets of heavy Majorana 
particles $N_j$ and $S_k$, $j,k=1,2,3$, with masses 
$m_{N_j} \ll m_{S_k}$. The double seesaw scenario allows 
the RH neutrinos $N_j$ to have masses naturally  
at the GeV-TeV scale.

In our analysis of the  $0\nu\beta\beta$
decay predictions of the model, we have 
considered the case of $m_{N_j} \sim (1 - 1000)$ GeV 
and ${\rm max}(m_{S_k}) \sim (1 - 10)$ TeV, $m_{N_j} \ll m_{S_k}$. 
Working with a specific version of the model 
which can be obtained by employing symmetry arguments and 
in which the 
Dirac mass terms  $M^\nu_{\rm D}$ and $M_{\rm RS}$
are diagonal, $M^\nu_{\rm D} = k_d\,{\bf I}$, $M_{\rm RS} = k_{rs}\,{\bf I}$,
$k_d$ and $k_{RS}$ being constant mass parameters and 
{\bf I} the $3\times 3$ unit matrix, 
we have studied in detail the new ``non-standard'' contributions to the 
$0\nu\beta\beta$ decay amplitude and half-life 
arising from diagrams with an exchange of virtual $N_j$ and $S_k$.
The self-consistency of the considered set-up 
requires that the lightest neutrino mass for the neutrino mass 
spectrum with normal ordering (NO), $m_1$, 
or with inverted ordering (IO), $m_3$, 
has to be not smaller that approximately $10^{-4}$ eV.
Moreover, the RH neutrino ($N_{\beta R}$) and sterile fermion  
($S_{\gamma L}$) mixings are determined 
by the light neutrino PMNS mixing matrix.   
In the analysis of the new non-standard contributions to the 
 $0\nu\beta\beta$ decay amplitude we took into account 
the values of the nuclear matrix elements (NMEs) 
 ${\cal M}^{0\nu}_N$ and ${\cal M}^{0\nu}_\nu$
 associated with, respectively, the light and heavy Majorana neutrino 
exchange contributions, calculated for the four isotopes 
$^{76}$Ge,  $^{82}$Se,  $^{130}$Te, $^{136}$Xe
by six different groups of authors using different 
methods of NME calculation (Table \ref{tab:nucl-matrix}). 
We made use of the fact that the ratio  
${\cal M}^{0\nu}_N/{\cal M}^{0\nu}_\nu$ 
reported by each of the six cited groups 
is essentially the same for the considered four isotopes -- 
it varies with the isotope by not more than 
$\sim 15\%$. For a given isotope 
the ratio of interest obtained by the six different 
methods of NME calculation 
varies by a factor of up to $\sim 3.5$. 
In view of this, we took into account the 
uncertainties in the NME calculations  by using 
the following reference range of the ratio 
${\cal M}^{0\nu}_N/{\cal M}^{0\nu}_\nu = 22.2 - 76.3$, 
which corresponds to  $^{76}$Ge.


 We analyzed in detail the properties 
of the new non-standard contributions to the 
$0\nu\beta\beta$ decay amplitude arising due to the exchange 
of virtual heavy Majorana fermions $N_j$ and $S_k$, 
parametrized as effective Majorana masses 
$m^N_{\beta\beta,R}$ (Eq. \ref{eqn:mnNOIO}) and $m^S_{\beta\beta,R}$ 
(Eq. \ref{eq:mSNO}), respectively. 
These analyses showed that 
both $|m^N_{\beta\beta,R}|$ and $|m^S_{\beta\beta,R}|$ 
are strongly enhanced at relatively small values of 
the lightest neutrino mass $m_{1(3)}\sim (10^{-4} - 8\times 10^{-3})$ eV.
The effect of this enhancement is particularly 
important in the case of NO neutrino mass spectrum. 
The non-standard contributions are so large 
at the indicated small values of  $m_{1(3)}$ that 
for $m_{1} \lesssim 2\times 10^{-4}$ eV 
in the NO case, they are strongly disfavored (if not ruled out) 
even for the minimal value 
of $M^N_{0\nu}/M^\nu_{0\nu} = 22.2$ 
by the existing upper limits from the KamLAND-Zen and GERDA 
experiments. In the IO case, the new contributions are 
also stronlgly disfavored  for $M^N_{0\nu}/M^\nu_{0\nu} = 76.3$; 
for  $M^N_{0\nu}/M^\nu_{0\nu} = 22.2$ they are 
disfavored for $\alpha = \beta = 0$, while for $\alpha = \pi$, $\beta = 0$, 
they are still compatible with the  KamLAND-Zen and GERDA 
conservative upper limits. 
We find, in general,  that in both NO and IO cases
the new  non-standard contributions due to $N_j$ and $S_k$ exchange 
are dominant over the standard light neutrino exchange contribution  at 
values of the lightest neutrino mass $m_{1(3)} \sim (10^{-4} - 10^{-2})$ eV:
$|m^{\nu+|N+S|}_{ee}| > (>>) |m^\nu_{ee}|$, 
where  $m^{\nu+|N+S|}_{ee}$ is the generalised effective Majorana mass (GEMM) 
which accounts for all contributions to the $0\nu\beta\beta$ decay amplitude 
(Eqs. (\ref{halflife:no-int}),  (\ref{eqn:mee_SR})
and   (\ref{eq:int})), and $m^\nu_{ee}$ is the effective Majorana 
mass associated with the standard light neutrino exchange contribution 
(Eq. (\ref{effmassnu})). The effective Majorana mass $|m^S_{\beta\beta,R}|$ 
associated with $S_k$ exchange contribution was shown to be practically 
independent of the Majorana phases $\alpha$ and $\beta$,  
while that due to exchange of  $N_j$, $|m^N_{\beta\beta,R}|$,
exhibits strong dependence on $\alpha$ and $\beta$ 
similar to $|m^\nu_{ee}|$. 

 For NO spectrum and $\alpha = \beta = 0$, 
the inequality  $|m^{\nu+|N+S|}_{ee}|> (>>) |m^\nu_{ee}|$
always holds in the interval of values of
 $10^{-4}~{\rm eV} \lesssim m_{1} \lesssim 8\times 10^{-3}$ eV 
where the non-standard contributions are significant.
In  this interval and for  $M^N_{0\nu}/M^\nu_{0\nu} = 76.3~(22.2)$, 
$|m^{\nu+|N+S|}_{ee}|\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  0.009~(0.007)$ eV.
With the exception of a very narrow interval 
around $m_{1} \cong 4.0~(2.5)\times 10^{-3}$ eV at which the 
quoted minimum of $|m^{\nu+|N+S|}_{ee}|$ takes place,  
we have $|m^{\nu+|N+S|}_{ee}|\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  0.010$ eV.  
In most of the considered intervals of values of 
$m_1$ the $0\nu\beta\beta$ decay  half-life 
$T^{\nu +|N_R + S_L|}_{1/2} \lesssim 10^{28}$ yrs. 

In the case of NO spectrum, $\alpha =\pi$, $\beta = 0$ and 
$M^N_{0\nu}/M^\nu_{0\nu} = 76.3~(22.2)$,  
we find that $|m^{\nu+|N+S|}_{ee}|\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  0.010$ eV, 
and  $T^{\nu +|N_R + S_L|}_{1/2} \lesssim 10^{28}$ yrs,
at $m_{1} \lesssim 2.0~(1.5)\times 10^{-3}$ eV.
For the minimal value of  $|m^{\nu+|N+S|}_{ee}|$ 
we get ${\rm min}(|m^{\nu+|N+S|}_{ee}|) \cong 9~(3)\times 10^{-4}$ eV.
It takes place at $m_{1} \cong 9.0~(8.0)\times 10^{-3}$ eV.
We note that  $|m^\nu_{ee}|$ goes through zero at 
$m_{1} \cong 2.26\times 10^{-3}$ eV, while 
$|m^{\nu+|N+S|}_{ee}|\cong 18.8~(5.3)\times 10^{-3}~$ eV 
at this value of $m_1$. Thus, the strong suppression of the 
$0\nu\beta\beta$ decay rate at $m_{1} \sim 2.26\times 10^{-3}$ eV 
and in the interval $m_1 \cong (1.5 - 8.0)\times 10^{-3}$ eV 
in the case of only standard contribution due to the exchange of 
light Majorana neutrinos $\nu_i$ 
(Fig. \ref{plot:onubb2}, upper left panel)
is avoided due to the new non-standard contributions. 
 
For IO spectrum, we find that 
$|m^{\nu+|N+S|}_{ee}|\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  0.015$ eV 
for $m_3 < 9\times 10^{-3}$ eV, where 
$|m^{\nu+|N+S|}_{ee}| >  |m^\nu_{ee}|$ 
for any of the considered values of $M^N_{0\nu}/M^\nu_{0\nu}$
and of $\alpha$ and $\beta$. The approximate 
equality  $|m^{\nu+|N+S|}_{ee}| \cong  |m^\nu_{ee}|$ 
holds at $m_3\ \raisebox{-.4ex}{\rlap{$\sim$}} \raisebox{.4ex}{$>$}\  10^{-2}$ eV.
For all considered values of the parameters 
the predicted half-life  
$T^{\nu +|N_R + S_L|}_{1/2} \lesssim 10^{28}$ yrs, while  
in the case of $\alpha = \beta = 0$, 
we have  $T^{\nu +|N_R + S_L|}_{1/2} \lesssim 2\times 10^{27}$ yrs. 

It follows from our results that in most of the parameter space of 
the considered model, the predictions for the 
$0\nu\beta\beta$ decay generalised effective Majorana 
mass and half-life are within the sensitivity range of the 
planned next generation of neutrinoless double beta decay 
experiments LEGEND-200 (LEGEND-1000), nEXO, KamlAND-Zen-II, CUPID, NEXT-HD 
(see \cite{Giuliani:2019uno,Agostini:2022zub} 
and references quoted therein).


\section*{Acknowledgements}
Purushottam Sahu would like to acknowledge the Ministry of Education, 
Government of India for financial support. 
PS also acknowledges the support from the Abdus Salam International Centre 
for Theoretical Physics (ICTP) under the ``ICTP
Sandwich Training Educational Programme (STEP)'' SMR.3676 and SMR.3799.
The work of S. T. P. was supported in part by the European
Union's Horizon 2020 research and innovation program under the 
Marie Sk\l{}odowska-Curie grant agreement No.~860881-HIDDeN, by the Italian 
INFN program on Theoretical Astroparticle Physics and by the World Premier 
International Research Center Initiative (WPI Initiative, MEXT), Japan.
STP would like to thank Kavli IPMU, University of Tokyo, 
where part of this study was performed for the kind hospitality. 


\section{Appendix :  \\
Derivation of neutrino masses and mixings in left-right double 
seesaw model (LRDSM)}
\label{app:lrdssm}
\subsection{LRDSM mass matrix}
We discuss here the implementation and derivation of double seesaw mechanism in the considered left-right symmetric model. The neutral fermions needed for LRDSM are active left-handed neutrinos, $\nu_L$, active right-handed neutrinos, $N_R$ and sterile neutrinos, $S_L$. The relevant mass terms are given by 
\begin{eqnarray}
\mathcal{L}_{\rm LRDSM} &=& \mathcal{L}_{M_D}+\mathcal{L}_{M_{RS}} 
+ \mathcal{L}_{M_S} \nonumber \\
\mathcal{L}_{M_D} &=& - \sum_{\alpha, \beta} \overline{\nu_{\alpha L}} [M_D]_{\alpha \beta} N_{\beta R} \mbox{+ h.c.}\nonumber \\
\mathcal{L}_{M_{RS}} &=&  \sum_{\alpha, \beta} \overline{S_{\alpha L}} [M_{RS}]_{\alpha \beta} N_{\beta R} \mbox{+ h.c.}
\nonumber \\
\mathcal{L}_{M_{S}} &=&  \frac{1}{2} \sum_{\alpha, \beta} \overline{S^c_{\alpha L}} [M_{S}]_{\alpha \beta} S_{\beta L} \mbox{+ h.c.}
\end{eqnarray}
The flavour states for active left-handed neutrinos $\nu_{\alpha L}$, right-handed neutrinos $N_{\beta R}$ and sterile neutrinos $S_{\gamma L}$ are defined as follows
\begin{eqnarray}
\nu_{\alpha L} = \begin{pmatrix}
\nu_{eL} \\ \nu_{\mu L} \\ \nu_{\tau L}
\end{pmatrix}\, , ~ 
N_{\beta R} = \begin{pmatrix}
N_{1 R} \\ N_{2 R} \\ N_{3 R}
\end{pmatrix}\, , ~
S_{\gamma L} = \begin{pmatrix}
S_{1 L} \\ S_{2 L} \\ S_{3 L}
\end{pmatrix}\,  
\end{eqnarray}
Similarly, their mass states can be written as,
\begin{eqnarray}
\nu_{i L} = \begin{pmatrix}
\nu_{1L} \\ \nu_{2 L} \\ \nu_{3L}
\end{pmatrix}\, , ~ 
N^c_{j R} = \begin{pmatrix}
N^c_{1 R} \\ N^c_{2 R} \\ N^c_{3 R}
\end{pmatrix}\, , ~ 
S_{k L} = \begin{pmatrix}
S_{1 L} \\ S_{2 L} \\ S_{3 L}
\end{pmatrix}\,
\end{eqnarray}


The $9 \times 9$ neutral lepton mass matrix in 
the basis $\left(\nu_L, N^c_R, S_L\right)$ is given by
\begin{equation}
 \mathcal{M}_{LRDSM}=
 \left[ 
\begin{array}{c | c} 
  \begin{array}{c c} 
     \blue{\bf 0} & M_D \\ 
     M^T_D & \blue{\bf 0}
  \end{array} & 
  \begin{array}{c} 
     \red{0} \\ M_{RS}
  \end{array} \\ 
\hline 
  \begin{array}{c c} 
     \,~\red{0}  &\quad  M^T_{RS}
  \end{array} & 
  \begin{array}{c} 
     \blue{M_S}
  \end{array} \\
 \end{array} 
\right]
\label{eqn:dss-a}       
\end{equation}
where each elements of the matrix is a $3 \times 3$ matrix. Here $M_D$ is the Dirac neutrino mass matrix connecting $\nu_L - N_R$, $M_{RS}$ is the mixing matrix in the $N_R-S_L$ sector, $M_S$ is the Majorana mass matrix for sterile neutrino $S_L$. 

In order to diagonalize the above mass matrix we have used the following masse hierarchy,
\begin{equation}
 M_D < M_{RS} < M_{S} .
\end{equation}

The diagonalisation of $\mathcal{M_{\rm LRDSM}}$ after changing it from flavour basis to mass basis is done by a generalized unitary transformation as, 
\begin{eqnarray}
&&\mid \Psi \rangle_{\rm flavor} = V\mid \Psi \rangle_{\rm mass}\\
&\mbox{or,}& \begin{pmatrix}
\nu_{\alpha L}\\ N^c_{\beta R}\\ S_{\gamma L}
\end{pmatrix}
= 
\begin{pmatrix}
V_{\alpha i}^{\nu \nu} & V_{\alpha j}^{\nu N} & V_{\alpha k}^{\nu S}\\V_{\beta i}^{ N \nu} & V_{\beta j}^{NN} & V_{\beta k}^{N S}\\V_{\gamma i}^{S \nu} & V_{\gamma j}^{S N} & V_{\gamma k}^{SS}
\end{pmatrix}               
\begin{pmatrix}
\nu_{i}\\ N^c_{j}\\ S_{k}
\end{pmatrix} \, \\
&&   V^{\dagger} \mathcal{M_{\rm LRDSM}} V^* 
= \mathcal{\widehat{M}}_{\rm LRDSM} \nonumber \\
&&\hspace*{2.2cm} = \mbox{diag} \left(m_{i},m_{N_j},m_{S_k} \right) \nonumber \\
&&\hspace*{2.2cm} = \mbox{diag} \left(m_{1},m_{2},m_{3},m_{N_1},m_{N_2},m_{N_3},m_{S_1},m_{S_2},m_{S_3} \right)
\end{eqnarray}
Here the indices $\alpha, \beta, \gamma$ run over three generations of light left-handed neutrinos, heavy right-handed neutrinos and sterile neutrinos respectively, whereas the indices $i,j,k$ run over corresponding mass states. 

\subsection{Block Diagonalization of double seesaw mass matrix in LRDSM}
 
\label{app:lrdsm}
Let's write the matrix in eq.\ref{eqn:dss-a} as,

\begin{eqnarray}
\mathcal{M}_{LRDSM}=\mathcal{M_{\nu}}=\begin{pmatrix}
              \mathcal{M}_L & \mathcal{M}_D  \\
              \mathcal{M}^T_D & \mathcal{M}_S
             \end{pmatrix} \,, ~ \mbox{where,} \nonumber \\
   \mathcal{M}_L = \begin{pmatrix}
                                  0 & M_D \\
                                  M^T_D  & 0
                                 \end{pmatrix}
                           , 
                           \mathcal{M}_D=  \begin{pmatrix}
                                  0 \\ M_{RS}
                                 \end{pmatrix},\mathcal{M}_S = M_S\,
\end{eqnarray}
The complete block diagonalization is achieved in two steps by recursively integrating out the heavier modes as
\begin{eqnarray}
\mathcal{W}_{1}^{\dagger} \mathcal{M_{\nu}} \mathcal{W}_{1}^* 
      =\mathcal{\hat M_{\nu}^\prime} &\mbox{and}& \mathcal{W}_{2}^{\dagger} \mathcal{\hat M_{\nu}^\prime} \mathcal{W}_{2}^* 
      =\mathcal{\hat{M}_{\nu}}
\end{eqnarray}
where $\mathcal{\hat M_{\nu}^\prime}$ is block diagonalised $9\times9$ matrix after integrating out the heaviest mode and $\mathcal{\hat{M}_{\nu}}$
is the block diagonalised $9\times9$ matrix after integrating out the next heaviest mode. The transformation matrix $\mathcal{W}_{1}$ can be written as a general unitary matrix in the form
\begin{eqnarray}
\mathcal{W}_{1}^* &=& \begin{pmatrix}
              \sqrt{1-\mathcal{B}\mathcal{B}^{\dagger}}& \mathcal{B}  \\
              -\mathcal{B}^{\dagger}& \sqrt{1-\mathcal{B}^{\dagger}\mathcal{B}}
             \end{pmatrix} 
\end{eqnarray} where $\mathcal{B}$ is a $ 6 \times 3 $ dimentional matrix.
\begin{eqnarray}
\sqrt{1-\mathcal{B}\mathcal{B}^{\dagger}} &=& 1 - \frac{1}{2}\mathcal{B}\mathcal{B}^{\dagger}-\frac{1}{8}\left(\mathcal{B}\mathcal{B}^{\dagger}\right)^2
+ \cdots \nonumber \\
\mathcal{B} = \sum \mathcal{B}_i 
\end{eqnarray}
At leading order it looks like 
\begin{eqnarray}
\sqrt{1-\mathcal{B}\mathcal{B}^{\dagger}} &\simeq& 1 - \frac{1}{2}\mathcal{B}\mathcal{B}^{\dagger}-\frac{1}{8}\left(\mathcal{B}_1\mathcal{B}_2^{\dagger}+\mathcal{B}_2\mathcal{B}_1^{\dagger}\right)
\end{eqnarray}
The form of mixing matrix $\mathcal{B}_1^{\dagger}$ and $\mathcal{B}_1$ is given by
\begin{eqnarray}
\mathcal{B}_1^{\dagger}=\mathcal{M}_S^{-1} \cdot \mathcal{M^T_D} &=& M_S^{-1}\cdot\begin{pmatrix}
                                  0 & M_{RS}^T
                                 \end{pmatrix} 
                                 = \begin{pmatrix}
                                  0 &M^{-1}_S \cdot M_{RS}^T
                                 \end{pmatrix}  \nonumber\\
\mathcal{B}_1 &=& \begin{pmatrix}
                                  0 \\
                                M_{RS}M^{-1}_S
                                 \end{pmatrix} \nonumber\\
\sqrt{1-\mathcal{B}\mathcal{B}^{\dagger}} &\simeq& \begin{pmatrix}
                                  1& 0\\
                                  0  & 1-\frac{1}{2}M_{RS}M^{-1}_S  \cdot M^{-1}_S M^T_{RS}
                                 \end{pmatrix}\nonumber\\
\sqrt{1-\mathcal{B}^{\dagger}\mathcal{B}} &\simeq& 1-\frac{1}{2} M^{-1}_S M^T_{RS} \cdot  M_{RS}M^{-1}_S                                 
\end{eqnarray}
Thus, the first block diagonalised mixing matrix $\mathcal{W}_1$ becomes,
\begin{eqnarray}
\mathcal{W}_1 = \begin{pmatrix}
                                  1 & 0 & 0\\
                                  0 &1-\frac{1}{2}M_{RS}M^{-1}_S \cdot M^{-1}_S  M_{RS}^T& M_{RS}M^{-1}_S \\
                                  0& -M^{-1}_SM_{RS}^T & 1-\frac{1}{2} M^{-1}_SM_{RS}^T \cdot M_{RS}M^{-1}_S
                                 \end{pmatrix}
 \end{eqnarray}

After this diagonalization, $\mathcal{\hat M_{\nu}^\prime}$ has the following form.
\begin{eqnarray}
\mathcal{\hat M_{\nu}^\prime} &=& \begin{pmatrix}
                                  \mathcal{M}_{eff} & 0\\
                                  0  & \mathcal{M}_S
                                 \end{pmatrix}
 \end{eqnarray}
where,
\begin{eqnarray}
\mathcal{M}_{eff} &=& \mathcal{M}_L - \mathcal{M}_D \mathcal{M}^{-1}_S \mathcal{M}^T_D \nonumber \\
&=&\begin{pmatrix}
                     0 & M_D \\
                     M^T_D   & 0 
                    \end{pmatrix} 
    -  \begin{pmatrix}
         0 \\
         M_{RS}
        \end{pmatrix} M^{-1}_S
                  \begin{pmatrix}
                   0 & M^T_{RS} 
                  \end{pmatrix} \nonumber \\
&=&   \begin{pmatrix}
       0    & M_D  \\
       M^T_D          &  -M_{RS} M^{-1}_S M^T_{RS}
       \end{pmatrix} \, 
\end{eqnarray}
$\mathcal{M}_{eff}$ can be further diagonalised by $\mathcal{W}_2$ as,
\begin{eqnarray}
\mathcal{S}^{\dagger} \mathcal{M}_{eff} \mathcal{S}^* &=&\begin{pmatrix}
       m_{\nu}   & 0 \\
       0       &  M_R
       \end{pmatrix},
\end{eqnarray}
where, 
\begin{eqnarray}
 m_\nu &=& - M_D \left( -M_{RS} M^{-1}_S M^T_{RS} \right)^{-1} M^T_D , \nonumber \\
 && M_R = -M_{RS} M^{-1}_S M^T_{RS}.
\end{eqnarray}
The transformation matrix $\mathcal{S}$ is
\begin{eqnarray}
\mathcal{S}^* &=& \begin{pmatrix}
              \sqrt{1-\mathcal{A}\mathcal{A}^{\dagger}}& \mathcal{A}  \\
              -\mathcal{A}^{\dagger}& \sqrt{1-\mathcal{A}^{\dagger}\mathcal{A}}
             \end{pmatrix}              
\end{eqnarray}
such that, 
\begin{eqnarray}
\mathcal{A}^{\dagger}&=&\left( -M_{RS} M^{-1}_S M^T_{RS} \right)^{-1} M_D \\
&& = -M^{-1^T}_{RS} M_S M^{-1}_{RS}M_{D} = X^{\dagger}
\end{eqnarray}
\begin{eqnarray}
\text{Thus,} ~\mathcal{W}_{2} &=& \begin{pmatrix}
  \mathcal{S}  & 0\\
   0  &  1
\end{pmatrix}\\ &=& \begin{pmatrix}
1-\frac{1}{2}X X^{\dagger}  & X &0\\
   -X^{\dagger}  & 1-\frac{1}{2}X^{\dagger}X & 0 \\
   0 &0 & 1

\end{pmatrix} \nonumber \\
\end{eqnarray}

\subsection{Complete diagonalization and physical neutrino masses}
After block diagonalization, the mass matrix for the three types of neutrinos are further diagonalized by respective unitary mixing matrices $U_\nu$, $U_N$, $U_S$ resulting in physical masses for the neutrinos as follows.

\begin{eqnarray}
U_{9 \times 9} = \begin{pmatrix}
{ U_\nu}_{3 \times 3}  &  { 0}_{3 \times 3}  & { 0}_{3 \times 3}  \\
{ 0}_{3 \times 3}  &  {U_N}_{3 \times 3}    & { 0}_{3 \times 3}  \\
{ 0}_{3 \times 3}  &  { 0}_{3 \times 3}  & {U_S}_{3 \times 3}
\end{pmatrix} 
\end{eqnarray}

\begin{eqnarray}
&&U^\dagger_\nu m_\nu U^*_\nu = \hat{m}_\nu = \mbox{diag}\left(m_{\nu_1}, m_{\nu_2}, m_{\nu_3} \right) \nonumber \\
&&U^\dagger_N M_N U^*_N = \hat{M}_N = \mbox{diag}\left(M_{N_1}, M_{N_2}, M_{N_3} \right) \nonumber \\
&&U^\dagger_S M_S U^*_S = \hat{M}_S = \mbox{diag}\left(M_{S_1}, M_{S_2}, M_{S_3} \right)
\end{eqnarray}

The complete mixing matrix  now becomes,
\begin{eqnarray}
V &=& \mathcal{W}_1\cdot \mathcal{W}_2 \cdot \mathcal{U}\nonumber\\
&=& \begin{pmatrix}
1  & 0 &0\\
   0  & 1-\frac{1}{2} Y Y^{\dagger} & Y \\
   0 &-Y^{\dagger} & 1-\frac{1}{2}  Y^{\dagger} Y
 
\end{pmatrix}
\cdot
\begin{pmatrix}
1-\frac{1}{2}X X^{\dagger}  & X &0\\
   -X^{\dagger}  & 1-\frac{1}{2}X^{\dagger}X & 0 \\
   0 &0 & 1
   \end{pmatrix} 
   \cdot \begin{pmatrix}
U_{\nu}  & 0&0\\
   0 & U_N& 0 \\
   0 &0 & U_S
   \end{pmatrix}  \nonumber\\
&=& \begin{pmatrix}
U_{\nu}\left(1-\frac{1}{2}X X^{\dagger} \right)  & U_{N}X &0\\
   -U{\nu}X^{\dagger}\left(1-\frac{1}{2}YY^{\dagger} \right)  & U_N \left(1-\frac{1}{2}X^{\dagger}X \right) \left(1-\frac{1}{2}YY^{\dagger} \right) & U_{S} Y \\
   U{\nu}X^{\dagger}Y^{\dagger} &-U_NY^{\dagger} & U_S\left(1-\frac{1}{2}Y^{\dagger}Y \right)
\end{pmatrix}   
\label{mixingmatrix1}
\end{eqnarray}

where $ X^{\dagger} = -M^{-1}_{RS} M_S M^{-1}_{RS}M_{D}$, $Y = M_{RS} M^{-1}_{S}$ and fixing the typical magnitudes for $ M_D \simeq $~0.1\,MeV, $ M_{RS} \simeq $~1\,TeV, $M_{S} \simeq$ 10\,TeV we get $X \simeq 10^{-6} $, $Y\simeq 0.1$.  Since $U_{\nu}$, $U_N$ and $U_S$ are of $\mathcal{O}(1)$, the matrix elements of $\mathcal{V}$ are approximated to be
\begin{eqnarray}
\begin{pmatrix}
\mathcal{V}_{\alpha i}^{\nu \nu} & \mathcal{V}_{\alpha j}^{\nu N} & \mathcal{V}_{\alpha k}^{\nu S}\\\mathcal{V}_{\beta i}^{ N \nu} & \mathcal{V}_{\beta j}^{NN} & \mathcal{V}_{\beta k}^{NS}\\\mathcal{V}_{\gamma i}^{S \nu} & \mathcal{V}_{\gamma j}^{SN} & \mathcal{V}_{\gamma k}^{SS}
\end{pmatrix}  \simeq \begin{pmatrix}
{\cal O}(1.0)       & {\cal O}(10^{-6})     & 0  \\
{\cal O}(10^{-6})      & {\cal O}(1.0)    & {\cal O}(0.1)    \\
{\cal O}(10^{-7})   & {\cal O}(0.1)    & {\cal O}(1.0) 
\end{pmatrix}
\end{eqnarray} 
which generates sizable contribution to neutrinoless double beta decay.




\clearpage
\bibliographystyle{utcaps_mod}
