# ρ exchange contribution to neutrinoless double beta decay

ρρ exchange contribution to neutrinoless double beta decay Namit Mahajan ∗ Theoretical Physics Division, Physical Research Laboratory, Navrangpura, Ahmedabad 380 009, India

We consider ρ meson contributions to neutrinoless double beta decay amplitude stemming fromthe hadronization of the short distance quark-electron currents. These contributions are evaluatedwithin vacuum dominance approximation. The one and two ρ exchange contributions aﬀect theFermi transition nuclear matrix element in a way that lead to near cancellations in the same chirality,left-left and right-right, short range amplitudes when these new contributions are combined withthe conventional short range amplitudes, while the left-right amplitude almost triples. This thennecessitates the inclusion of ρ exchange amplitudes in any phenomenological study, like in left-righttheories. Experiments have ﬁrmly established that the neu-trinos, which are massless within the Standard Model(SM) of particle physics, have non-zero, albeit tiny mass,and diﬀerent ﬂavours of neutrinos mix with each other(see [1] for the current best ﬁt values of the parame-ters). Further, neutrinos are electrically neutral, allow-ing them to be their own anti-particles ie can be Ma-jorana in nature [2]. Neutrinoless double beta decay(0 ν β ), ( A, Z ) → ( A, Z + 2) + 2 e − , provides an un-ambiguous way of establishing the Majorana nature ofthe neutrinos, and also the lepton number violation [3].Theoretically as well, 0 ν β decay is heralded as a use-ful probe of physics beyond SM, having particular rele-vance for neutrino masses and mass hierarchy. For an in-complete list discussing 0 ν β decay phenomenology andrelated signatuures see e.g. [4]. The search for neu-trinoless double beta decay thus constitutes an impor-tant endevour. Experimentally, studies have been car-ried out or planned on several nuclei ([5]- [27]). Onlyone of the experiments [5] (HM) has claimed observationof 0 ν β signal in G e . The half-life at 68% conﬁdencelevel is: T ν / ( G e ) = 2 . +0 . − . × y r . A combina-tion of the ﬁrst results from Kamland-Zen and EXO-200,both using X e , yielded a lower limit on the half-life T ν / ( X e ) > . × y r which is at variance with theHM claim, and so is the GERDA result.Neutrinoless double beta decay process can proceed viathe light neutrinos, ν 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 speciﬁc models like R-parityviolating supersymmetric theories or theories with lepto-quarks (see [28] and references therein for a quick reviewof some of the essential theoretical and experimental is-sues). The short range part due to the heavy physicsis due to intermediate particles with masses much largerthan the relevant scale of the process ∼ O (GeV), allow-ing for the heavier degrees of freedom to be systemati-cally integrated out, leaving behind a series of operatorsbuilt out of low energy ﬁelds, the up and down quarks ∗ Electronic address: [email protected] and electrons, weighted by the short distance coeﬃcients,called Wilson coeﬃcients (denoted by C A below). Thisprovides a very convenient and systematic framework toevaluate the decay amplitude in terms of short distancecoeﬃcients which encode all the information about thehigh energy physics. In this process of integrating outthe heavy degrees of freedom, the operators and thusthe eﬀective Lagrangian obtained is at the typical scaleof the heavy particles. Using then the renormalizationgroup equations (RGEs), perturbative QCD eﬀects canbe computed. These QCD corrections have been shownto be very signiﬁcant [29], in particular due to the colourmismatched operators (see also [30]). In this way, thehigh energy particle physics input gets separated fromthe low energy dynamics contained in the nuclear ma-trix elements (NMEs) of the quark level operators sand-wiched between the nucleon states (see [31] to get anidea of diﬀerent approaches to calculate these NMEs).These NMEs are usually the source of large uncertaintyto 0 ν β predictions, and in principle one could comparepredictions for various nuclei undergoing 0 ν β transitionto understand and eventually reduce the sensitivity onNMEs.To be concrete, we begin by considering the short rangequark-electron operators (denoted by L q e below in thetext) and as an example, consider a heavy right handedneutrino, mass M N , and SM gauge group. The resultingquark level 0 ν β amplitude takes the form A ∼ ( V ud T ei ) M W M N (cid:124) (cid:123)(cid:122) (cid:125) G ¯ uγ µ (1 − γ ) d ¯ uγ µ (1 − γ ) d (cid:124) (cid:123)(cid:122) (cid:125) J q,µ J µq ¯ e (1 + γ ) e c (cid:124) (cid:123)(cid:122) (cid:125) j e where V ud and T ei are the quark and leptonic mixingmatrix elements while G = G F /M N contains the shortdistance physics, with G F being the Fermi constant. Thephysical 0 ν β amplitude is obtained by sandwiching thequark level operators between the initial and ﬁnal nuclearstates, | i (cid:105) and | f (cid:105) , ﬁnally evaluated in terms of the NMEs: A ν β = (cid:104) f | i H e ff | i (cid:105) ∼ G (cid:104) f |J q,µ J µq | i (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) NME j e The short distance or high energy physics cleanly sep-arates from the low energy matrix elements. The low a r X i v : . [ h e p - ph ] F e b energy eﬀective Lagrangian is expressed as a sum ofoperators, O A weighted by the Wilson coeﬃcents C A : L e ff = G A C A O A , where we have allowed for more thanone G for more complicated theories. In the exam-ple considered above, there is only one operator O = J q,µ J µq j e = ¯ u i γ µ (1 − γ ) d i ¯ u j γ µ (1 − γ ) d j ¯ e (1+ γ ) e c ( i, j denoting the colour indices) and the corresponding Wil-son coeﬃcient C = 1. In other models like SUSY withR-parity violation or leptoquarks, Fierz transformationshave to be employed to bring the operators in form sim-ilar to above. The Lorentz and Dirac structure of thequark level operator involved decides which NME entersthe 0 ν β rate. Perturbative QCD corrections don’t justcorrect the Wilson coeﬃcients but also lead to colour mis-matched operators (typically with strength that is 1 /N c ( N c = 3 being the number of colours) of the colour al-lowed operators) which are then Fierz transformed andcan give diﬀerent operators and thereby bringing a hostof diﬀerent NMEs which would not have been expectedotherwise.In the usual treatment of the short range part (see[28]), one considers following simpliﬁcations: (i) assumethat the ﬁnal electrons are emitted in the S-wave ie thelong wavelength approximation is employed; (ii) closureapproximation ie the energy of the virtual neutrino ismuch larger than the nulcear excitation energy, therebyallowing to sum over the intermediate set of states withgreat ease; (iii) (non-relativistic) impulse approximationwhich allows to write the matrix elements of product ofthe quark currents between the nuclear states in terms ofsay the Fermi and Gammow-Teller matrix elements. Forthe case of quark currents being V-A form, the aboveprocedure implies that J q,µ ( (cid:126)x ) J µq ( (cid:126)x ) = (cid:88) n,m τ n + τ m + δ ( (cid:126)x − (cid:126)r n ) δ ( (cid:126)x − (cid:126)r m )( g V ( q ) − g A ( q ) (cid:126)σ n · (cid:126)σ m ) (1)which lead to the Fermi (vector part of the current) andGammow-Teller (axial-vector part of the current) nuclearmatrix elements, M F and M GT respectively M F = (cid:104) Ψ f | (cid:88) n,m H ( r n , r m , ¯ E ) τ n + τ m + | Ψ i (cid:105) (2) M GT = (cid:104) Ψ f | (cid:88) n,m H ( r n , r m , ¯ E ) τ n + τ m + (cid:126)σ n · (cid:126)σ m | Ψ i (cid:105) (3)In the above expressions, σ and τ are Pauli matrices and τ + = ( τ + iτ ) /

2; indices n, m run over all the nucleons inthe nucleus and H ( r n , r m , ¯ E ) denotes the heavy neutrinopotential. Further, g V (0) = 1 and g A (0) (cid:39) .

27. Ψ i,f arethe initial and ﬁnal state nuclear wave-functions. It isworth mentioning that in writing these matrix elementsstarting from the quark level currents, the dipole andanapole terms are not shown as they are much smallerthan those displayed above. For the present purpose, we shall choose to neglect them but they can be systemati-cally included. See [31] for details and updated numericalvalues of the NMEs calculated within diﬀerent schemes.The diﬀerences in the numerical values of NMEs are atthe heart of large uncertainties in the predictions.Considering only vector and axial quark currents, thefollowing set of short distance operators are to be con-sidered (showing only the colour matched operators; i, j in the operators below denote the colour indices) O LL = ¯ u i γ µ (1 − γ ) d i ¯ u j γ µ (1 − γ ) d j ¯ e (1 + γ ) e c O RR = ¯ u i γ µ (1 + γ ) d i ¯ u j γ µ (1 + γ ) d j ¯ e (1 + γ ) e c O LR = ¯ u i γ µ (1 − γ ) d i ¯ u j γ µ (1 + γ ) d j ¯ e (1 + γ ) e c O RL = ¯ u i γ µ (1 + γ ) d i ¯ u j γ µ (1 − γ ) d j ¯ e (1 + γ ) e c (4)with their associated Wilson coeﬃcients, 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 eﬀective ﬁeld theory (EFT) point of view, therewill be scalar-pseudoscalar and tensor operators (andtheir combinations) as well. For deﬁniteness, we focuson the operators in Eq.(4). It turns out that the ShortRange (SR) NMEs depend on whether the quark chiral-ities are same or diﬀerent: M LLSR = M RRSR = g V M F − g A M GT M LR + RLSR = g V M F + g A M GT (5)It is to be noted that M GT ∼ − (3-4) M F for all the nu-clei and diﬀerent determinations of NMEs considered (seefor example [31]). Below, for simplicity and convenience,we shall assume that the two terms on the right hand sideof the above equations diﬀer by a factor ∼

5, after takinginto account g V and g A . Since the discussion will notbe speciﬁc to a particular nucleus, we shall assume thisfor all the nuclei. Therefore, one has the approximateresults: M LLSR = M RRSR ∼ M F and M LR + RLSR ∼ − M F ,thereby yielding the following approximate forms for thecorresponding contributions to the amplitude: A LL/RRSR ∼ C LL/RR M F A LR/RLSR ∼ − C LR/RL M F (6)A given theory of lepton number violation implies aset of quark-lepton level operators comprising the eﬀec-tive Lagrangian. Using the RGEs and including eﬀects ofoperator mixing, the ﬁnal amplitude for the 0 ν β can bewritten in terms of short distance coeﬃcients and variousNMEs. The rate thus computed can be then contrastedwith the experimental limits on the half life of the neutri-noless double beta decay process for the speciﬁc nucleusand stringent limits are obtained on the parameters ofthe theory. Alternatively, an EFT point of view could beadopted and all the relevant operators are then written atthe hadronic scale, using which the amplitude and thusthe decay rate is computed, which then leads to con-straints on diﬀerent eﬀective coeﬃcients. One of thesetwo is followed in the phenomenological studies.The above is not the entire story. Till now, the0 ν β amplitude is calculated by directly evaluating thenuclear matrix elements from the quark currents. How-ever, it was shown in [32] that there are additional con-tributions wherein these quark currents ﬁrst hadroniseinto pions, and then one can evaluate the two pion ex-change contributions which were shown to be the dom-inant ones. These authors used the on-shell matchingconditions to match the matrix elements of quark-leptoncurrents between the pion states on to the hadron level ef-fective terms. Ref.[33] provided a systematic EFT set-upto match the quark-electron operators to a chiral eﬀec-tive theory with two pion and electron ( ππee ) vertices aswell as two nucleon, one pion and electron ( N N πee ) ver-tices and four nucleon and electron (

N N N N ee ) vertices,where π , N and e denote pion, nucleon and electron re-spectively. It also discussed inclusion of next to leadingand next to next to leading order contributions. It wasconﬁrmed that the two pion exchange contributions in-deed dominate though their calculation employed NaiveDimensional Analysis (NDA) for the power counting incontrast to the Vacuum Insertion Approximation (VIA)as employed in the on-shell matching in [32]. That VIAmay miss some fraction of the total contribution is notan unexpected satatement since in VIA, only the contri-bution due to vacuum is retained while that due to otherstates is neglected. However, VIA provides a simple andquick approximation to estimate these new eﬀects. Theone and two pion contributions are non-local as they in-volve pion propagators. Recall that m π + ∼

140 M eV and the typical momentum, (cid:126)q , ﬂowing through the prop-agators is O (200 M eV ). This approach has been furtherdeveloped and reﬁned in [34]. That the two pion contri-bution can play a rather important role in phenomenolog-ical studies has been re-emphasized in [35] in the contextof left-right symmetric theories. It is worth pointing outthat almost all phenomenological studies on 0 ν β do notinclude these contributions and as recently pointed outin [35], the inferences and limits drawn could be signiﬁ-cantly altered once these are included.Given that there are important contributions to0 ν β amplitude beyond those obtained by directly tak-ing the matrix elements of the short distance quark cur-rents between the nuclear states emanating from one andtwo pion exchange diagrams, one is led to ask if thereare additional contributions beyond the pionic ones. Inparticular, does ρ -meson exchange, analogous to pion ex-change, yield signiﬁcant contributions? We now focus ourattention on new contributions coming from hadroniza-tion of quark currents into ρ mesons. Relevant diagramsare shown in Fig.(1). Compared to the pion exchange,there are two diﬀerences: (i) m ρ (= 770 M eV ) >> | (cid:126)q | ( ∼

200 M eV ), (ii) m π << m N ( ∼

940 M eV ) ∼ . m ρ .Statement (i) means that it is a good approximationto neglect | q | in comparison with m ρ , both in the nu-merator and denominator of the ρ -meson propagator: − i ( g µν − q µ q ν /m ρ ) / ( q − m ρ ) → ig µν /m ρ ie at this levelof approximation, the ρ contribution is being treated as a local one. It is to be noted that this additional local con-tribution should be distinguished from with theconven-tional short range one, since having the explicit ρ -mesonspeciﬁc mass and couplings will play an important role,as we see below. FIG. 1: Hadronic level diagrams (drawn using JaxoDraw[39]) with the intermediate dotted lines denoting a ρ meson.(a) NNNNee ; (b) Two ρ exchange; (c) NNρee at one ver-tex and other vertex being ρNN (possibly including parityviolating terms) and the permuted diagram with the verticesexchanged.

The parity conserving strong ρN N vertices are de-scribed by the phenomenological Lagrangian (see for ex-ample [36]) L ρNN = g ρ ¯ N (cid:18) γ µ + i χ ρ m N σ µν q ν (cid:19) (cid:126)τ · (cid:126)ρ µ N (7)where again the nucleon anomalous magnetic momentterm, with strength χ ρ (also denoted as µ N in literature),is neglected as q << m N . There exist various determina-tions of g ρ and χ ρ but the combination g ρ (1 + χ ρ ) takesroughly the same value ≈

21. In what follows, we choose g ρ ∼ . L ρρee = G F m ρ m N ( m ρ a ρ (cid:126)ρ α · (cid:126)ρ α ) ¯ e (1 + γ ) e c (8)Using this vertex, along with the parity conserving ρN N vertex and approximating the ρ propagators as above,the amplitude for diagram in Fig. 1(b) takes the form A ρ = G F m ρ m N g ρ a ρ [ ¯ u p γ µ u n ][ ¯ u p γ µ u n ] ¯ e (1 + γ ) e c (9)To obtain a ρ , on-shell matching condition is employed: (cid:104) ρ, e |L qe | ρ (cid:105) = (cid:104) ρ, e |L ρρee | ρ (cid:105) (10)where to facilitate a quick comparison, the quark-electronLagrangian is rewritten as L qe = G F m N (cid:88) A C A O A ¯ e (1 + γ ) e c (11)with A = LL, RR, LR + RL . To complete the on-shellmatching, the last step needed is to make use of vacuum,dominance or VIA. This is achieved through the use offollowing deﬁning relations (see [38]): (cid:104) | ¯ uγ µ | ρ ( P ) (cid:105) = f ρ m ρ (cid:15) µ (cid:104) | ¯ uσ µν | ρ ( P ) (cid:105) = if Tρ ( (cid:15) µ P ν − (cid:15) ν P µ ) (12)with f ρ = 198 ± f Tρ = 160 ±

10 MeV. On-shellmatching, assuming VIA, yields a ρ = − f ρ m ρ (cid:88) A C A ∼ − . (cid:88) A C A (13)Next, consider the one ρ exchange diagram (c) inFig.(1) and its permuted one. The N N ρee interactionterm is written as L NNρee = G F m ρ m N a ρ ¯ N ( γ µ (cid:126)τ · (cid:126)ρ µ ) N ¯ e (1 + γ ) e c (14)resulting in one ρ exchange contribution to the0 ν β amplitude given by A ρ = − G F m N g ρ a ρ [ ¯ u p γ µ u n ][ ¯ u p γ µ u n ] ¯ e (1 + γ ) e c (15)Following the same steps as above for on-shell match-ing and factorizing the quark level product of currentsusing VIA as: (cid:104) p |J q,µ J µq | ρn (cid:105) = (cid:104) p |J q,µ | n (cid:105)(cid:104) |J µq | ρ (cid:105) , oneobtains after taking into account the combinatoric factors a ρ = 4 f ρ g ρ mρ g V (cid:88) A C A ∼ . (cid:88) A C A (16)where g V (0) = 1 has been used.The inclusion of the ρ exchange diagrams thus pro-duces an extra contribution A V IAρ to 0 ν β A V IAρ = A ρ + A ρ (17) ∼ − G F m N (7 (cid:88) A C A )[ ¯ u p γ µ u n ][ ¯ u p γ µ u n ] ¯ e (1 + γ ) e c One immediately notices from the Dirac structure of theamplitudes in Eq.(9) and Eq.(15) that ρ exchange onlyresults in Fermi transition matrix element. The net eﬀectdue to A ρ and A ρ is to modify the LL , RR and LR/RL short range amplitudes, Eq.(6): A LL/RRSR → C LL/RR M F − C LL/RR M F ∼ A LR/RLSR → − C LR/RL M F − C LR/RL M F = − C LR/RL M F (18) This constitutes the main result of this study, namelyan almost complete cancellation brought about due to ρ exchange contributions for the LL and RR short rangeamplitudes, A LLSR and A RRSR , and the mixed left-right am-plitude, A LR + LRSR getting almost tripled. These resultsclearly show that ρ exchange diagrams do impact quitesigniﬁcantly and have the potential to completely changethe phenomenological predictions, and thus call for thephenomenological analyses to be redone.In this article, we have studied the ρ exchange contri-butions to 0 ν β when the quarks in the short distanceproduct of currents hadronize into ρ mesons, which aremassive compared to the momentum ﬂowing throughthe internal lines. This allows the ρ propagators to beshrunk to a point and employing VIA, the one and two ρ meson contributions to neutrinoless double beta de-cay amplitudes have been estimated. It is found thatfor the vector and axial-vector currents, as consideredhere, the ρ exchange results in a contribution whichbrings an almost complete cancellation when combinedwith the conventional short range left-left and right-right0 ν β amplitudes, while the mixed chirality left-right am-plitude gets enhnaced by a factor ∼

3. This marks theﬁrst step in highlighting the importance of including ρ exchange diagrams. A similar analysis can be carried forshort distance operators with Dirac structure other thanthe vector/axial-vector considered in this paper. The two ρ exchange contribution due to tensor currents is quicklyseen to lead to very similar results. Acknowledging theshort comings of VIA, a better approach would be tomap the quark-electron operators onto a set of operatorsin the chiral Lagrangian where the ρ is properly incorpo-rated as a dynamical degree of freedom, interacting withthe nucleons. This, like the pion-nucleon chiral EFT [33],will enable proper power counting once ρ meson is also in-cluded, such that possible hadronization both into pionsand ρ is incorporated in a systematic fashion. Such a chi-ral EFT can then be used for phenomenological analysissince these hadron level contributions (local here for the ρ while non-local for pion exchange) have signiﬁcant im-pact. Beyond this step of lowest order evaluation withinsuch a chiral EFT, lies systematic calculation of higherorder contributions. As mentioned before, neutrinolessdouble beta decay has been heralded as a harbinger ofnew physics, and it is of utmost importance and rele-vance to have as precise predictions (for a given NME)for the half life to infer the underlying mechanism of lep-ton number violation. [1] I. Esteban, M. C. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni and T. Schwetz, JHEP , 106(2019); M. Tanabashi et al. [Particle Data Group],Phys. Rev. D , no.3, 030001 (2018) P. F. deSalas, D. V. Forero, S. Gariazzo, P. Mart´ınez-Mirav´e,O. Mena, C. A. Ternes, M. T´ortola and J. W. F. Valle,[arXiv:2006.11237 [hep-ph]].[2] E. Majorana, Nuovo Cim. , 171 (1937).[3] W. H. Furry, Phys. Rev. , 1184 (1939); J. Schechterand J. W. F. Valle, Phys. Rev. D , 2951 (1982).[4] W. -Y. Keung and G. Senjanovic, Phys. Rev. Lett. ,1427 (1983); R. N. Mohapatra, Phys. Rev. D , 3457(1986); M. Hirsch, H. V. Klapdor-Kleingrothaus andS. G. Kovalenko, Phys. Rev. Lett. , 17-20 (1995);M. Hirsch, H. V. Klapdor-Kleingrothaus and S. G. Ko-valenko, Phys. Rev. D , 1329-1348 (1996); M. Hirsch,H. V. Klapdor-Kleingrothaus and S. G. Kovalenko, Phys.Lett. B , 17 (1996); M. Hirsch, H. V. Klapdor-Kleingrothaus and S. G. Kovalenko, Phys. Rev. D ,1947-1961 (1998); M. Hirsch and J. W. F. Valle, Nucl.Phys. B , 60-78 (1999); A. Faessler, S. Kovalenkoand F. Simkovic, Phys. Rev. D , 055004 (1998);H. V. Klapdor-Kleingrothaus and U. Sarkar, Mod. Phys.Lett. A , 2469 (2001); H. V. Klapdor-Kleingrothausand U. Sarkar, Mod. Phys. Lett. A , 2243 (2003);V. Cirigliano, A. Kurylov, M. J. Ramsey-Musolf andP. Vogel, Phys. Rev. Lett. , 231802 (2004); G. Prezeau,Phys. Lett. B , 93-97 (2006); F. Deppisch andH. Pas, Phys. Rev. Lett. , 232501 (2007); B. C. Al-lanach, C. H. Kom and H. Pas, Phys. Rev. Lett. ,091801 (2009); M. Blennow, E. Fernandez-Martinez,J. Lopez-Pavon and J. Menendez, JHEP , 096(2010); A. Maiezza, M. Nemevsek, F. Nesti and G. Sen-janovic, Phys. Rev. D , 055022 (2010); A. Ibarra,E. Molinaro and S. T. Petcov, JHEP , 108 (2010);V. Tello, M. Nemevsek, F. Nesti, G. Senjanovic andF. Vissani, Phys. Rev. Lett. , 151801 (2011); M. Ne-mevsek, F. Nesti, G. Senjanovic and Y. Zhang, Phys.Rev. D , 115014 (2011); J. Bergstrom, A. Merle andT. Ohlsson, JHEP , 122 (2011); J. Chakrabortty,H. Z. Devi, S. Goswami and S. Patra, JHEP ,008 (2012); M. Mitra, G. Senjanovic and F. Vissani,Nucl. Phys. B , 26 (2012); M. Nemevsek, F. Nesti,G. Senjanovic and V. Tello, arXiv:1112.3061 [hep-ph];R. L. Awasthi, M. K. Parida and S. Patra, JHEP ,122 (2013); S. Pascoli, M. Mitra and S. Wong, Phys.Rev. D , no.9, 093005 (2014); J. Barry and W. Rode-johann, JHEP , 153 (2013); P. S. Bhupal Dev,S. Goswami, M. Mitra and W. Rodejohann, Phys. Rev.D , 091301 (2013); J. C. Helo, M. Hirsch, H. P¨asand S. G. Kovalenko, Phys. Rev. D , 073011 (2013);C. H. Lee, P. S. Bhupal Dev and R. N. Mohapatra,Phys. Rev. D , no.9, 093010 (2013); P. S. BhupalDev, S. Goswami and M. Mitra, Phys. Rev. D , no.11,113004 (2015); F. F. Deppisch, T. E. Gonzalo, S. Patra,N. Sahu and U. Sarkar, Phys. Rev. D , no.1, 015018(2015); S. Dell’Oro, S. Marcocci, M. Viel and F. Vis-sani, JCAP , 023 (2015); T. Peng, M. J. Ramsey-Musolf and P. Winslow, Phys. Rev. D , no.9, 093002(2016); G. Bambhaniya, P. S. B. Dev, S. Goswami andM. Mitra, JHEP , 046 (2016); M. L. Graesser, JHEP , 099 (2017); M. Duerr, M. Lindner and A. Merle,JHEP , 091 (2011); F. F. Deppisch, C. Hati, S. Patra,P. Pritimita and U. Sarkar, Phys. Rev. D , no.3, 035005(2018); H. Borgohain and M. K. Das, Int. J. Theor. Phys. , no.9, 2911-2934 (2017); D. Borah, A. Dasgupta andS. Patra, Int. J. Mod. Phys. A , no.35, 1850198 (2018);P. D. Bolton, F. F. Deppisch and P. S. Bhupal Dev, JHEP , 170 (2020); M. Frank, C. Majumdar, P. Poulose,S. Senapati and U. A. Yajnik, Phys. Rev. D , no.7,075020 (2020).[5] H. V. Klapdor-Kleingrothaus and I. V. Krivosheina,Mod. Phys. Lett. A , 1547 (2006).[6] A. Gando et al. [KamLAND-Zen Collaboration], Phys.Rev. Lett. , no. 6, 062502 (2013).[7] M. Auger et al. [EXO Collaboration], Phys. Rev. Lett. , 032505 (2012).[8] M. Agostini et al. [GERDA], Phys. Rev. Lett. , no.12,122503 (2013).[9] R. Arnold et al. [NEMO-3], Phys. Rev. D , no.11,111101 (2014).[10] J. B. Albert et al. [EXO-200], Nature , 229-234(2014).[11] K. Alfonso et al. [CUORE], Phys. Rev. Lett. , no.10,102502 (2015).[12] A. Gando et al. [KamLAND-Zen], Phys. Rev. Lett. ,no.8, 082503 (2016).[13] M. Agostini, M. Allardt, A. M. Bakalyarov, M. Balata,I. Barabanov, L. Baudis, C. Bauer, E. Bellotti, S. Bel-ogurov and S. T. Belyaev, et al. Nature , 47 (2017).[14] C. Alduino et al. [CUORE], Phys. Rev. Lett. , no.13,132501 (2018).[15] C. E. Aalseth et al. [Majorana], Phys. Rev. Lett. ,no.13, 132502 (2018).[16] M. Agostini et al. [GERDA], Phys. Rev. Lett. , no.13,132503 (2018).[17] O. Azzolini et al. [CUPID], Eur. Phys. J. C , no.11,888 (2018).[18] V. Alenkov, H. W. Bae, J. Beyer, R. S. Boiko,K. Boonin, O. Buzanov, N. Chanthima, M. K. Cheoun,D. M. Chernyak and J. S. Choe, et al. Eur. Phys. J. C , no.9, 791 (2019).[19] S. I. Alvis et al. [Majorana], Phys. Rev. C , no.2,025501 (2019).[20] O. Azzolini et al. [CUPID], Phys. Rev. Lett. , no.3,032501 (2019) doi:10.1103/PhysRevLett.123.032501[arXiv:1906.05001 [nucl-ex]].[21] D. Q. Adams et al. [CUORE], Phys. Rev. Lett. ,no.12, 122501 (2020).[22] O. Azzolini, J. W. Beeman, F. Bellini, M. Beretta, M. Bi-assoni, C. Broﬀerio, C. Bucci, S. Capelli, L. Cardani andE. Celi, et al. Eur. Phys. J. C , no.8, 702 (2020).[23] M. Agostini et al. [GERDA], Phys. Rev. Lett. ,252502 (2020).[24] E. Armengaud et al. [CUPID], [arXiv:2011.13243 [nucl-ex]].[25] J. Paton [SNO+], [arXiv:1904.01418 [hep-ex]].[26] E. Armengaud, C. Augier, A. S. Barabash, F. Bellini,G. Benato, A. Benˆoıt, M. Beretta, L. Berg´e, J. Billardand Y. A. Borovlev, et al. Eur. Phys. J. C , no.1, 44(2020).[27] J. B. Albert et al. [nEXO], Phys. Rev. C , no.6, 065503 (2018).[28] M. Doi, T. Kotani and E. Takasugi, Prog. Theor. Phys.Suppl. , 1 (1985); T. Tomoda, Rept. Prog. Phys. , 53 (1991); W. Rodejohann, Int. J. Mod. Phys. E , 1833 (2011); J. J. Gomez-Cadenas, J. Martin-Albo,M. Mezzetto, F. Monrabal and M. Sorel, Riv. Nuovo Cim. , 29 (2012); J. D. Vergados, H. Ejiri and F. Simkovic,Rept. Prog. Phys. , 106301 (2012); F. F. Deppisch,M. Hirsch and H. Pas, J. Phys. G , 124007 (2012);S. M. Bilenky and C. Giunti, Int. J. Mod. Phys. A ,no.04n05, 1530001 (2015); J. D. Vergados, H. Ejiri andF. ˇSimkovic, Int. J. Mod. Phys. E , no.11, 1630007(2016); C. Aalseth, H. Back, L. J. Dauwe, D. Dean,G. Drexlin, Y. Efremenko, H. Ejiri, S. Elliott, J. Engeland B. Fujikawa, et al. [arXiv:hep-ph/0412300 [hep-ph]];M. J. Dolinski, A. W. P. Poon and W. Rodejohann, Ann.Rev. Nucl. Part. Sci. , 219-251 (2019);[29] N. Mahajan, Phys. Rev. Lett. , 031804 (2014);N. Mahajan, Phys. Rev. D , no.3, 035015 (2014).[30] M. Gonz´alez, M. Hirsch and S. G. Kovalenko, Phys.Rev. D , no.1, 013017 (2016) [erratum: Phys. Rev.D , no.9, 099907 (2018)]; C. Arbel´aez, M. Gonz´alez,M. Hirsch and S. Kovalenko, Phys. Rev. D , no.9,096014 (2016) [erratum: Phys. Rev. D , no.9, 099904(2018)]; C. Arbel´aez, M. Gonz´alez, S. Kovalenko andM. Hirsch, Phys. Rev. D , no.1, 015010 (2017);M. Gonz´alez, M. Hirsch and S. Kovalenko, Phys. Rev.D , no.11, 115005 (2018); C. Ayala, G. Cvetic andL. Gonzalez, Phys. Rev. D , no.9, 094003 (2020).[31] F. Simkovic, A. Faessler, V. Rodin, P. Vogel and J. En-gel, Phys. Rev. C , 045503 (2008); A. Faessler,G. L. Fogli, E. Lisi, V. Rodin, A. M. Rotunno andF. Simkovic, Phys. Rev. D , 053001 (2009); J. Menen-dez, A. Poves, E. Caurier and F. Nowacki, Nucl.Phys. A , 139 (2009); P. K. Rath, R. Chandra,K. Chaturvedi, P. K. Raina and J. G. Hirsch, Phys.Rev. C , 064310 (2010); P. K. Rath, R. Chandra,P. K. Raina, K. Chaturvedi and J. G. Hirsch, Phys.Rev. C , 014308 (2012); A. Meroni, S. T. Petcovand F. Simkovic, JHEP , 025 (2013); J. Barea,J. Kotila and F. Iachello, Phys. Rev. C , 014315(2013); M. T. Mustonen and J. Engel, arXiv:1301.6997[nucl-th]; F. Simkovic, V. Rodin, A. Faessler and P. Vo-gel, Phys. Rev. C , 045501 (2013); R. Sahu andV. K. B. Kota, Int. J. Mod. Phys. E , no.03, 1550022(2015); J. Hyv¨arinen and J. Suhonen, Phys. Rev. C , no.2, 024613 (2015); J. Barea, J. Kotila and F. Iachello,Phys. Rev. C , no.1, 014315 (2013); J. Barea, J. Kotilaand F. Iachello, Phys. Rev. C , no.3, 034304 (2015);L. S. Song, J. M. Yao, P. Ring and J. Meng, Phys. Rev.C , no.2, 024305 (2017); C. F. Jiao, J. Engel andJ. D. Holt, Phys. Rev. C , no.5, 054310 (2017); S. Pas-tore, J. Carlson, V. Cirigliano, W. Dekens, E. Mereghettiand R. B. Wiringa, Phys. Rev. C , no.1, 014606 (2018);D. L. Fang, A. Faessler and F. Simkovic, Phys. Rev. C ,no.4, 045503 (2018); L. Graf, F. F. Deppisch, F. Iachelloand J. Kotila, Phys. Rev. D , no.9, 095023 (2018);F. F. Deppisch, L. Graf, F. Iachello and J. Kotila, Phys.Rev. D , no.9, 095016 (2020)[32] A. Faessler, S. Kovalenko, F. Simkovic andJ. Schwieger, Phys. Rev. Lett. , 183-186 (1997)doi:10.1103/PhysRevLett.78.183; A. Faessler, S. Ko-valenko and F. Simkovic, Phys. Rev. D , 115004(1998).[33] G. Prezeau, M. Ramsey-Musolf and P. Vogel, Phys. Rev.D , 034016 (2003).[34] V. Cirigliano, W. Dekens, J. de Vries, M. L. Graesserand E. Mereghetti, JHEP , 082 (2017); V. Cirigliano,W. Dekens, E. Mereghetti and A. Walker-Loud, Phys.Rev. C , no.6, 065501 (2018) [erratum: Phys. Rev.C , no.1, 019903 (2019)]; V. Cirigliano, W. Dekens,J. De Vries, M. L. Graesser, E. Mereghetti, S. Pas-tore and U. Van Kolck, Phys. Rev. Lett. , no.20,202001 (2018); V. Cirigliano, W. Dekens, J. de Vries,M. L. Graesser and E. Mereghetti, JHEP , 097 (2018);V. Cirigliano, W. Dekens, J. De Vries, M. L. Graesser,E. Mereghetti, S. Pastore, M. Piarulli, U. Van Kolck andR. B. Wiringa, Phys. Rev. C , no.5, 055504 (2019).[35] G. Li, M. Ramsey-Musolf and J. C. Vasquez,[arXiv:2009.01257 [hep-ph]].[36] E. Fischbach and D. Tadic, Phys. Rept. , 123-186(1973).[37] S. L. Zhu, C. M. Maekawa, B. R. Holstein, M. J. Ramsey-Musolf and U. van Kolck, Nucl. Phys. A , 435-498(2005).[38] P. Ball, V. M. Braun, Y. Koike and K. Tanaka, Nucl.Phys. B , 323-382 (1998)[39] D. Binosi and L. Theussl, Comput. Phys. Commun. ,76 (2004) [arXiv:hep-ph/0309015]; D. Binosi, J. Collins,C. Kaufhold and L. Theussl, Comput. Phys. Commun.180