Experiments have firmly established that the neutrinos, which are massless within the Standard Model (SM) of particle physics, have non-zero, albeit tiny mass, and different flavours of neutrinos mix with each other (see [1,2,3] for the current best fit values of the parameters). Further, neutrinos are electrically neutral, allowing them to be their own anti-particles ie can be Majorana in nature [4]. Neutrinoless double beta decay (\(0\nu 2\beta \,\)), \((A,Z)\rightarrow (A,Z+2) + 2e^-\), provides an unambiguous way of establishing the Majorana nature of the neutrinos, and also the lepton number violation [5, 6]. Theoretically as well, \(0\nu 2\beta \,\)decay is heralded as a useful probe of physics beyond SM, having particular relevance for neutrino masses and mass hierarchy. For an incomplete list discussing \(0\nu 2\beta \,\)decay phenomenology and related signatuures see e.g. [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48]. The search for neutrinoless double beta decay thus constitutes an important endevour. Experimentally, extensive studies have been carried out or are planned on several nuclei [49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72]. Currently, the best limits on the half life of nuclei are: \(T_{1/2}(^{76}Ge > 1.8\times 10^{26}\,\, \mathrm {yr.\,\, at\,\,} 90\%\,\, \textrm{CL}\) [67], and \(T_{1/2}(^{136}Xe > 2.3\times 10^{26}\,\, \mathrm {yr.\,\, at\,\,} 90\%\,\, \textrm{CL}\) [72]. These limits are expected to improve by an order or two in future.

Neutrinoless double beta decay process can proceed via the light neutrinos, \(\nu _i\)’s, (the so called long range part) or due to heavy degrees of freedom (the short range part) for example through exchange of heavy neutrinos, \(N_i\)’s, or other heavy particles in specific models like R-parity violating supersymmetric theories or theories with leptoquarks (see [73,74,75,76,77,78,79,80,81,82] and references therein for a quick review of some of the essential theoretical and experimental issues). The short range part due to the heavy physics is due to intermediate particles with masses much larger than the relevant scale of the process \(\sim {\mathcal {O}}\)(GeV), allowing for the heavier degrees of freedom to be systematically integrated out, leaving behind a series of operators built out of low energy fields, the up and down quarks and electrons, weighted by the short distance coefficients, called Wilson coefficients (denoted by \(C_A\) below). This provides a very convenient and systematic framework to evaluate the decay amplitude in terms of short distance coefficients which encode all the information about the high energy physics. In this process of integrating out the heavy degrees of freedom, the operators and thus the effective Lagrangian obtained is at the typical scale of the heavy particles. Using then the renormalization group equations (RGEs), perturbative QCD effects can be computed. These QCD corrections have been shown to be very significant [83, 84], in particular due to the colour mismatched operators (see also [85,86,87,88,89]). In this way, the high energy particle physics input gets separated from the low energy dynamics contained in the nuclear matrix elements (NMEs) of the quark level operators sandwiched between the nucleon states (see [90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108] to get an idea of different approaches to calculate these NMEs). These NMEs are usually the source of large uncertainty to \(0\nu 2\beta \,\)predictions, and in principle one could compare predictions for various nuclei undergoing \(0\nu 2\beta \,\)transition to understand and eventually reduce the sensitivity on NMEs.

To be concrete, we begin by considering the short range quark-electron operators (denoted by \({\mathcal {L}}_{\mathrm qe}\) below in the text) and as an example, consider a heavy right handed neutrino, mass \(M_N\), and SM gauge group. The resulting quark level \(0\nu 2\beta \,\)amplitude takes the form

$$\begin{aligned} {\mathcal {A}}\sim & {} \underbrace{\frac{(V_{ud} T_{ei})^2}{M_W^4M_N}}_{G}\underbrace{\bar{u}\gamma _{\mu }(1-\gamma _5)d\, \bar{u}\gamma ^{\mu }(1-\gamma _5)d}_{{\mathcal {J}}_{q,\mu }{\mathcal {J}}_q^{\mu }}\, \underbrace{\bar{e}(1+\gamma _5)e^c}_{j_e} \end{aligned}$$

where \(V_{ud}\) and \(T_{ei}\) are the quark and leptonic mixing matrix elements while \(G = G_F^2/M_N\) contains the short distance physics, with \(G_F\) being the Fermi constant. The physical \(0\nu 2\beta \,\)amplitude is obtained by sandwiching the quark level operators between the initial and final nuclear states, \(\vert i\rangle \) and \(\vert f\rangle \), finally evaluated in terms of the NMEs:

$$\begin{aligned} {\mathcal {A}}_{0\nu 2\beta } = \langle f\vert i{\mathcal {H}}_{\mathrm eff}\vert i\rangle \sim G\, \underbrace{\langle f\vert {\mathcal {J}}_{q,\mu }{\mathcal {J}}_q^{\mu }\vert i\rangle }_{\varvec{NME}}\,j_e \end{aligned}$$

The short distance or high energy physics cleanly separates from the low energy matrix elements. The low energy effective Lagrangian is expressed as a sum of operators, \(O_A\) weighted by the Wilson coefficents \(C_A\): \({\mathcal {L}}_{\mathrm eff} = G_A C_A O_A\), where we have allowed for more than one G for more complicated theories. In the example considered above, there is only one operator \(O_1 = {\mathcal {J}}_{q,\mu }{\mathcal {J}}_q^{\mu }\,j_e = \bar{u_i}\gamma _{\mu }(1-\gamma _5)d_i\,\bar{u_j}\gamma ^{\mu }(1-\gamma _5)d_j\, \bar{e}(1+\gamma _5)e^c\) (i, j denoting the colour indices) and the corresponding Wilson coefficient \(C_1=1\). In other models like SUSY with R-parity violation or leptoquarks, Fierz transformations have to be employed to bring the operators in form similar to above. The Lorentz and Dirac structure of the quark level operator involved decides which NME enters the \(0\nu 2\beta \,\)rate. Perturbative QCD corrections don’t just correct the Wilson coefficients but also lead to colour mismatched operators (typically with strength that is \(1/N_c\) (\(N_c=3\) being the number of colours) of the colour allowed operators) which are then Fierz transformed and can give different operators and thereby bringing a host of different NMEs which would not have been expected otherwise.

In the usual treatment of the short range part (see [73,74,75,76,77,78,79,80,81,82]), one considers following simplifications: (i) assume that the final electrons are emitted in the S-wave ie the long wavelength approximation is employed; (ii) closure approximation ie the energy of the virtual neutrino is much larger than the nulcear excitation energy, thereby allowing to sum over the intermediate set of states with great ease; (iii) (non-relativistic) impulse approximation which allows to write the matrix elements of product of the quark currents between the nuclear states in terms of say the Fermi and Gammow–Teller matrix elements. For the case of quark currents being V-A form, the above procedure implies that

$$\begin{aligned} {\mathcal {J}}_{q,\mu }(\mathbf {x_1}){\mathcal {J}}_q^{\mu }(\mathbf {x_2})= & {} \sum _{n,m}\tau _+^n\tau _+^m\delta (\mathbf {x_1}-\mathbf {r_n}) \delta (\mathbf {x_2}-\mathbf {r_m}) \nonumber \\{} & {} \times (g_V^2(q^2)-g_A^2(q^2)\mathbf {\sigma ^n}\cdot \mathbf {\sigma ^m}) \end{aligned}$$
(1)

which lead to the Fermi (vector part of the current) and Gammow–Teller (axial-vector part of the current) nuclear matrix elements, \({\mathcal {M}}_{F}\) and \({\mathcal {M}}_{GT}\) respectively

$$\begin{aligned} {\mathcal {M}}_{F}= & {} \left\langle \Psi _f\vert \sum _{n,m} H(r_n,r_m,\bar{E})\tau _+^n\tau _+^m\vert \Psi _i\right\rangle \end{aligned}$$
(2)
$$\begin{aligned} {\mathcal {M}}_{GT}= & {} \left\langle \Psi _f\vert \sum _{n,m} H(r_n,r_m,\bar{E})\tau _+^n\tau _+^m \mathbf {\sigma ^n}\cdot \mathbf {\sigma ^m}\vert \Psi _i\right\rangle \end{aligned}$$
(3)

In the above expressions, \(\sigma \) and \(\tau \) are Pauli matrices and \(\tau _+ = (\tau _1+i\tau _2)/2\); indices n, m run over all the nucleons in the nucleus and \(H(r_n,r_m,\bar{E})\) denotes the heavy neutrino potential. Further, \(g_V(0)=1\) and \(g_A(0)\simeq 1.27\). \(\Psi _{i,f}\) are the initial and final state nuclear wave-functions. It is worth mentioning that in writing these matrix elements starting from the quark level currents, the dipole and anapole terms are not shown as they are much smaller than those displayed above. For the present purpose, we shall choose to neglect them but they can be systematically included. See [90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108] for details and updated numerical values of the NMEs calculated within different schemes. The differences in the numerical values of NMEs are at the heart of large uncertainties in the predictions.

Considering only vector and axial quark currents, the following set of short distance operators are to be considered (showing only the colour matched operators; i, j in the operators below denote the colour indices)

$$\begin{aligned} O^{LL}= & {} \bar{u_i}\gamma _{\mu }(1-\gamma _5)d_i\,\bar{u_j}\gamma ^{\mu }(1-\gamma _5)d_j\,\bar{e}(1+\gamma _5)e^c \nonumber \\ O^{RR}= & {} \bar{u_i}\gamma _{\mu }(1+\gamma _5)d_i\,\bar{u_j}\gamma ^{\mu }(1+\gamma _5)d_j\,\bar{e}(1+\gamma _5)e^c \nonumber \\ O^{LR}= & {} \bar{u_i}\gamma _{\mu }(1-\gamma _5)d_i\,\bar{u_j}\gamma ^{\mu }(1+\gamma _5)d_j\,\bar{e}(1+\gamma _5)e^c \nonumber \\ O^{RL}= & {} \bar{u_i}\gamma _{\mu }(1+\gamma _5)d_i\,\bar{u_j}\gamma ^{\mu }(1-\gamma _5)d_j\,\bar{e}(1+\gamma _5)e^c \end{aligned}$$
(4)

with their associated Wilson coefficients, \(C^{LL,RR,LR,RL}\). These are the set of operators say in left–right theories, with LR, RL stemming from the heavy-light W mixing. From an effective field theory (EFT) point of view, there will be scalar-pseudoscalar and tensor operators (and their combinations) as well. For definiteness, we focus on the operators in Eq. (4). It turns out that the Short Range (SR) NMEs depend on whether the quark chiralities are same or different:

$$\begin{aligned} {\mathcal {M}}_{SR}^{LL} = {\mathcal {M}}_{SR}^{RR}= & {} g_V^2{\mathcal {M}}_{F} - g_A^2{\mathcal {M}}_{GT} \nonumber \\ {\mathcal {M}}_{SR}^{LR+RL}= & {} g_V^2{\mathcal {M}}_{F} + g_A^2{\mathcal {M}}_{GT} \end{aligned}$$
(5)

It is to be noted that there is a very large spread in the determinations of the NMEs. In particular, what is relevant here, \(-{\mathcal {M}}_{GT}\) can be larger than \({\mathcal {M}}_{F}\) by a factor (see for example [90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108]) to almost an order of magnitude (see [109, 110]). Below, for simplicity and convenience, we shall assume that the two terms on the right hand side of the above equations differ by a factor \(\sim 5\), after taking into account \(g_V\) and \(g_A\). Since the discussion will not be specific to a particular nucleus, we shall assume this for all the nuclei. Therefore, one has the approximate results: \({\mathcal {M}}_{SR}^{LL} = {\mathcal {M}}_{SR}^{RR} \sim 6{\mathcal {M}}_{F}\) and \({\mathcal {M}}_{SR}^{LR+RL} \sim -4{\mathcal {M}}_{F}\), thereby yielding the following approximate forms for the corresponding contributions to the amplitude:

$$\begin{aligned} {\mathcal {A}}_{SR}^{LL/RR}\sim & {} 6\,C^{LL/RR}\,{\mathcal {M}}_{F} \nonumber \\ {\mathcal {A}}_{SR}^{LR/RL}\sim & {} -4\,C^{LR/RL}\,{\mathcal {M}}_{F} \end{aligned}$$
(6)

It is to be kept in mind that the above relations are obtained by making a choice for the NMEs (\({\mathcal {M}}_{GT} \sim -5{\mathcal {M}}_{F}\)) to bring out the essence of the argument but when actually computing rates extreme care muct be exercised since the extent to which this relation holds will crucially depend on the specific set of NMEs adopted. Below, whenever a statement to such an effect is made, it is with this approximate relation in mind.

A given theory of lepton number violation implies a set of quark-lepton level operators comprising the effective Lagrangian. Using the RGEs and including effects of operator mixing, the final amplitude for the \(0\nu 2\beta \,\)can be written in terms of short distance coefficients and various NMEs. The rate thus computed can be then contrasted with the experimental limits on the half life of the neutrinoless double beta decay process for the specific nucleus and stringent limits are obtained on the parameters of the theory. Alternatively, an EFT point of view could be adopted and all the relevant operators are then written at the hadronic scale, using which the amplitude and thus the decay rate is computed, which then leads to constraints on different effective coefficients. One of these two is followed in the phenomenological studies. For the purpose of illustration and to be concrete, we’ve considered the short range mechanism originating due to the right handed neutrinos and having SM like interactions. In a different model, there could be other interactions leading to newer operators. A similar analysis would be required. In case there are more than one mechanism present for the neutrinoless double beta decay to proceed, there can be a competition between the different mechanisms. Similarly, there would in general be a need to add the contribution due to the light neutrino exchange. Since the aim here is not a detailed numerical study in a given model, we continue the discussion along the lines adopted thus far ie consider the short range operators listed in Eq. (4).

The above is however not the entire story. Till now, the \(0\nu 2\beta \,\)amplitude is calculated by directly evaluating the nuclear matrix elements from the quark currents. However, it was shown in [111, 112] that there are additional contributions wherein these quark currents first hadronise into pions, and then one can evaluate the two pion exchange contributions which were shown to be the dominant ones. See also [113]. These authors used the on-shell matching conditions to match the matrix elements of quark-lepton currents between the pion states on to the hadron level effective terms. Ref. [114] provided a systematic EFT set-up to match the quark-electron operators to a chiral effective theory with two pion and electron (\(\pi \pi ee\)) vertices as well as two nucleon, one pion and electron (\(NN\pi ee\)) vertices and four nucleon and electron (NNNNee) vertices, where \(\pi \), N and e denote pion, nucleon and electron respectively. It also discussed inclusion of next to leading and next to next to leading order contributions. It was confirmed that the two pion exchange contributions indeed dominate though their calculation employed Naive Dimensional Analysis (NDA) for the power counting in contrast to the Vacuum Insertion Approximation (VIA) as employed in the on-shell matching in [111, 112]. That VIA may miss some fraction of the total contribution is not an unexpected satatement since in VIA, only the contribution due to vacuum is retained while that due to other states is neglected. However, VIA provides a simple and quick approximation to estimate these new effects. The one and two pion contributions are non-local as they involve pion propagators. Recall that \(m_{\pi ^+} \sim 140\,{\textrm{MeV}}\) and the typical momentum, \(\vec {q}\), flowing through the propagators is \({\mathcal {O}}(200\,{\textrm{MeV}})\). On the other hand, the contribution due to NNNNee vertices is a local one. This approach has been further developed and refined in [115,116,117,118,119]. That the two pion contribution can play a rather important role in phenomenological studies has been re-emphasized in [120] in the context of left–right symmetric theories. It is worth pointing out that almost all phenomenological studies on \(0\nu 2\beta \,\)do not include these contributions and as recently pointed out in [120], the inferences and limits drawn could be significantly altered once these are included.

Given that there are important contributions to \(0\nu 2\beta \,\)amplitude beyond those obtained by directly taking the matrix elements of the short distance quark currents between the nuclear states emanating from one and two pion exchange diagrams, one is led to ask if there are additional contributions beyond the pionic ones. In particular, does \(\rho \)-meson exchange, analogous to pion exchange, yield significant contributions? We now focus our attention on new contributions coming from hadronization of quark currents into \(\rho \) mesons. Relevant diagrams are shown in Fig. 1. Compared to the pion exchange, there are two differences: (i) \(m_{\rho }\, (= 770\, {\textrm{MeV}})>> \vert \vec {q}\vert \, (\sim 200\, {\textrm{MeV}})\), (ii) \(m_{\pi }<< m_N \,(\sim 940\, {\textrm{MeV}}) \sim 1.2m_{\rho }\). Statement (i) means that it is a good approximation to neglect \(\vert q\vert \) in comparison with \(m_{\rho }\), both in the numerator and denominator of the \(\rho \)-meson propagator: \(-i(g_{\mu \nu }-q_{\mu }q_{\nu }/m_{\rho }^2)/(q^2-m_{\rho }^2) \rightarrow ig_{\mu \nu }/m_{\rho }^2\) ie at this level of approximation, the \(\rho \) contribution is being treated as a local one. It is to be noted that this additional local contribution should be distinguished from with theconventional short range one, since having the explicit \(\rho \)-meson specific mass and couplings will play an important role, as we see below.

Fig. 1
figure 1

Hadronic level diagrams (drawn using JaxoDraw [124, 125]) with the intermediate dotted lines denoting a \(\rho \) meson. a NNNNee; b two \(\rho \) exchange; c \(NN\rho ee\) at one vertex and other vertex being \(\rho NN\) (possibly including parity violating terms) and the permuted diagram with the vertices exchanged

Full size image

The parity conserving strong \(\rho NN\) vertices are described by the phenomenological Lagrangian (see for example [121])

$$\begin{aligned} {\mathcal {L}}_{\rho NN} = g_{\rho }\bar{N}\left( \gamma _{\mu } + i\frac{\chi _{\rho }}{2m_N}\sigma _{\mu \nu }q^{\nu }\right) \vec {\tau }\cdot \vec {\rho }^{\mu }N \end{aligned}$$
(7)

where again the nucleon anomalous magnetic moment term, with strength \(\chi _{\rho }\) (also denoted as \(\mu _N\) in literature), is neglected as \(q<< m_N\). There exist various determinations of \(g_{\rho }\) and \(\chi _{\rho }\) but the combination \(g_{\rho }(1+\chi _{\rho })\) takes roughly the same value \(\approx 21\). In what follows, we choose \(g_{\rho } \sim 4.5\) (see [122] and references therein). Parity violating terms can also be written. The parity violating couplings are found to be significantly smaller compared to parity conserving ones and therefore, in the present context such terms are not included. It is worthwhile to point out that the coupling analogous to \(g_{\rho }\) for the \(\omega \)-meson is large: \(g_{\omega } \sim 14\) [122]. However, since \(\omega \) is charge neutral, it does not contribute to the \(0\nu 2\beta \,\)amplitude.

Consider first the diagram (b) in Fig. 1. Analogous to the pion case in [111, 112], we write the relevant term in the hadronic Lagrangian as:

$$\begin{aligned} {\mathcal {L}}_{\rho \rho ee} = \frac{G_F^2m_{\rho }^2}{2 m_N}(m_{\rho }^2 a_{2\rho }\vec {\rho }^{\alpha }\cdot \vec {\rho }_{\alpha }) \,\bar{e}(1+\gamma _5)e^c \end{aligned}$$
(8)

Using this vertex, along with the parity conserving \(\rho NN\) vertex and approximating the \(\rho \) propagators as above, the amplitude for diagram in Fig. 1b takes the form

$$\begin{aligned} {\mathcal {A}}_{2\rho } = i\frac{G_F^2m_{\rho }^2}{2 m_N} g_{\rho }^2a_{2\rho }[\bar{u_p}\gamma _{\mu }u_n][\bar{u_p}\gamma ^{\mu }u_n]\,\bar{e}(1+\gamma _5)e^c \end{aligned}$$
(9)

To obtain \(a_{2\rho }\), on-shell matching condition is employed:

$$\begin{aligned} \langle \rho ,2e\vert {\mathcal {L}}_{qe}\vert \rho \rangle = \langle \rho ,2e\vert {\mathcal {L}}_{\rho \rho ee}\vert \rho \rangle \end{aligned}$$
(10)

where to facilitate a quick comparison, the quark-electron Lagrangian is rewritten as

$$\begin{aligned} {\mathcal {L}}_{qe} = \frac{G_F^2}{2 m_N} \sum _A\,C^{A}O^{A}\,\bar{e}(1+\gamma _5)e^c \end{aligned}$$
(11)

with \(A = LL,\, RR,\, LR+RL\). To complete the on-shell matching, the last step needed is to make use of vacuum, dominance or VIA. This is achieved through the use of following defining relations (see [123]):

$$\begin{aligned} \langle 0\vert \bar{u}\gamma _{\mu }\vert \rho (P)\rangle= & {} f_{\rho }m_{\rho }\epsilon _{\mu } \nonumber \\ \langle 0\vert \bar{u}\sigma _{\mu \nu }\vert \rho (P)\rangle= & {} if^T_{\rho }(\epsilon _{\mu }P_{\nu }-\epsilon _{\nu }P_{\mu }) \end{aligned}$$
(12)

with \(f_{\rho } = 198\pm 7\) MeV and \(f^T_{\rho } = 160\pm 10\) MeV. \(\epsilon _{\mu }\) in the above equations denotes the \(\rho \)-meson polarization vector. On-shell matching, assuming VIA, yields

$$\begin{aligned} a_{2\rho } = -\frac{f_{\rho }^2}{m_{\rho }^2}\sum _A C^A \end{aligned}$$
(13)

Plugging this in Eq. (9), the final form of two \(\rho \) exchange diagram after accounting for the combinatoric factors reads

$$\begin{aligned} {\mathcal {A}}_{2\rho }= & {} -i\frac{G_F^2}{2 m_N}\left( 2\frac{g_{\rho }^2f_{\rho }^2}{m_{\rho }^2}\sum _A C^A\right) [\bar{u_p}\gamma _{\mu }u_n][\bar{u_p}\gamma ^{\mu }u_n]\nonumber \\{} & {} \times \bar{e}(1+\gamma _5)e^c \nonumber \\\sim & {} -i\frac{G_F^2}{2 m_N}\left( 2.5\sum _A C^A\right) [\bar{u_p}\gamma _{\mu }u_n][\bar{u_p}\gamma ^{\mu }u_n]\nonumber \\{} & {} \times \bar{e}(1+\gamma _5)e^c \end{aligned}$$
(14)

Next, consider the one \(\rho \) exchange diagram (c) in Fig. 1 and its permuted one. The \(NN\rho ee\) interaction term is written as

$$\begin{aligned} {\mathcal {L}}_{NN\rho ee} = \frac{G_F^2m_{\rho }^2}{2 m_N}a_{1\rho }\,\bar{N}(\gamma _{\mu }\vec {\tau }\cdot \vec {\rho }^{\mu })N \,\bar{e}(1+\gamma _5)e^c \end{aligned}$$
(15)

resulting in one \(\rho \) exchange contribution to the \(0\nu 2\beta \,\)amplitude given by

$$\begin{aligned} {\mathcal {A}}_{1\rho } = -i\frac{G_F^2}{2 m_N} g_{\rho }a_{1\rho }[\bar{u_p}\gamma _{\mu }u_n][\bar{u_p}\gamma ^{\mu }u_n]\,\bar{e}(1+\gamma _5)e^c \end{aligned}$$
(16)

Following the same steps as above for on-shell matching and factorizing the quark level product of currents using VIA as: \(\langle p\vert {\mathcal {J}}_{q,\mu }{\mathcal {J}}_q^{\mu }\vert \rho n\rangle = \langle p\vert {\mathcal {J}}_{q,\mu }\vert n\rangle \langle 0\vert {\mathcal {J}}_q^{\mu }\vert \rho \rangle = g_V (\bar{u}_p\gamma _{\mu }u_n) m_{\rho }f_{\rho }\epsilon ^{\mu }\), one obtains

$$\begin{aligned} a_{1\rho } = \frac{f_{\rho }g_V}{m{\rho }} \sum _A C^A \end{aligned}$$
(17)

Plugged back, one obtains the final amplitude due to one \(\rho \) meson exchange after incorporating the combinatoric factors as

$$\begin{aligned} {\mathcal {A}}_{1\rho }= & {} -i\frac{G_F^2}{2 m_N} \left( \frac{4 g_{\rho }f_{\rho }g_V}{m{\rho }}\sum _A C^A\right) [\bar{u_p}\gamma _{\mu }u_n][\bar{u_p}\gamma ^{\mu }u_n]\nonumber \\{} & {} \times \bar{e}(1+\gamma _5)e^c \nonumber \\\sim & {} -i\frac{G_F^2}{2 m_N}\left( 4.5 \sum _A C^A\right) [\bar{u_p}\gamma _{\mu }u_n][\bar{u_p}\gamma ^{\mu }u_n]\nonumber \\{} & {} \times \bar{e}(1+\gamma _5)e^c \end{aligned}$$
(18)

where \(g_V(0) = 1\) has been used.

The inclusion of the \(\rho \) exchange diagrams thus produces an extra contribution\({\mathcal {A}}^{VIA}_{\rho }\) to \(0\nu 2\beta \,\)

$$\begin{aligned} {\mathcal {A}}^{VIA}_{\rho }= & {} {\mathcal {A}}_{1\rho } + {\mathcal {A}}_{2\rho }\nonumber \\\sim & {} -i\frac{G_F^2}{2 m_N}\left( 7\sum _A C^A\right) [\bar{u_p}\gamma _{\mu }u_n][\bar{u_p}\gamma ^{\mu }u_n]\,\bar{e}(1+\gamma _5)e^c \nonumber \\ \end{aligned}$$
(19)

One immediately notices from the Dirac structure of the amplitudes in Eqs. (9) and (16) that \(\rho \) exchange only results in Fermi transition matrix element. The net effect due to \({\mathcal {A}}_{1\rho }\) and \({\mathcal {A}}_{2\rho }\) is to modify the LL, RR and LR/RL short range amplitudes, Eq. (6):

$$\begin{aligned} {\mathcal {A}}_{SR}^{LL/RR}\rightarrow & {} 6\,C^{LL/RR}\,{\mathcal {M}}_{F} - 7\,C^{LL/RR} {\mathcal {M}}_{F} \sim 0 \nonumber \\ {\mathcal {A}}_{SR}^{LR/RL}\rightarrow & {} -4\,C^{LR/RL}\,{\mathcal {M}}_{F} - 7\,C^{LR/RL} {\mathcal {M}}_{F} \nonumber \\= & {} -11\,C^{LR/RL} {\mathcal {M}}_{F} \end{aligned}$$
(20)

This constitutes the main result of this study, namely an almost complete cancellation brought about due to \(\rho \) exchange contributions for the LL and RR short range amplitudes, \({\mathcal {A}}_{SR}^{LL}\) and \({\mathcal {A}}_{SR}^{RR}\), and the mixed left–right amplitude, \({\mathcal {A}}_{SR}^{LR+LR}\) getting almost tripled. These results tend to indicate that \(\rho \) exchange diagrams can impact the phenomenological predictions significantly. The almost complete cancellation and three fold enhancement observed above hinge on the choices of various couplings and parameters adopted, and caution must be exercised when actually making use of these infereces since the extent to which destructive or constructive effects are observed is bound to change. While there are uncertainties in the values of these couplings and parameters, the broad conclusions are expected to hold though the extent to which the cancellation or enhancement, as the case may be, will change. The NMEs themselves have a large spread but as is clear from the expressions above, the cancellation or enhancement is brought by the numerical factors multiplying these NMEs. Thus, the spread or large uncertainty in NMEs is not going to change this observation. What will change with the spread in the NME values is the final numerical value of amplitude or the half life. This is not the focus of the present work but is definitely an important aspect to be considered while providing precise predictions for the \(0\nu 2\beta \,\)rates.

Before closing, it may be worthwhile to discuss the differences between the contributions that stem from hadronization of the quarks into pions versus rho mesons. As briefly alluded to above in the article, pion mass being comparable to the scale of momentum transfer makes the pion (one pion and two pion exchange) contributions non-local. These contributions are governed by appropriate powers of \(g_{\rho }\) and \(\pi \pi ee\) or \(NN\pi ee\) vertices. The inclusion of these non-local contributions can enhance some of the otherwise suppressed contributions (see eg [120]). In contrast, due to heavier mass, the rho-meson contributions turn out to be effectively local and the significant contribution brought in by such terms is due to multiplicative factors that bring in appropriate factors of \(g_{\rho }\) along with other combinatoric factors (see Eqs. (14) and (18)). Another noteworthy point is that for the \((V-A)(V-A)\) interactions considered here as a concrete example, the \(\rho \)-meson contributes only to the Fermi transition and thus has a distinctive impact. This will change if different interactions are considered. At this point it may be useful to mention that while the pion exchange diagrams are expected to be dominant compared to the rho exchange ones, the importance lies in the fcat that the rho exchange diagrams being effectively local in nature impact the local contributions. For example, \({\mathcal {A}}_{SR}^{LL/RR} \) drastically diminish in magnitude once rho exchange pieces are accounted for. The pionic long distance, non-local contributions will thus provide the amplitude for the short range operator. Also, there will be a competition between the final pion and rho mediated pieces and a priori it is not clear or obvious which will eventually dominate. It will be worth and useful to consider all the short range operators at the quark electron level and work out the impact of \(\rho \)-meson contributions, especially for some specific beyond the standard model scenarios like left–right symmetric model. While a complete study within concrete models including both the pion and \(\rho \)-meson contributions is left for a future study, it may be instructive to have a crude estimate of the expected strength of \(\rho \)-meson contribution relative to the pion contribution. To this end, we consider the LR operator that gives the leading order double-pion contribution (see Eq. (48) in [114]). A rough estimate using two pion contribution as a representative pion exchange contribution and making use of the \(\rho \) meson contribution obtained above, we get (with representative values of the NMEs used both for the pion and \(\rho \) case):

$$\begin{aligned} \vert \frac{{\mathcal {A}}^{(\rho )}}{{\mathcal {A}}^{(2\pi )}}\vert \sim 0.1 \end{aligned}$$
(21)

This is certainly a very crude estimate but serves to show that for certain operators the \(\rho \) exchange contribution may not be completely insignificant. This then strengthens the case for a more systematic study of such contributions. The above calculations and estimate comparing the \(\rho \) contribution to the pion contribution reveal two pieces: while the \(\rho \) contribution at the amplitude level can be few tens of percent compared to pion contribution, the \(\rho \) contribution (for the choices of couplings etc assumed) can result in large cancellation for the LL or RR type operators and enhancement for mixed operators LR. Thus, say for the LL or RR operators, it just leaves the pion contribution as the sole contribution coming from the dimension nine purely short range quark-lepton level operators.

In this article, we have studied the \(\rho \) exchange contributions to \(0\nu 2\beta \,\)when the quarks in the short distance product of currents hadronize into \(\rho \) mesons, which are massive compared to the momentum flowing through the internal lines. This allows the \(\rho \) propagators to be shrunk to a point and employing VIA, the one and two \(\rho \) meson contributions to neutrinoless double beta decay amplitudes have been estimated. It is found that for the vector and axial-vector currents, as considered here, the \(\rho \) exchange results in a contribution which brings an almost complete cancellation when combined with the conventional short range left–left and right–right \(0\nu 2\beta \,\)amplitudes, while the mixed chirality left–right amplitude gets enhnaced by a factor \(\sim 3\). This marks the first step in highlighting the importance of including \(\rho \) exchange diagrams. A similar analysis can be carried for short distance operators with Dirac structure other than the vector/axial-vector considered in this paper. The two \(\rho \) exchange contribution due to tensor currents is quickly seen to lead to very similar results. Acknowledging the short comings of VIA, a better approach would be to map the quark-electron operators onto a set of operators in the chiral Lagrangian where the \(\rho \) is properly incorporated as a dynamical degree of freedom, interacting with the nucleons. This, like the pion-nucleon chiral EFT [114], will enable proper power counting once \(\rho \) meson is also included, such that possible hadronization both into pions and \(\rho \) is incorporated in a systematic fashion. Since rho-meson is much heavier and quite close to the chiral perturbation scale, one possibility to have a consistent framework would be similar to heavy hadron chiral perturbation theory (see [126,127,128]) and include the nucleons. More specifically, the vector mesons can be treated as static heavy fields with fixed four-velocity, \(v^{\mu },\, v^2=1\). Then, like for a heavy quark, these heavy fields can be written in terms of velocity dependent fields by factoring out the mechanical part: \(exp(-i\mu v.x)\), with \(\mu \) denoting a typical octect orr singlet mass scale. See [129] for some details on constructing chiral perturbation theory for vector mesons. The task for \(0\nu 2\beta \,\)case is more involved. Such an EFT can then be used for phenomenological analysis since these hadron level contributions (local here for the \(\rho \) while non-local for pion exchange) have significant impact. Beyond this step of lowest order evaluation within such a chiral EFT, lies systematic calculation of higher order contributions. As mentioned before, neutrinoless double beta decay has been heralded as a harbinger of new physics, and it is of utmost importance and relevance to have as precise predictions (for a given NME) for the half life to infer the underlying mechanism of lepton number violation.