Analysis of light neutrino exchange and short-range mechanisms in 0νββ decay

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Analysis of light neutrino exchange and short-range mechanisms in 0νββ decay © 2020 the Authors Published version Deppisch, Frank F.; Graf, Lukas; Iachello, Francesco; Kotila, Jenni Deppisch, F. F., Graf, L., Iachello, F., & Kotila, J. (2020). Analysis of light neutrino exchange and short-range mechanisms in 0νββ decay. Physical Review D, 102(9), Article 095016. https://doi.org/10.1103/physrevd.102.095016 2020 Analysis of light neutrino exchange and short-range mechanisms in 0νββ decay Frank F. Deppisch ,1,* Lukas Graf ,2,†Francesco Iachello,3,‡and Jenni Kotila4,3,§ 1Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, United Kingdom 2Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany 3Center for Theoretical Physics, Sloane Physics Laboratory, Yale University, New Haven, Connecticut 06520-8120, USA 4Finnish Institute for Educational Research, University of Jyväskylä, P.O. Box 35, FI-40014 Jyväskylä, Finland (Received 3 October 2020; accepted 20 October 2020; published 16 November 2020) Neutrinoless double beta decay (0νββ) is a crucial test for lepton number violation. Observation of this process would have fundamental implications for neutrino physics, theories beyond the Standard Model and cosmology. Focusing on so-called short-range operators of 0νββ and their potential interplay with the standard light Majorana neutrino exchange, we present the first complete calculation of the relevant nuclear matrix elements, performed within the interacting boson model (IBM-2). Furthermore, we calculate the relevant phase space factors using exact Dirac electron wave functions, taking into account the finite nuclear size and screening by the electron cloud. The obtained numerical results are presented together with up-to-date limits on the standard mass mechanism and effective 0νββ short-range operators in the interacting boson model framework. Finally, we interpret the limits in the particle physics scenarios incorporating heavy sterile neutrinos, left-right symmetry and R-parity violating supersymmetry. DOI: 10.1103/PhysRevD.102.095016 I. INTRODUCTION The nature of neutrinos and especially the origin of their masses are a crucial open question. While the Standard Model (SM) successfully explains the masses of the charged fermions it must be extended to incorporate neutrino masses. It would either require the presence of sterile neutrino states or effective lepton number violating (LNV) interactions. The first scenario allows the generation of Dirac neutrino masses analogous to those of the charged fermions. While certainly feasible, given the stringent upper limits mν≲Oð0.1ÞeV on the absolute neutrino masses from Tritium decay [1,2] and cosmological observations [3], tiny Higgs Yukawa couplings are required. Also, total lepton number Lwill no longer be an accidental symmetry. Unless Lsymmetry is imposed, the sterile neutrinos would acquire an LNV Majorana mass. The most popular example for such a scenario is the seesaw mechanism where the sterile neutrinos have such a large Majorana mass M≈1014 GeV naturally leading to light neutrino masses mν≈0.1eV [4–8]. High-scale seesaw mechanisms, or more generally scenarios where Lis broken at very high scales, are not the only way to generate light Majorana neutrino masses; other possibilities include LNV at low scales in secluded sectors, at a higher loop order, and when allowing higherdimensional effective interactions. If L-breaking occurs close to the electroweak scale, higher-dimensional LNV operators can be important. From a phenomenological point of view, searching for processes that violate total Lthus play a crucial role in neutrino and beyond-the-SM (BSM) physics. We here focus on the search for 0νββ decay as the most sensitive approach to probe Majorana neutrino masses. Currently, the most stringent limit on the 0νββ decay half-life T1=2is set in the Germanium isotope 76 32Ge [9], T1=2ð76GeÞ≡T1=2ð76 32Ge →76 34Se þe−e−Þ>1.8×1026 yr: ð1Þ However, Majorana neutrino masses are not the only contribution from BSM physics to 0νββ decay. We can generally consider the 0νββ decay rate by expressing high *[email protected] †[email protected]pg.de ‡[email protected] §[email protected] Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Funded by SCOAP3. PHYSICAL REVIEW D 102, 095016 (2020) 2470-0010=2020=102(9)=095016(27) 095016-1 Published by the American Physical Society scale new physics contributions in terms of effective lowenergy operators [10–13]. This only assumes that there are no exotic particles beyond the SM below the 0νββ energy scale of mF≈100 MeV. In this paper, we concentrate on so called short-range operators and their interplay with the standard light Majorana neutrino mass mechanism. As context, we provide a brief overview of the possible mechanisms for 0νββ decay which can be categorized in two main classes: (i) Long-range transitions via exchange of a light neutrino. This includes the so-called standard mass mechanism in Fig. 1(a), which is only possible if the neutrino is identical to its own antiparticle; i.e., if it is a Majorana fermion. The 0νββ decay rate can be estimated as Γ0νββ mν∼m2 νG4 Fm2 FQ5 ββ ∼ ðmν=0.1eVÞ2ð1026 yrÞ−1. Here, GFis the SM Fermi coupling and the phase space scales as Q5 ββ with the kinetic energy release Qββ ¼Oð1MeVÞfor typical double beta decays. Specifically, the mass mechanism of 0νββ decay is sensitive to the effective neutrino mass mββ ¼PiU2 eimνi, summing over the light Majorana neutrino masses mνiweighted by the square of the charged-current leptonic mixing matrix elements Uei. The inverse 0νββ decay half-life in a given isotope is then conventionally expressed as T−1 1=2¼jmββj2 m2 e GνjMνj2;ð2Þ with the phase space factor (PSF) Gνand the nuclear matrix element (NME) Mν. The normalization with respect to the electron mass meyields a small dimensionless parameter jϵνj¼jmββj=me. The current bound in Eq. (1) sets a limit jmββj≲ 79–180 meV at 90% confidence level (C.L.) for an unquenched axial coupling gA¼1.27 [9], with the uncertainty mainly due to the NMEs in different nuclear models. Future experiments will probe jmββj≈20 meV [14], corresponding to the minimal value for inversely ordered neutrinos. In BSM scenarios, a neutrino mass insertion is not necessarily required, cf. Fig. 1(b). In such cases, the decay rate is estimated as Γ0νββ LR ∼v2Λ−6 O7G2 Fm4 FQ5 ββ∼ ð105GeV=ΛO7Þ6×ð1026 yrÞ−1, with the SM Higgs vacuum expectation value (VEV) v¼246 GeV and the scale ΛO7of the exotic dim-7 operator. Such long-range mechanisms via the exchange of light Majorana neutrinos with interactions beyond the SM have received considerable attention [15–20], as the suppression at dim-7 is still fairly low and 0νββ decay is sensitive to high scales. We note, though, that due to the neutrino helicity-flip intrinsic in the operator, typical mechanisms are suppressed by the light neutrino masses. It is generically difficult to have a dim-7 operator where the exotic long-range contribution dominates over the standard mass mechanism [21], though it can be achieved in ultraviolet complete theories with a modestly suppressed standard contribution [22–25]. (ii) Short-range contributions where all mediating particles are heavier than mF≈100 MeV, cf. Fig. 1(c), are represented as contact interactions with six external fermions. These are the main focus of our analysis and they are generated by dim-9 and higher odd-dimensional operators. For a dim-9 operator, the decay rate can be estimated as Γ0νββ SR ∼ Λ−10 O9m6 FQ5 ββ ∼ð5TeV=ΛO9Þ10ð1026 yrÞ−1, with the operator scale ΛO9. The inverse 0νββ decay halflife triggered by such a mechanism is expressed similarly to Eq. (2) as T−1 1=2¼jϵIj2GIjMIj2, with the PSF GIand NME MI, both depending on the Lorentz structure of the effective operator. The coupling constant ϵIparametrizes the particle physics dynamics, i.e., the masses of the heavy states integrated out and their couplings. While such shortrange contributions do not involve the exchange of light neutrinos at all, they still require the breaking of lepton number and the SM neutrinos will be of Majorana type. The short-range and standard mass mechanisms are thus expected to compete but the relative strength is highly dependent on the underlying model. A detailed analytic derivation of the relevant NMEs for short-range operators was provided in our previous (a) (b) (c) FIG. 1. Contributions to 0νββ decay from effective LNVoperators. (a) Standard light neutrino exchange via dim-5 operator; (b) longrange contribution via dim-7 operator; (c) short-range contribution via dim-9 operator. Adapted from [26]. DEPPISCH, GRAF, IACHELLO, and KOTILA PHYS. REV. D 102, 095016 (2020) 095016-2 paper [27], where we included additional NMEs that become important when the latest values of the nucleon form factors are taken into account. Moreover, we calculated PSFs using the exact radial wave functions and we presented the single electron energy and angular correlation distributions for the exotic short-range 0νββ decay mechanisms. In the present paper, we numerically evaluate all relevant NMEs within the interacting boson model (IBM-2) framework. This will allow us to set upper limits on the effective couplings ϵIwhere we will highlight the exchange of heavy sterile neutrinos as an important example. Within the same framework, we also provide updated NMEs for the standard light neutrino exchange and we analyse its interference with short-range mechanisms. The NMEs for the 0νββ transitions are generally difficult to calculate and the limits derived are affected for any contribution. Detailed treatments using different nuclear structure model approaches can be found in [28–37]. Despite tremendous efforts to improve the nuclear theory calculation, the latest matrix elements obtained using various approaches differ in many cases by factors of ∼ð2–3Þ. The paper is organized as follows. We summarize the effective short-range Lagrangian at the quark level in Sec. II together with examples of underlying particle physics scenarios. The calculation of the 0νββ NMEs in the IBM-2 NME framework is outlined in Sec. III and that of the PSFs in Sec. IV. We then present our numerical results in Sec. Vwhere we provide up-to-date limits on the standard mass mechanism and effective short-range 0νββ operators. Section VI concludes our discussion with a summary and an outlook. II. SHORT-RANGE LNV OPERATORS AND NEUTRINO MASS MODELS In general, new physics where lepton number is broken at a a high scale will induce SM effective operators of dimension five, seven, nine, and higher [38,39]. After electroweak symmetry breaking, this will give rise to longand short-range contributions to 0νββ decay as outlined in the introduction, cf. Fig. 1. In this work we focus on shortrange contributions and their potential interplay with the standard mass mechanism. A. Effective Lagrangian The general effective short-range interaction Lagrangian can be written in terms of five different Lorentz-invariant classes of fermion current products [11], LSR ¼G2 Fcos2θC 2mpX C1;C2;cϵχ 1JC1JC2jcþϵχ 2Jμν C1JC2;μνjc þϵχ 3Jμ C1JC2;μjcþϵχ 4Jμ C1JC2;μνjνþϵχ 5Jμ C1JC2jμ þH:c:; ð3Þ where the sum is over all unique combinations χ¼ fC1;C 2ð;cÞg of chiralities C1;C 2;c¼R,Lof the quark and electron currents involved, JR;L ¼¯ uað1γ5Þda;J μ R;L ¼¯ uaγμð1γ5Þda; Jμν R;L ¼¯ uaσμνð1γ5Þda;ð4Þ jR;L ¼¯ eð1∓γ5Þec;j μ¼¯ eγμγ5ec:ð5Þ Here, the four-component Dirac spinor operators representing the up quark, down quark, and electron are denoted by u,d, and e, respectively. Quark SUð3ÞCcolor indices are denoted by a, and each quark current forms a color singlet in our parametrization. As the lepton current must violate lepton number by two units, the charge conjugate electron field ecappears there. Note that the chirality assignment in jR;L is flipped, i.e., the index Lis associated with 1þγ5. This is due to the appearance of the charge-conjugated electron field and, for example, the operator ¯ eð1þγ5Þec describes the creation of two left-handed electrons. Furthermore, the usual definition σμν ¼i 2½γμ;γνis used. The normalization of the Lagrangian by the factor G2 Fcos2θC=ð2mpÞwith the Fermi constant GF, the SM Cabibbo angle θCand the proton mass mpis conventional and results in dimensionless couplings ϵχ i. In principle, each unique current combination will be associated with a separate coupling, ϵχ i¼ϵC1C2ðcÞ i. Note that, in Ref. [11], the Lagrangian is defined without the factor cos2θC.We chose to include it as the resulting PSFs can be defined in the same way as that for standard light neutrino exchange, cf. Sec. IV C. Not all possible combinations of chiralities have to be considered in the Lagrangian Eq. (3), as redundancies and cancellations occur. First, the identity ½¯ uσμνð1þγ5Þd½¯ uσμνð1−γ5Þd ≡½¯ uσμνð1−γ5Þd½¯ uσμνð1þγ5Þd¼0ð6Þ implies that terms corresponding to ϵRLL 2,ϵLRL 2,ϵRLR 2, and ϵLRR 2trivially vanish. Second, the Pauli exclusion principle dictates that ¯ eγμec¼0and ¯ eσμνð1γ5Þec¼0, and thus any operator containing vector, tensor, or axial-tensor electron currents can be omitted. Altogether, the shortrange operators in Eq. (3) contain 24 independent ninedimensional operators invariant under the broken SM gauge group SUð3ÞC×Uð1ÞQ[27]. B. Example new physics scenarios with short-range contributions To illustrate the generation of different short-range contributions, we consider three well-known scenarios beyond the SM. ANALYSIS OF LIGHT NEUTRINO EXCHANGE AND SHORT- …PHYS. REV. D 102, 095016 (2020) 095016-3 1. Light and heavy neutrinos As discussed in the introduction, the exchange of light active Majorana neutrinos is the most prominent mechanism for 0νββ decay. As a long-range contribution, it is not represented in the Lagrangian Eq. (3) but arises from the SM charged current L¼GFcos θC ffiffiffi 2 p½¯ uγμð1−γ5ÞdX 3 i¼1 Uei½¯ eγμð1−γ5Þνi þH:c:: ð7Þ The sum is over the three SM neutrino mass eigenstates νi, constructed as the Majorana spinors νi¼νi;L þνc i;L from the SM active left-handed neutrinos νi;L and their charge conjugates. This gives rise to the mass mechanism of 0νββ decay sensitive to the effective Majorana neutrino mass mββ ¼X 3 i¼1 U2 eimνi:ð8Þ The 0νββ decay half-life in a given isotope is then conventionally expressed as in Eq. (2). One of the most attractive extensions of the SM involves adding fermionic states νi;S (i¼1;…;n N) that are sterile under the SM gauge interactions. They can thus acquire (Dirac or Majorana type) masses without spoiling the SM gauge invariance and eventually mix with the SM neutrinos after electroweak symmetry breaking. We can again form Majorana states by constructing Ni¼νi;S þνc i;S. The sterile states participate in the leptonic charged current due to mixing with the active neutrinos, L¼GFcos θC ffiffiffi 2 p½¯ uγμð1−γ5ÞdX nN i¼1 VeNi½¯ eγμð1−γ5ÞNi þH:c:; ð9Þ where VeNiare the elements of the active-sterile mixing matrix. If the sterile neutrinos are much lighter than the nuclear physics scale pF≈100 MeV, their contributions to 0νββ decay will be completely analogous to that of the active neutrinos and they can be included in Eq. (2) by replacing mββ →mββ þX nN i¼1 V2 eNimNi;ðmNi≪100 MeVÞ:ð10Þ Note that the Uei, and hence mββ, as well as the VeNiare in general complex numbers and cancellations can occur. In fact, if the Majorana states Niare solely responsible for the light neutrinos masses in a Seesaw scenario, the active and sterile contributions cancel to zero. If instead the sterile states are much heavier than the nuclear physics scale, mNi≫100 MeV, they can be integrated out, resulting in a contribution of the type Jμ LJL;μjLand the associated coupling ϵLLL 3is matched with the underlying physics parameters as ϵLLL 3¼X nN i¼1 V2 eNi mp mNi ;ðmNi≫100 MeVÞ:ð11Þ Note that the above considerations apply for sterile neutrinos that are Majorana fermions. This includes quasi-Dirac states that can be described by pairs of Majorana neutrinos (N1,N2) with a small mass splitting jmN1−mN2j≪mN1;2and a relative CP phase of π=2, VeN2¼iVeN1⇒V2 eN2¼−V2 eN1. In the limit of Dirac sterile neutrinos with mN1¼mN2, the contributions to 0νββ decay cancel. 2. Left-right symmetry The minimal left-right symmetric model (LRSM) is based on the extended gauge symmetry SUð3ÞC× SUð2ÞL×SUð2ÞR×Uð1ÞB−L[40–42]. It has a rich neutrino and 0νββ decay phenomenology as it naturally contains right-handed Majorana neutrinos Ni(i¼1,2,3)thatare charged under the SUð2ÞRpart of the gauge group, forming a doublet together with the right-handed leptons. This gives rise to right-handed charged currents, L¼g2 RcosθR C 8m2 WR½¯ uγμð1þγ5ÞdX 3 i¼1 UR ei½¯ eγμð1þγ5ÞNiþH:c:; ð12Þ mediated by a right-handed WRboson with the gauge coupling strength gRof the SUð2ÞRgroup. The angle θR C and the mixing matrix URare the right-handed equivalents of the Cabibbo angle and the Pontecorvo–Maki–Nakagawa– Sakata matrix, respectively. The LRSM gauge group is understood to be spontaneously broken to that of the SM at a high scale giving masses to the right-handed WRboson and neutrinos Ni. In turn, the active SM neutrino acquire masses via mixing with the heavy neutrinos (seesaw type I) as well as via the VEVof an electroweak triplet Higgs scalar present in the model (seesaw type II). Hence, the standard light neutrino and the sterile heavy neutrino contribution described above are generally present. In addition, the equivalent diagram with a heavy neutrino and two WRbosons contributes, giving rise to the short-range operator Jμ RJR;μjRwith jR¼¯ eð1−γ5Þec associated with ϵRRR 3matched to the underlying physics parameters as DEPPISCH, GRAF, IACHELLO, and KOTILA PHYS. REV. D 102, 095016 (2020) 095016-4 ϵRRR 3¼g2 R g2f2 LR X 3 i¼1ðUR eiÞ2mp mNi ;with fLR ¼gR g cos θR C cos θC m2 W m2 WR ;ð13Þ where gis the SM SUð2ÞLgauge coupling strength. Note that the contribution is not suppressed by the small, lightheavy neutrino mixing but instead by the expectedly high WRmass mWR. The right-handed mixing matrix URis approximately unitary with elements of order one, although cancellations due to complex phases can occur. The SM Wand the WRboson are also expected to mix with an angle as large as sin θW LR ≲gRm2 W=ðgm2 WRÞ. This permits the right-handed lepton current to couple with a left-handed quark current mediated by the SM Wgiving rise to the contributions ϵLRR 3¼ϵRLR 3¼sin θW LR fLR ϵRRR 3;ϵLLR 3¼sin2θW LR f2 LR ϵRRR 3: ð14Þ With the Wmixing taking the generic value sin θW LR ≈ gRm2 W=ðgm2 WRÞ≈fLR, all three effective couplings are of the same order. As mentioned, the LRSM also has the standard contribution from mββ and the sterile neutrino contribution ϵLLL 3in Eq. (11). In addition, the LRSM in principle also gives rise to the remaining contributions of type ϵ3, namely ϵLRL 3¼ϵRLL 3and ϵRRL 3but these are suppressed by both the light-heavy neutrino mixing and the high WRmass. Furthermore, the LRSM gives rise to additional long-range contributions that are not directly suppressed by the light neutrino masses. Finally, the LRSM has contributions from the electroweak triplet scalars ΔL;R that acquire VEVs vR,vL∼v2=vR during the spontaneous symmetry breaking, where vRis the breaking scale of the left-right symmetry. This gives rise to a diagram to 0νββ decay mediated by two WRbosons and the doubly charged scalars Δ−− L;R. Taking into account the W boson mixing, the contributions are ϵRRR 3¼g2 R g2f2 LR X 3 i¼1ðUR eiÞ2mpmNi m2 Δþþ R ; ϵLRR 3¼ϵRLR 3¼sin θW LR fLR ϵRRR 3; ϵLLR 3¼sin2θW LR f2 LR ϵRRR 3;ð15Þ analogous to Eqs. (13) and (14). Here, the heavy neutrino masses mNiappear because the couplings of the triplet Higgs to the gauge boson and electrons are proportional to vRand the heavy neutrino Yukawa coupling, mN∼yNvR. Likewise, there are contributions from the left-handed Δþþ L but they are additionally suppressed by the light neutrino masses (instead of mNi) and thus negligible. 3. R-parity violating supersymmetry As the final example of an ultraviolet-complete theory, we consider the minimal supersymmetric Standard Model with R-parity violation [43,44]. Without explicitly imposing invariance under the discrete Rsymmetry where each field carries the multiplicative quantum number R¼ ð−1Þ3BþLþ2S, with the baryon number B, total lepton number L, and spin S, the minimal supersymmetric Standard Model allows for the R-parity breaking terms, W⊃λijkLiLj¯ Ekþλ0 ijkLiQj¯ Dkþλ00 ijk ¯ Ui¯ Dj¯ Dk;ð16Þ in the superpotential. Here, the indices i,j,kdenote flavor generations of the superfields L,¯ E,Q,¯ D, and ¯ U, associated with the SM weak lepton doublet L, the lepton singlet ec, the quark doublet Qand the quark singlets dc,uc. Shortrange contributions to 0νββ are induced by the second term in Eq. (16), namely that associated with λ0 111 for the first lepton and quark generations [45]. They arise from diagrams with intermediate, heavy neutralinos, gluinos, squarks, and sleptons. The corresponding short-range Lagrangian is [46] LSR ⊃G2 Fcos2θC 2mpðϵRRL 1JRJRþϵRRL 2Jμν RJR;μνÞjL;ð17Þ i.e., a subset of the general short-range Lagrangian in Eq. (7) with scalar and tensor quark currents. The effective couplings ϵRRL 1and ϵRRL 2are generally functions of all supersymmetric particle masses and couplings involved. We here follow the assumptions of gluino dominance [46] where the diagrams involving gluinos and squarks contribute, ϵRRL 1¼8παsλ02 111 9cos2θC G−2 F m4 ˜ q mp m˜ g ;ϵRRL 2¼− 1 8ϵRRL 1:ð18Þ Here we also assume degeneracy of squark masses m˜q¼ m˜ uL¼m˜ dRin line with Ref. [46]. In addition, m˜ gis the gluino mass and αs¼0.127 is the strong fine structure constant at mW. Note that the gluino dominance assumption is based on the relevant NME values and limits on supersymmetry particle masses from other sources and may thus not be appropriate in light of new results. We nevertheless adopt it for simplicity and to compare with Ref. [46]. ANALYSIS OF LIGHT NEUTRINO EXCHANGE AND SHORT- …PHYS. REV. D 102, 095016 (2020) 095016-5 III. DETERMINATION OF NUCLEAR MATRIX ELEMENTS The NMEs for short-range mechanisms have been analytically derived in [27]. We follow the approach therein and summarize the basic formalism using nucleon form factors. A. Nucleon form factors The nucleon matrix elements of the colonr-singlet quark currents in Eq. (3) have the structure [47] hpj¯ uð1γ5Þdjni¼¯ Nτþ½FSðq2ÞFP0ðq2Þγ5N0;ð19Þ hpj¯ uγμð1γ5Þdjni¼¯ NτþFVðq2Þγμ−iFWðq2Þ 2mp σμνqνN0 ¯ NτþFAðq2Þγμγ5−FPðq2Þ 2mp γ5qμN0; ð20Þ hpj¯ uσμνð1γ5Þdjni¼¯ NτþJμν i 2ϵμνρσJρσN0;ð21Þ where τþdenotes the isospin-raising operator which converts a neutron into a proton, and the tensor Jμν in Eq. (21) is defined as Jμν ¼FT1ðq2Þσμν þiFT2ðq2Þ mpðγμqν−γνqμÞ þFT3ðq2Þ m2 pðσμρqρqν−σνρqρqμÞ:ð22Þ The above matrix elements generally depend on the neutron and proton momenta pn¼pN0and pp¼pN, respectively. The nucleon form factors are then functions of the momentum transfer q¼pp−pn. The most general parametrization of the vector current in Eq. (20) would include also induced scalar and axial-tensor terms—these can be, however, safely neglected, since they vanish in the isospin-symmetric limit and they are not enhanced by any other effects [48]. The momentum dependence in Eqs. (19)–(21) is encoded in the nucleon form factors FXðq2Þwith X¼S; P0;V;W;A;P;T1;T2;T3, usually parametrized in the so-called dipole form, FXðq2Þ¼gX=ð1þq2=m2 XÞ2. Here, the so called charge gXrepresents the value of the form factor at zero momentum transfer, gX≡FXð0Þ, and the scale mXdetermines the shape of the form factor. We apply this parametrization to all form factors except for the pseudoscalar form factors FP0ðq2Þand FPðq2Þ, which are enhanced by the pion resonance. The form factors with their corresponding parametrizations and charges are given by FSðq2Þ¼ gS ð1þq2=m2 VÞ2;g S¼1.0½49;ð23Þ FP0ðq2Þ¼ gP0 ð1þq2=m2 VÞ2 1 1þq2=m2 π ;g P0¼349 ½49; ð24Þ FVðq2Þ¼ gV ð1þq2=m2 VÞ2;g V¼1.0;ð25Þ FWðq2Þ¼ gW ð1þq2=m2 VÞ2;g W¼3.7;ð26Þ FAðq2Þ¼ gA ð1þq2=m2 AÞ2;g A¼1.269;ð27Þ FPðq2Þ¼ gP ð1þq2=m2 AÞ2 1 1þq2=m2 π ; gP¼4gA m2 p m2 π1−m2 π m2 A¼231 ½50;ð28Þ FTiðq2Þ¼ gTi ð1þq2=m2 VÞ2;g T1;2;3¼1.0;−3.3;1.34 ½47: ð29Þ The shape parameters are mV¼0.84 GeV, mA¼1.09 GeV [49] and the pion mass is mπ¼0.138 GeV. The form factors FVðq2Þ,FWðq2Þ, and FAðq2Þcan be determined experimentally and the parametrizations shown above provide a good description in the range 0≤jqj≤ 200 MeV of interest in 0νββ decay. On the other hand, as it is not possible to directly obtain the induced pseudoscalar form factor from experiment, we use the parametrization suggested in Ref. [50], which is based on the partially conserved axial-vector current hypothesis. The corresponding value of the free gPcharge agrees with the recent chiral perturbation theory analysis [51], which yields the value gP¼233. The value is also consistent with measurements of muon capture. With the muon mass mμ¼0.105 GeV, the resulting value of FPð−0.88m2 μÞ¼8.0 agrees well with the measured value of FPð−0.88m2 μÞ¼ 8.06 0.55 [52]. The scalar and pseudoscalar charges, gS and gP0, come from recent lattice QCD calculations [53].As there is not much information on the q2dependence of the corresponding form factors, we use the dipole parametrization, which, in the Breit frame, is the Fourier transform of the matter distribution. In the case of the pseudoscalar form factor we also include the monopole factor 1=ð1þq2=m2 πÞ used in chiral perturbation theory. As for the tensor form factors, only FT1enters our calculations. The value of the corresponding charge gT1quoted by Ref. [53] reads 0.987 0.055. We emphasize that the charges in Eqs. (23)–(29) are applicable at the free nucleon level. When calculating the 0νββ decay NMEs, we will use an DEPPISCH, GRAF, IACHELLO, and KOTILA PHYS. REV. D 102, 095016 (2020) 095016-6 effective axial-vector charge gA¼1.0and, consequently, an induced pseudoscalar charge gPðgA¼1.0Þ¼182 to approximately account for quenching in the nuclear medium. B. Nuclear matrix elements The five different types of quark current products appearing in Eq. (3) are mapped to the nucleon matrix elements according to Eqs. (19)–(21). By virtue of a nonrelativistic expansion and the closure approximation, the resulting product of nucleon matrix elements is then mapped to the nuclear matrix element between the final and initial 0þnuclear states involved in the 0νββ decay. This procedure is described in Ref. [27], and we here summarize the definition of NMEs involved. One should note that in the following expressions the relative sign between Gamow-Teller (GT) and tensor (T) terms is different than in our previous papers [28–30,54] and other available literature taking into account tensor terms using the formulation in [50]. The confusion about the relative sign arises from Eqs. (13) and (22) in [50], where in Eq. (13) a minus sign is used in front of the tensor term, while in Eq. (22) the plus sign is used. The tensor term contributes very little to the standard long-range mechanism, but, in the case of short-range mechanisms, it has a notable effect. Thus we have checked the derivation and concluded that the following signs should be used. The NMEs for the five short-range operators will generally depend on the chiralities of the two quark currents involved. For the first three operators associated with ϵχ 1,ϵχ 2, and ϵχ 3, the two quark currents are of the same type. Consequently, three possible combinations occur corresponding to the chiralities RR,LL,andðRL þLRÞ=2.It turns out that the resulting NMEs only depend on whether the quark chiralities are equal (RR,LL)ordifferent ðRL þLRÞ=2, represented by the upper and lower sign, respectively, in the expressions M1¼g2 SMFg2 P0 12 ðM0P0P0 GT þM0P0P0 TÞ;ð30Þ M2¼−2g2 T1 MT1T1 GT ;ð31Þ M3¼g2 VMFþðgVþgWÞ2 12 ð−2M0WW GT þM0WW TÞ ∓g2 AMAA GT −gAgP 6ðM0AP GT þM0AP TÞ þg2 P 48 ðM00PP GT þM00PP TÞ:ð32Þ For the operators associated with ϵχ 4and ϵχ 5, the two quark currents involved have different Lorentz structures and thus all four possible combinations of chiralities have to be considered in principle: RR,LL,RL,andLR. Again, it turns out that the NMEs only distinguish between the case where the quark chiralities are the same (RR,LL →upper sign) or different (RL,LR →lower sign), M4¼∓igAgT1MAT1 GT −gPgT1 12 ðM0PT1 GT þM0PT1 TÞ;ð33Þ M5¼gVgSMFgAgP0 12 ð ˜ MAP0 GT þ ˜ MAP0 TÞ −gPgP0 24 ðM0q0PP0 GT þM0q0PP0 TÞ:ð34Þ In the above expressions, we have explicitly factored the form factor charges gX¼FXð0Þ.Theqdependence arising from the product of the reduced form factors FXðq2Þ=gXis still to be included in the various matrix elements appearing in Eqs. (30)–(34).TheindividualFermi(MF), GamowTeller (MGT), and tensor (MT) NMEs along with the associated reduced form factor products ˜ hðq2Þare given in Table I. The numerical values of these NME will be given in Sec. III C, but we would like to note that the socalled recoil NMEs, ˜ MAP GT and ˜ MAP T, and the NMEs explicitly depending on the temporal momentum transfer q0,M0q0PP GT ,M0q0PP Tare difficult to evaluate exactly. We instead assume that the sum of nucleon spatial momenta is Q¼paþpb≈q[15–17], approximately applicable in an average sense considering that the NME is calculated summing over all nucleons in the nucleus. Similarly, we take the average value q0∼q2=mp≈10 MeV [17] for the temporal component of the momentum transfer. This allows us to reduce the corresponding NMEs as indicated in Table I. In addition to the product of the reduced nucleon form factors, the NMEs listed in Table Ialso contain the socalled neutrino potential describing the qdependence of the underlying particle physics mediator of 0νββ decay. Here we follow the formulation of [29,50] where the two-body transition operator is constructed in momentum space as the product of the neutrino potential vðqÞtimes the product of the reduced form factors ˜ hðq2Þ. In the case of the shortrange mechanisms we consider here, the neutrino potential is especially simple; as pointlike operators, they are described by a Dirac delta function in configuration space, δðra−rbÞ, hence in momentum space it is a q-independent constant. Following the usual normalization the short-range neutrino potential is [29,50] vðq2Þ¼2 π 1 memp :ð35Þ We also consider the standard light neutrino exchange mechanism with the NME ANALYSIS OF LIGHT NEUTRINO EXCHANGE AND SHORT- …PHYS. REV. D 102, 095016 (2020) 095016-7 Mν¼g2 VMF−g2 AMAA GT þgAgP 6ðM0AP GT þM0AP TÞ þðgVþgWÞ2 12 ð−2M0WW GT þM0WW TÞ −g2 P 48 ðM00PP GT þM00PP TÞ:ð36Þ Note that this is fully analogous to M3in Eq. (32) in the case where the quark currents have the same chirality, but the crucial difference is that the NMEs in Eq. (36) are calculated with the appropriate neutrino potential in momentum space [29], vðqÞ¼2 π 1 qðqþ˜ AÞ:ð37Þ Here, the neutrino mass has been neglected in comparison with the neutrino momentum q∼100 MeV, and ˜ Ais the closure energy, taken from Ref. [55] or estimated by the systematics, ˜ A¼1.12 ffiffiffiffi A pMeV. This describes the longrange exchange of an essentially massless neutrino mediating 0νββ decay in this case. As noted earlier, the relative sign between the GT and Tterms in Eq. (36) is different than in our previous papers [28–30,54] and other literature. Our derivation of the NMEs performed within the phenomenological framework of the nucleon form factors can be compared with an alternative way which has been developed in the literature over recent years. It is based on chiral effective field theory [56], i.e., the effective theory describing interactions at low energy in terms of baryons, mesons, photons, and leptons [20,24,57,58]. In this approach the process of hadronization is replaced by a perturbative expansion in terms of q=Λχreflecting the approximate chiral symmetry of QCD, where Λχ≃mp≈ 1GeV is the chiral symmetry breaking scale. The chiral Lagrangian on which the corresponding calculation is TABLE I. Double beta decay Fermi (MF), Gamow-Teller (MGT), and tensor (MT) NMEs appearing in Eqs. (30)–(34), with the associated reduced form factor product ˜ hðq2Þ. The NMEs are calculated using the functions h∘ðq2Þ¼vðq2Þ ˜ h∘ðq2Þenhanced by the neutrino potential Eq. (35) for short-range mechanisms and standard light neutrino exchange, Eq. (37). The subscript Xstands for X¼V;W;T1, for which the same form factor shape parameter mVapplies. The Pauli matrices in the space of the spins of the individual nucleons a,bare represented as σa;b and the tensor NMEs are calculated over Sab ¼3ðσa·qÞðσb·qÞ−ðσa·σbÞ. NME ˜ h∘ðq2Þ MF¼hhXXðq2Þi ˜ hXXðq2Þ¼ 1 ð1þq2=m2 VÞ4 M0P0P0 GT ¼hq2 m2 phPPðq2Þðσa·σbÞi ˜ hPPðq2Þ¼ 1 ð1þq2=m2 AÞ4 1 ð1þq2=m2 πÞ2 M0P0P0 T¼hq2 m2 phPPðq2ÞSabi ˜ hPPðq2Þ MT1T1 GT ¼hhXXðq2Þðσa·σbÞi ˜ hXXðq2Þ M0WW GT ¼hq2 m2 phXXðq2Þðσa·σbÞi ˜ hXXðq2Þ M0WW T¼hq2 m2 phXXðq2ÞSabi ˜ hXXðq2Þ MAA GT ¼hhAAðq2Þðσa·σbÞi ˜ hAAðq2Þ¼ 1 ð1þq2=m2 AÞ4 M0AP GT ¼hq2 m2 phAPðq2Þðσa·σbÞi ˜ hAPðq2Þ¼ 1 ð1þq2=m2 AÞ4 1 1þq2=m2 π M0AP T¼hq2 m2 phAPðq2ÞSabi ˜ hAPðq2Þ M00PP GT ¼hq4 m4 phPPðq2Þðσa·σbÞi ˜ hPPðq2Þ M00PP T¼hq4 m4 phPPðq2ÞSabi ˜ hPPðq2Þ MAT1 GT ¼hhAXðq2Þðσa·σbÞi ˜ hAXðq2Þ¼ 1 ð1þq2=m2 VÞ2 1 ð1þq2=m2 AÞ2 M0PT1 GT ¼hq2 m2 phXPðq2Þðσa·σbÞi ˜ hXPðq2Þ¼ 1 ð1þq2=m2 VÞ2 1 ð1þq2=m2 AÞ2 1 1þq2=m2 π M0PT1 T¼hq2 m2 phXPðq2ÞSabi ˜ hXPðq2Þ ˜ MAP0 GT ¼h Q·q m2 phAPðq2Þðσa·σbÞi≈M0AP GT ˜ hAPðq2Þ ˜ MAP0 T¼h Q·q m2 phAPðq2ÞSabi≈M0AP T ˜ hAPðq2Þ M0q0PP0 GT ¼h q0q2 m3 phPPðq2Þðσa·σbÞi≈10−2M0P0P0 GT ˜ hPPðq2Þ M0q0PP0 T¼h q0q2 m3 phPPðq2ÞSabi≈10−2M0P0P0 T ˜ hPPðq2Þ DEPPISCH, GRAF, IACHELLO, and KOTILA PHYS. REV. D 102, 095016 (2020) 095016-8 and their ratio K¼B=A, the angular distribution reads dΓ dcos θ¼Γ 2ð1þKcos θÞ:ð54Þ With the given information we determine the single electron distribution dΓ=dE1and the angular correlation αðE1Þfor the three relevant phase space factors that occur for short-range operators: fð0;1Þ 11þ(for mechanisms i¼1,2,3 with a scalar electron current), fð0;1Þ 66 (for mechanisms i¼4, 5 with an axial-vector electron current), fð0Þ 16 (for interference between the two classes), and fð0Þ 11−(for interference between i¼1, 2, 3 of different lepton chirality). As already noted, fð1Þ 16 vanishes, as does fð1Þ 11−. The electron phase space distributions fð0;1Þ 11þalso apply for the standard mass mechanism, calculated in the closure approximation. The resulting single energy distribution and angular correlation were already presented in Ref. [27] for several isotopes, but in Fig. 4(left) we illustrate the normalized single energy distributions for 76Ge as functions of the kinetic energy Ekin 1¼E1−meof one of the electrons; i.e., the range is from zero up to Qββ value. As can be seen, the term fð0Þ 11−produces an energy distribution virtually indistinguishable from that of fð0Þ 16 . All mechanisms produce a hill-like shaped energy distribution and observing the single energy spectrum is not expected to help distinguish between the standard mass mechanism (corresponding to fð0Þ 11þ) and any of the short-range mechanisms. The angular correlation αðEkin 1Þ, shown in Fig. 4(right), can distinguish between different mechanisms, namely short-range mechanisms of type i¼4, 5 produce electrons that are emitted collinearly whereas for i¼1, 2, 3 and the standard mass mechanism, they are dominantly back to back. As mentioned, the factors fð1Þ 16 and fð1Þ 11−vanish. There is therefore no change of the angular correlation due to interference and the angular correlation is an incoherent sum over contributions. C. Total decay rate Finally, we can integrate over the whole electron phase space to determine the total decay rate Γand the decay halflife T1=2, Γ¼ln 2 T1=2¼2CZQββþme me dE1wðE1ÞaðE1Þ:ð55Þ To facilitate calculation of the total rate under the presence of one or more mechanisms, we define the integrated PSFs [76] Gð0;1Þ ij ¼2C ln 2 gð0;1Þ ij 4R2 AZQββþme me dE1wðE1Þ ×fð0;1Þ ij ðE1;Q ββ þ2me−E1Þ;ð56Þ with gð0;1Þ 11¼1,gð0;1Þ 66 ¼1=16,gð0Þ 16 ¼1=4,gð1Þ 16 ¼0. The factor 1=R2 Ahas been introduced to conform with our convention where the NMEs are made dimensionless by multiplying with the nuclear radius RA. The numerical values of the PSFs Gð0;1Þ ij are given in Table VII, in units of 10−15 yr−1. As mentioned earlier, both Gð1Þ 16 and Gð1Þ 11− vanish, corresponding to the absence of interference in the angular part bðE1Þ. With the above PSFs, the inverse 0νββ decay half-life can be written 0.0 0.5 1.0 1.5 2.0 0.2 0.3 0.4 0.5 0.6 0.7 0.5 1.0 1.5 2.0 1.0 0.5 0.0 0.5 1.0 FIG. 4. Normalized single electron energy distributions Γ−1dΓ=dEkin 1(left) and angular correlation αðEkin 1Þ(right) for 76Ge as a function of the kinetic energy Ekin 1¼E1−me. Shown are the phase space factors for Eq. (50) in the former and for Eq. (51) in the latter. ANALYSIS OF LIGHT NEUTRINO EXCHANGE AND SHORT- …PHYS. REV. D 102, 095016 (2020) 095016-15 T−1 1=2¼Gð0Þ 11þX 3 I¼1 ϵL IMIþϵνMν 2 þGð0Þ 11þX 3 I¼1 ϵR IMI 2 þGð0Þ 66 X 5 I¼4 ϵIMI 2 þGð0Þ 11−×2ReX 3 I¼1 ϵL IMIþϵνMνX 3 I¼1 ϵR IMI þGð0Þ 16 ×2ReX 3 I¼1 ϵL IMI−X 3 I¼1 ϵR IMIþϵνMνX 5 I¼4 ϵIMI:ð57Þ Expressed in this way, the inverse half-life now only depends on the NMEs in Tables VI and IV (last column), the PSFs in Table VII and the coefficients ϵI,ϵν¼mββ=me encapsulating the particle physics aspects. V. RESULTS A. Bounds on the effective neutrino mass With ϵν¼mββ=meand the other short-range ϵIset to zero, Eq. (57) simplifies to the well know formula for light neutrino exchange, T−1 1=2¼jmββj2 m2 e Gð0Þ 11þjMνj2:ð58Þ Using the updated NME values for the light neutrino exchange mechanism shown in Table IV (last column) we can set new limits on the effective 0νββ mass jmββj.For isotopes with existing experimental bounds on the 0νββ decay rate, the resulting limits at 90% C.L. are summarized in Table VIII. As mentioned, the axial coupling is set to gA¼1.0. Generally, the limits have improved compared to the previous analysis [30]. This is a consequence of the better experimental limits for 76Ge, 82Se, 130Te, and 136Xe as TABLE VII. PSFs in units of 10−15 yr−1used in the calculation of the total decay rate for the standard light neutrino exchange and short-range mechanisms. The PSFs corresponding to fð1Þ 11−and fð1Þ 16 vanish. Isotope Gð0Þ 11þGð0Þ 11−Gð0Þ 66 Gð0Þ 16 Gð1Þ 11 Gð1Þ 66 [10−15 yr−1] 76Ge 2.360 −0.280 1.320 0.870 −1.954 0.977 82Se 10.19 −0.712 5.450 2.925 −9.079 4.539 96Zr 20.58 −1.190 10.88 5.403 −21.62 9.335 100Mo 15.91 −1.053 8.482 4.456 −14.25 7.125 110Pd 4.807 −0.541 2.674 1.730 −4.014 2.007 116Cd 16.69 −1.187 8.938 4.843 −19.37 7.414 124Sn 9.028 −0.843 4.935 2.976 −7.760 3.880 128Te 0.585 −0.156 0.371 0.313 −0.390 0.195 130Te 14.20 −1.142 7.672 4.367 −12.45 6.223 134Xe 0.597 −0.164 0.380 0.323 −0.394 0.197 136Xe 14.56 −1.197 7.876 4.524 −12.72 6.361 148Nd 10.07 −1.084 5.579 3.548 −14.19 4.246 150Nd 62.98 −3.125 33.05 15.44 −57.83 28.91 154Sm 3.005 −0.539 1.772 1.338 −2.291 1.145 160Gd 9.526 −1.129 5.321 3.506 −7.917 3.958 198Pt 7.513 −1.305 4.409 3.278 −5.844 2.922 232Th 13.87 −2.419 8.144 6.019 −10.92 5.457 238U 33.45 −4.176 18.81 12.46 −28.02 14.01 TABLE VIII. Upper limits on the effective 0νββ mass jmββjand the short-range ϵIcouplings in units of 10−10 from current experimental bounds Texp 1=2at 90% C.L., assuming a single contribution at a time and gA¼1.0. The chiralities of the involved quark currents are specified: the label XX stands for the case when both chiralities are the same, XX ¼RR; LL, and XY applies if the chiralities are different, XY ¼RL; LR. The limit on ϵ4applies for all chirality combinations. Isotope Texp 1=2[yr] jmββjjϵXX 1jjϵXY 1jjϵXX 2jjϵXX 3jjϵXY 3jjϵ4jjϵXX 5jjϵXY 5j [meV] ½10−10 76Ge 1.8×1026 [9] 118 2.90 2.84 88.4 77.1 154 130 102 68.1 82Se 2.4×1024 [77] 599 15.9 15.5 445 375 768 654 764 440 96Zr 9.2×1021 [78] 9130 85.5 84.8 5640 8510 12600 11300 1200 1110 100Mo 1.1×1024 [79] 733 6.10 6.04 401 608 901 774 84.1 77.5 116Cd 2.2×1023 [80] 2720 22.3 22.1 1430 2090 3170 2800 321 294 128Te 1.1×1023 [81] 13300 283 277 9300 8080 17300 12100 7630 5390 130Te 3.2×1025 [82] 252 5.38 5.27 178 153 336 270 158 112 136Xe 1.1×1026 [83] 114 2.50 2.45 83.4 72.5 157 127 74 52.4 150Nd 2.0×1022 [84] 3830 45.5 45.1 2730 3590 6190 5240 659 596 DEPPISCH, GRAF, IACHELLO, and KOTILA PHYS. REV. D 102, 095016 (2020) 095016-16 well as of the updated single particle energies in the NMEs for 76Ge, 82Se, 96Zr, and 150Nd. In Fig. 5, we compare the existing limit and future sensitivities in a plot correlating the 0νββ mass jmββjwith the sum of neutrino masses Σmν¼mν1þmν2þmν3for the standard picture of three active neutrinos. The shaded regions indicate, as usual, the allowed parameter space for normally (NO) and inversely (IO) ordered neutrino spectra by varying over the Majorana CP phases, where we take the best fit values of the oscillation angles and masssquared differences as given in [85]. Using our NMEs, the currently best limit is set by the KamLAND-Zen collaboration T1=2ð136XeÞ>1.1×1026 yr [83] resulting in jmββj<114 meV at 90% C.L. for gA¼1.0. The recent final result from GERDA with T1=2ð76GeÞ>1.8×1026 yr [9] corresponds to an essentially equal limit of jmββj< 118 meV at 90% C.L. In Fig. 5we also illustrate the corresponding limit assuming no quenching with gA¼ 1.27, giving jmββj<76 meV. In addition to the current limit we also show two examples of prospective sensitivities T1=2ð100MoÞ¼5×1026 yr expected at AMoRE-II [86] and T1=2ð76GeÞ¼1028 yr for LEGEND-1000 [87]. The latter will probe the full IO regime and a large chunk of the NO regime. Neutrino masses are also probed by the cosmological effect of the relic neutrino background on the cosmic microwave background and the structure of the universe. Current observations are compatible with no effect arising from neutrino masses setting stringent limits on Σmνdown to Σmν<150 meV at 90% C.L. [88] The limit generally depends on the neutrino ordering due to different priors in the statistical analysis and it is affected by the choice of the astrophysical data. It can also be weakened if an underlying cosmological model other than the standard minimal ΛCDM is used. In Fig. 5we show the most conservative limits arising from a choice of cosmological models surveyed in Ref. [88]. Namely, Σmν<280 meV (NO) and Σmν<290 meV (IO) at 95% C.L. arise in the ΛCDM with nonzero neutrino masses and a free scaling of the so-called weak lensing amplitude Alens (ΛCDMþ ΣmνþAlens). These limits correspond to jmββj<89 meV (NO) and jmββj<101 meV. B. Bounds on effective short-range mechanisms We can likewise assume that only a single short-range contribution is present by setting all other coefficients to zero and specifically assuming that the standard light neutrino contribution is negligible. Equation (57) then reduces to T−1 1=2¼jϵIj2GIjMIj2;ð59Þ with the PSF GIand NME MIdepending on the type of contribution. From the current nonobservation of 0νββ decay we can then set upper limits on the effective ϵI couplings. These are also shown in Table VIII, using our calculated PSFs and NMEs with gA¼1.0. Different chiralities of the quark currents in the operators lead to different bounds as indicated, where ϵXX idenotes the case where the chiralities of the two quark currents are equal, XX ¼ LL; RR, whereas ϵXY iindicates that they are different, XY ¼RL; LR.Forϵ2, the quark currents are required to be equal, cf. Eq. (6),andforϵ4, the bounds do not depend on the choice of the quark chiralities. Considering that a single ϵIcontributes at a time, the limits do not depend on the lepton chirality as the corresponding PSFs are independent of it. Numerically, the best limits for all ϵIare currently derived from the KamLAND-Zen constraint T1=2ð136XeÞ> 1.1×1026 yr, except for ϵXY 3where the GERDA constraint is slightly better. In any case, the KamLAND-Zen and GERDA bounds result in essentially equally stringent limits in most of the cases, and they are of the order 10−10 to 10−8.Forϵ1and ϵ5, in addition to the improved experimental bounds, the limits are the most stringent due to enhanced values of the nucleon current charges, specifically the large value of the intrinsic pseudoscalar 50 100 150 200 250 300 350 400 0 20 40 60 80 100 120 FIG. 5. Relation between the 0νββ mass jmββjand the sum of neutrino masses Σmνfor normally (NO) and inversely ordered (IO) neutrinos with the oscillation parameters fixed to the current best fit values. The dark shaded regions denote the parameter space allowed by the limits on Σmνat 95% C.L. from cosmological searches. The horizontal bars indicate the current upper limit on jmββjand future sensitivities of 0νββ decay searches with an unquenched (gA¼1.27, bottom edge) and quenched (gA¼1.0, top edge) value of the axial coupling. ANALYSIS OF LIGHT NEUTRINO EXCHANGE AND SHORT- …PHYS. REV. D 102, 095016 (2020) 095016-17 charge gP0, see Eq. (24). In case of ϵ3the sign of the tensor nuclear matrix elements also plays an important role. The limits in Table VIII on the effective couplings apply at the QCD scale ΛQCD ≈1GeV. As described in [27] following [89,90] one can instead define the couplings at the electroweak scale mW¼80.4GeV and evolve them to ΛQCD, where the appropriate bound can be set employing the experimental limit on the 0νββ decay half-life. Because different operators mix radiatively, a single contribution at mWmay generally induce several contributions at ΛQCD. The limits obtained in this way can be compared more directly with constraints derived from collider experiments. The resulting bounds on the couplings cI¼ϵIðmWÞat mW, including QCD running effects, are displayed in Table IX. Note that the limit on jϵ4jsplits into two different values jcXX 4jand jcXY 4j, since the different quark current chiralities affect the running. Numerically, the limits can be weaker or stronger than those at ΛQCD due to the overlapping effect of the QCD corrections and the mixing of operators. The already stringent limits on ϵXX 1and ϵXY 1improve further at mWand cXY 1is the most strongly constrained coupling by KamLAND-Zen. On the other hand, the limit on cXX 2is relatively much weaker than that on ϵXX 2. This is an effect of the renormalization group mixing with cXX 1and partial cancellation with this induced term. The effective short-range operator couplings can be interpreted in terms of effective New Physics operator scales ΛIwhere we simply match 1 Λ5 I¼G2 Fcos2θC 2mp cI;ð60Þ using the couplings cIdefined at the electroweak scale. In Fig. 6we illustrate the current bounds and expected future sensitivities in 76Ge (blue), 100Mo (orange), and 136Xe (green). The colored bars indicate the lower bound on the given operator scale where the darkest shade corresponds to the current limit and the two increasingly lighter shades represent expected future sensitivities. For the three isotopes, the setups are these: (i) T1=2ð76GeÞ=ð1026 yrÞ¼1.8 (GERDA [9], current), 10 (LEGEND-200 [87]), 100 (LEGEND-1000 [87]); (ii) T1=2ð100MoÞ=ð1026 yrÞ¼0.011 (NEMO-3 [79], current), 5 (AMoRE-II [86]), 10 (CUPID [91]); (iii) T1=2ð136XeÞ=ð1026 yrÞ¼1.1(KamLAND-Zen400 [83], current), 5 (KamLAND-Zen-800 [92]), 9.2 (nEXO [93]). As before, we assume only one short-range contribution to be present at mWand we neglect any contribution from light neutrino exchange. As can be seen, the strong limits on cXX 1and cXY 1probe operator scales up to 18 TeV. The weakest limits, applying to c2;3;4, still probe scales of order 4–6 TeV. C. Interference between light neutrino exchange and short-range mechanisms So far we have only considered one mechanism (operator) to be present at a given time, either the light neutrino exchange or one of the short-range operators. We now discuss the effect of two or more mechanisms operating at the same time. A large number of combinations are of course possible, but at least the standard light neutrino contribution is expected to be present at some level in any case. This is because any new physics scenario that generates a ΔL¼2short-range operator is also expected to generate Majorana neutrino masses at a level to explain neutrino oscillations. Therefore, it is reasonable to look into the interference of one of the nonstandard short-range mechanisms with the standard light neutrino exchange. We here discuss a few illustrative scenarios. We first consider the interference with the operator associated with ϵLLL 3. As we have seen in Sec. II B 1,it is triggered by heavy sterile neutrinos. Under the presence of ϵνand ϵLLL 3, Eq. (57) simplifies to T−1 1=2¼Gð0Þ 11þ mββ me MνþϵLLL 3MLL 3 2 :ð61Þ TABLE IX. As Table VIII, but for the short-range couplings cI¼ϵIðmWÞin units of 10−10, defined at the scale mW¼80.4GeV and omitting jmββj. Compared to Table VIII, the limits on c4depend on whether the quark currents have the same (XX) or opposite (XY) chirality. Isotope Texp 1=2[yr] jcXX 1jjcXY 1jjcXX 2jjcXX 3jjcXY 3jjcXX 4jjcXY 4jjcXX 5jjcXY 5j ½10−10 76Ge 1.8×1026 [9] 1.42 0.948 611 101 177 286 185 50.3 22.9 82Se 2.4×1024 [77] 7.74 5.19 2630 494 882 1450 934 361 148 96Zr 9.2×1021 [78] 42.9 28.5 26900 11200 14500 17300 16100 616 372 100Mo 1.1×1024 [79] 3.06 2.03 1930 800 1040 1200 1110 43.1 26.1 116Cd 2.2×1023 [80] 11.2 7.40 7390 2760 3650 4470 4000 165 98.9 128Te 1.1×1023 [81] 139 92.6 76800 10600 19900 26400 17300 3820 1810 130Te 3.2×1025 [82] 2.64 1.76 1490 202 387 589 386 79.2 37.5 136Xe 1.1×1026 [83] 1.23 0.819 717 95.4 180 277 181 37.2 17.6 150Nd 2.0×1022 [84] 22.8 15.1 18200 4720 7120 8720 7490 337 201 DEPPISCH, GRAF, IACHELLO, and KOTILA PHYS. REV. D 102, 095016 (2020) 095016-18 Because light neutrino exchange and the operator associated with ϵLLL 3have the same leptonic structure, the two contributions add coherently. The same behavior occurs for all operators of type ϵ1;2;3with a left-handed leptonic current. Depending on the complex phases of the NMEs and the particle physics parameters mββ,ϵLLL 3, the interference can be constructive or destructive. The NMEs are conventionally defined to be real with values given in Sec. III. In the given scenario, both NMEs are negative. We can choose mββ to be real and positive and the interference is described by the relative phase of ϵLLL 3. The largest effect then arises when ϵLLL 3is real and positive (constructive) or negative (destructive). Specifically, if ϵLLL 3¼−jmββj=meðMν=MLL 3Þ, both contributions cancel each other. The general constraints on the ðjmββj;ϵLLL 3Þparameter space are depicted in Fig. 7(left). As discussed, we take both mββ and ϵLLL 3to be relatively real with mββ >0by convention. The light shaded areas are allowed given the current limits from 0νββ decays searches in 76Ge and 136Xe, whereas the dark shaded area denotes the sensitivity from future searches at T1=2ð76GeÞ¼1028 yr. The combination of contributions in Eq. (61) leads to the linear relation between the variables and no independent limits can be set on them. From cosmological observations we can infer the upper limit jmββj<101 meV, see Fig. 5, and neither jmββj nor jϵLLL 3jcan be arbitrarily large given this additional constraint. Thus allowing a contribution 0≤jmββj< 101 meV from light neutrino exchange, ϵLLL 3is currently constrained to the interval −137 ×10−10 <ϵLLL 3< 72.5×10−10, compared to jϵLLL 3j<72.5×10−10 in the case it is the only contribution. In Fig. 7(right), we show the equivalent plot in the ðjmββj;ΛLLL 3Þparameter plane, where the effective operator scale is defined through 1=ðΛLLL 3Þ5¼ G2 Fcos2θCϵLLL 3=ð2mpÞ. The current experimental constraints give jΛLLL 3j≳4.5TeV. If mββ is not restricted further independently, e.g., by inference from an improved measurement of Σmν, the future constraint T1=2ð76GeÞ¼ 1028 yr will still allow ΛLLL 3≈−4.8TeV due to destructive interference. In the case of the interference between the standard light neutrino contribution with one operator of the type ϵ1;2;3with a right-handed lepton current or of the type ϵ4;5, the overlap is suppressed by the interference between the different lepton currents involved. We here discuss the example ϵRR 5in which case Eq. (57) simplifies to T−1 1=2¼Gð0Þ 11þjMνj2jmββj2 m2 eþGð0Þ 66 jMRR 5j2jϵRR 5j2 þ2Gð0Þ 16 ðMνMRR 5ÞRemββ me ϵRR 5; ¼Ajmββj2þBjϵRR 5j2−2CjmββjjϵRR 5jcosðα−βÞ: ð62Þ Here, A¼Gð0Þ 11þjMνj2=m2 e,B¼Gð0Þ 66 jMRR 5j2,C¼Gð0Þ 16 jMνj jMRR 5j=meare positive coefficients (we consider the NMEs to be real with Mν,MRR 5having opposite signs, cf. Table VI), and α,βare the complex phases of mββ, ϵRR 5, respectively. We again consider that the relative phase between mββ and ϵRR 5is α−β¼0;πin which case Eq. (62) is a quadratic function in jmββjand ϵRR 5and for a 2 4 6 8 10 12 14 16 18 20 TeV 1 XX 1 XY 2 XX 3 XX 3 XY 4 XX 4 XY 5 XX 5 XY 76Ge 100Mo 136Xe current future 0 FIG. 6. Lower limits on the effective short-range operator scales ΛIdefined at mWand assuming all other contributions are zero. The limits are from the current bounds (dark shade) and two future sensitivities (lighter shades) in 76Ge at ð1.8;10;100Þ×1026 yr (left, blue), 100Mo at ð0.011;5;10Þ×1026 yr (middle, orange) and 136Xe at ð1.1;5;9.2Þ×1026 yr (right, green). ANALYSIS OF LIGHT NEUTRINO EXCHANGE AND SHORT- …PHYS. REV. D 102, 095016 (2020) 095016-19 given value of T−1 1=2represent an ellipse. This is shown in Fig. 8(left) where the tilting is determined by the size of the PSF Gð0Þ 16 relative to Gð0Þ 11þand Gð0Þ 66 . The currently most stringent constraint is set in 136Xe but the limit on ϵRR 5from 100Mo is competitive despite the much weaker half-life limit. This is a consequence of enhanced NME MRR 5in 100Mo, see Table VI. Figure 8(right) shows the equivalent plot for the effective operator scale ΛRR 5. As can be seen in Table VII, the PSFs Gð0Þ 16 applicable to all contributions of type ϵ4;5are generally quite sizeable resulting in a comparatively strong interference. On the other hand, the PSF Gð0Þ 11−regulates the interference with operators of type ϵR 1;2;3with right-handed lepton currents, see Eq. (57), which is suppressed by the electron mass compared to the beta decay Qββ value. D. Constraints on new physics scenarios The above constraints on the effective neutrino mass and short-range operator couplings can be interpreted in terms of the new physics scenarios introduced in Sec. II B. 0 20 40 60 80 100 120 140 200 100 0 100 200 mmeV 3 LLL 10 10 m290 meV (IO) T 1 2 Ge 1.8 10 26 yr T 1 2 Xe 1.1 10 26 yr T 1 2 Ge 10 28 yr 0 20 40 60 80 100 120 140 10 5 0 5 10 mmeV 3 LLL TeV m290 meV (IO) T 1 2 Ge 1.8 10 26 yr T 1 2 Xe 1.1 10 26 yr T12 Ge 1028yr FIG. 7. Constraints on the effective 0νββ decay mass jmββjand the short-range operator coupling ϵLLL 3(left) as well as the associated operator scale ΛLLL 3(right). All other effective couplings are set to zero. The highlighted regions denote the allowed parameter space from the current limits T1=2ð76GeÞ>1.8×1026 yr (light blue) and T1=2ð136XeÞ>1.1×1026 yr (light green), as well as the future sensitivity T1=2ð76GeÞ¼1028 yr (dashed blue). The grey shaded area on the right is excluded assuming the bound on the sum of neutrino masses of Σmν>290 meV from cosmological observations for an inverse neutrino mass ordering. 0 20 40 60 80 100 120 140 100 50 0 50 100 m290 meV IO T 1 2 Mo 1.1 10 24 yr T 1 2 Xe 1.1 10 26 yr T 1 2 Ge 1.8 10 26 yr T12 Ge 1028yr T12 Mo 1027yr 0 20 40 60 80 100 120 140 10 5 0 5 10 m290 meV IO T 1 2 Mo 1.1 10 24 yr T 1 2 Xe 1.1 10 26 yr T12 Ge 1.8 1026yr T12 Ge 1028yr T12 Mo 1027yr FIG. 8. As Fig. 7, but for the short-range operator coupling ϵRR 5and associated scale ΛRR 5, also showing constraints from the current limit and future sensitivity in 100Mo. DEPPISCH, GRAF, IACHELLO, and KOTILA PHYS. REV. D 102, 095016 (2020) 095016-20 1. Light and heavy sterile neutrinos In the sterile neutrino case discussed in Sec. II B 1,we consider the simplified scenario where a single sterile neutrino of mass mNwith mixing VeN to the electron neutrino contributes to 0νββ decay. The limiting cases where the sterile neutrino is much lighter and heavier than 100 MeV were discussed in Sec. II B 1. Currently, the most stringent limit in Table VIII on 0νββ decay contributions of heavy sterile neutrinos is set in 136Xe, ϵLLL 3<72.5×10−10 ⇒X nN i¼1 V2 eNi mNi−1 >1.3×108GeV; ð63Þ assuming that the contributions from the light SM neutrinos are negligible. To approximately incorporate the intermediate range mN≈100 MeV as well, we use the interpolation [63,94] T−1 1=2¼Gð0Þ 11þjMLL 3j2mpmN hq2iþm2 N2 jVeNj4;ð64Þ with the average momentum transfer hq2i¼mpmejMLL 3= Mνj. In Fig. 9, we show the current limit and future sensitivity in the ðmN;jVeNj2Þparameter space. The region above the 0νββ bottom-most contours indicated are ruled out by the corresponding observation, assuming that the sterile neutrino is of a Majorana nature. We compare the 0νββ decay constraints with other searches for sterile neutrinos which are being pursued in neutrino oscillations, single beta decays, meson decays, at colliders and in electroweak precision measurements. The most recent searches are generally summarized in Ref. [95] and collider signatures are reviewed in Refs. [96,97]. The shaded area is excluded by current data and the dashed lines give examples of sensitivities in future searches. This includes the Tritium decay experiment KATRIN [98], searches for long-lived particles (LLP) (the shape is mainly determined by the planned DUNE [99],SHiP[100], and FCC-ee collider [101]) and high energy colliders FCC-hh [102], ILC [103], and CLIC [104,105]. As can be seen, future 0νββ decay searches at a level of T1=2≈1028 yr will be able to probe mixing strengths expected for light neutrino neutrino mass generation via the seesaw mechanism, mν¼ jVeNj2mN0.01 eV for mN≲100 MeV. Likewise, 0νββ decay searches probe very heavy Majorana neutrinos with masses up to mN≈106GeV where electroweak precision measurements can otherwise set comparatively weak limits of order jVeNj2≲10−3. We stress that this strong sensitivity of 0νββ decay searches applies to purely Majorana neutrinos, which are difficult to reconcile with the lightness of active neutrinos for mN≳1GeV. For sterile neutrinos with such masses it is more natural that they form quasi-Dirac states where LNV is suppressed by a small mass splitting. In Fig. 9, we also show the sensitivity towards such a scenario where two Majorana neutrinos with a relative mass splitting of ΔmN=mN¼10−4form a quasi-Dirac pair, partially cancelling their contributions to 0νββ decay. While the sensitivity is strongly reduced, 0νββ decay searches are 10−910−610−31 103106 10−10 10−8 10−6 10−4 10−2 1 mNGeV VeN 2 m V eN 2 m N 0.01 eV KATRIN LLP ILC CLIC FCC hh Current Direct Searches T 1 2 Ge 10 28 yr T 1 2 Xe 1.1 10 26 yr Quasi Dirac ( m N m N 10 4 ) Majorana FIG. 9. Upper limit on the active-sterile neutrino mixing strength jVeN j2as a function of the sterile neutrino mass mNfrom current 0νββ decay searches (solid curves) and at future sensitivities with T1=2¼1028 yr (dashed curves). The sterile neutrino is assumed to be of Majorana or quasi-Dirac nature as indicated and contributions from light neutrinos are neglected. The blue shaded area is excluded by current data from neutrino oscillations, beta decays, meson decays, colliders, and electroweak precision measurements. The dashed contours indicate the estimated future sensitivity in Tritium decays (KATRIN), long-lived particle (LLP) searches and at colliders (FCChh, ILC, CLIC). The diagonal line gives the seesaw relation of light neutrino mass generation, mν¼jVeN j2mN¼0.01 eV. ANALYSIS OF LIGHT NEUTRINO EXCHANGE AND SHORT- …PHYS. REV. D 102, 095016 (2020) 095016-21 still competitive at this level for mN≈1MeV and mN≈100 GeV. 2. Left-right symmetry In Fig. 10, we show the limits from 0νββ decay searches on the right-handed WRboson mass mWRand the heavy neutrino mass mNin the LRSM introduced in Sec. II B 2. Here, we consider a simplified scenario with one lepton generation, i.e., a single heavy neutrino Nand UR e1¼1.We also choose the so-called manifest left-right symmetric case with gR¼g,cosθR C¼cos θCand take mΔ−− R¼mWRfor the mass of the doubly charged triplet Higgs. The solid and dashed 0νββ curves give the lower limit on mWRwhere we additionally neglect the Wboson mixing, sin θW LR ¼0,thus ϵRRR 3is the only contribution. The rise of the 0νββ curves to the right of mN≈103GeV in Fig. 10 results from the doubly charged Higgs contribution in Eq. (15) increasing linearly with mN. Note, though, that too large values of mN, compared to mWR, are not natural as they would require nonperturbative Yukawa couplings with the triplet Higgs. The 0νββ decay limits are compared with the direct limits from the LHC and the future SHiP experiment. The current LHC limits arise from dijet, eþEmiss [106] and eejj signatures [107]. The future LHC limits are estimated for 300 fb−1of luminosity and are taken from [106]. The dijet and eþEmiss signatures are largely independent of the heavy neutrino mass in the applicable kinematic regimes and are sensitive to mWR≈4–7TeV. On the other hand, the SHiP experiment would probe heavy neutrinos produced mainly in Dmeson decays and the strong sensitivity around mN≈1GeV shown is taken from Ref. [108]. As can be seen, 0νββ decay searches are especially sensitive for mN≲20 GeV. Note that we only consider heavy neutrino masses as light as mN¼100 MeV where the short-range contribution assumption is reasonable. For mN≈100 MeV we incorporate the approximation in Eq. (64); masses around and below this scale can be incorporated using an analysis of the relevant dim-7 operators [18–20,109] and by including the mass dependence of the neutrino potential [63,94,110]. Both the LHC and SHiP limits were derived assuming negligible Wboson mixing; those based on the lifetime of the heavy neutrino will be affected and need to be reassessed. We nevertheless also include the sensitivity of future 0νββ decay searches for sin θW LR ¼m2 W=m2 WR, i.e., the generically maximal value expected, where all three operators ϵRRR 3,ϵLLR 3,ϵLRR 3contribute. Future searches are then expected to be sensitive up to mWR≈26 TeV. 3. R-parity violating supersymmetry Assuming gluino dominance, R-parity violating supersymmetry will induce the contributions in Eq. (18). Neglecting any other contributions, including those from light neutrinos, Eq. (57) simplifies to T−1 1=2¼Gð0Þ 11þð1.95MRR 1−2.88MRR 2Þ2 ×8παsλ02 111 9cos2θC G−2 F m4 ˜ q mp m˜g2 ;ð65Þ where the numerical factors in front of the NMEs take into account the effect of QCD running; i.e., we here interpret the coupling strength λ0 111 at mW. Using the current KamLAND-Zen bound T1=2ð136XeÞ>1.1×1026 yr, this can be translated into an upper limit on λ0 111, λ0 111 <7.0×10−3m˜ q 1TeV2m˜ g 1TeV1=2 :ð66Þ This compares to the limit λ0 111 <7.2×10−3in [111]1for the same squark and gluino masses and the above KamLAND-Zen bound. Somewhat surprisingly, the limits are thus of a very similar size; whereas, in our case, the strong sensitivity is predominantly due to the enhanced 0.1 1 10 102103104105 10 5 30 15 20 25 LHC (300 fb 1) LHC SHiP T 1 2 Xe 1.1 10 26 yr T 1 2 Ge 10 28 yr T12 Ge 1028yr (sin LR W) FIG. 10. Lower limit on the right-handed WRboson mass mWR as a function of the right-handed neutrino mass mNfrom current 0νββ decay searches (solid curve) and the corresponding future sensitivities with T1=2ð76GeÞ¼1028 yr (dashed curve) in the LRSM with negligible Wboson mixing. The dotted curve indicates the future sensitivity on the scenario where the Wboson mixing is sin θW LR ¼m2 W=m2 WR. The blue shaded area is excluded by current data from the LHC and the dashed contours indicate the estimated future sensitivity at the LHC with 300 fb−1and at SHiP. 1Reference [111] contains updated experimental constraints and includes the effects of QCD running to the scale 1 TeV compared to [46]. DEPPISCH, GRAF, IACHELLO, and KOTILA PHYS. REV. D 102, 095016 (2020) 095016-22 value of the NME MRR 1resulting from the large pseudoscalar form factor in Eq. (24), in Ref. [111] it is an effect of the QCD running and operator mixing. If 0νββ decay is not observed in future experiments with a sensitivity approaching T1=2ð100MoÞ¼1027 yr, the limit will improve to λ0 111 <2.0×10−3m˜q 1TeV2m˜ g 1TeV1=2 :ð67Þ This is mainly a result of the strong sensitivity to ϵRR 1 especially in 100Mo, see Sec. VB. As mentioned, the derived limit is based on the assumption of gluino dominance. It will be important to reevaluate the impact of 0νββ decay searches on the R-parity violating supersymmetry in light of the new results and the current constraints from direct searches for supersymmetric particles. VI. SUMMARY AND CONCLUSION Signatures of total lepton number violation are crucial if we want to understand the origin of neutrino masses, which constitute a key open issue in particle physics. Neutrinoless double beta decay has so far been the only practical means to probe light Majorana neutrino masses at scales indicated by neutrino oscillations. In addition it is sensitive to new physics contributions from exotic particles and interactions coupling to first-generation quarks and electrons. Within an EFT framework, 0νββ decay searches strongly constrain contributions of that form. In this work we have concentrated on short-range contributions which result from integrating out exotic particles much heavier than the energy scale mF≈100 MeV of double beta decay, leading to effective dimension-nine operators of the form Λ−5¯ u¯ udd¯ e¯ e. In addition, we update calculations for the standard light neutrino exchange mechanism to analyze the interplay with short-range contributions. We have presented a first complete numerical evaluation of the NMEs needed for the description of short-range nonstandard mechanisms of 0νββ decay. The calculation is performed within the framework of IBM-2 with restoration of the isospin properties of the Fermi transition operator. We also use updated single particle energies extracted from experimental data on nuclei with one nucleon removed or added from shell closure. We include additional NMEs that become important when the latest values of the nucleon form factors are taken into account. However, the main difference to previous calculations is in the sign of the tensor NMEs; the present derivation gives a sign of the tensor term MT, which is opposite to that in, e.g., [30]. This change has little effect on the standard mechanism, for which MTis small ≈1%, but it is sizeable for short-range mechanisms. As noted, we have performed our calculation in the phenomenological framework of the interacting boson model, using nucleon currents in the impulse approximation including higher-order terms in the nucleon momentum transfer determined in [27]. We model pion-mediated modes via enhanced pseudoscalar nucleon form factors informed by partially conserved axial-vector current and lattice QCD calculations. In our numerical results we consider a possible quenching of the axial-vector coupling by choosing gA¼1.0compared to the unquenched value gA¼1.27.We follow this classical approach in contrast to ab initio methods based on chiral EFT interactions [56]. Such formulations promise the determination of NMEs with controllable errors, e.g., may address part of the quenching problem [112]. Calculations of the standard light neutrino exchange NME Mνfollowing this approach have become possible for the lightest double beta decay isotopes 48Ca [113–115],76Ge, and 82Se [114], indicating noticeably smaller values than those from phenomenological models such as IBM-2, see [116] for a recent review. If confirmed, this will require an understanding for such a deviation as well further studies to apply ab initio methods to heavier nuclei. NMEs should ideally be verified experimentally by employing single and double charge exchange reactions [71,72].ChiralEFT techniques have been used to reveal a potentially sizeable short-range contribution in standard light neutrino exchange [57] and to calculate exotic contributions [20,24]. In addition to the NMEs calculated in our approach we also present the full set of leptonic PSFs for all relevant isotopes, determined numerically including effects from the finite nuclear size and electron cloud screening corrections. This allows us to set updated limits on the effective couplings of all possible short-range operators contributing to 0νββ decay. Considering one operator at a time, the current limits correspond to operator scales ranging between 3 to 10 TeV, where the strongest sensitivity is achieved for operators enhanced by pion-mediated corrections, in agreement with previous analyses [24,117–119], in our case arising from enhanced pseudoscalar form factors. We further illustrate the interplay between different contributions by considering the interference between the standard light neutrino exchange with one short-range contribution ϵIthus setting constraints on the combined parameter space ðmββ;ϵIÞ. Finally, we apply the effective operator framework to three example new physics scenarios, namely the SM with sterile neutrinos, left-right symmetry, and R-parity violating supersymmetry. Here, we set updated constraints on simplified parameter spaces and compare them with limits coming from other searches. Searches for lepton number violating signatures, with 0νββ decay as the most prominent example, are crucial for our understanding of neutrinos and physics beyond the SM in general. Given that no clear sign of new physics has been seen so far, short-range operators as those considered in this work provide a model-agnostic means to probe the presence of lepton number violating physics. Due to the strong suppression—the 0νββ decay rate scales as ∝Λ−10—future experimental advances increasing the sensitivity by up to ANALYSIS OF LIGHT NEUTRINO EXCHANGE AND SHORT- …PHYS. REV. D 102, 095016 (2020) 095016-23 two orders of magnitude to half-lives T0νββ 1=2≈1027−28 yr will only result in modest improvements in constraining Λ, see Fig. 6. Detailed analyses such as our work and [24] are still important as these operator scales Λ≈4–18 TeV are in a regime relevant for the LHC and potential future colliders. If an exotic short-range contribution were to be observed, it would indicate that light neutrino masses have their origin around the TeV scale. It would also have profound consequences on possible explanations of the matterantimatter asymmetry of the Universe, with the observation of nonstandard 0νββ decay contributions disfavoring baryogenesis mechanisms operating above the electroweak scale [26,120]. ACKNOWLEDGMENTS The authors would like to thank Jose Barea for providing the code to calculate standard mechanisms of double beta decay in IBM-2 and Patrick Bolton for sharing the direct sterile neutrino search limits and sensitivities. The authors would also like to thank Martin Hirsch for a careful reading of the manuscript and useful discussions. This work was supported in part by the U.S. Department of Energy (Grant No. DE-FG-02-91ER-40608) and the UK Royal Society International Exchange program. The work of J. K. was supported by the Academy of Finland (Grants No. 314733 and No. 320062). L. G. and F. F. D. acknowledge support from the UK Science and Technology Facilities Council (STFC) via a Consolidated Grant (Reference ST/P00072X/1). APPENDIX A: PARAMETERS OF THE IBM-2 HAMILTONIAN A detailed description of the IBM-2 Hamiltonian is given in [60,121]. For most nuclei, the Hamiltonian parameters are taken from the literature [122–135]. The new calculations are done using the program NPBOS [121]. TABLE X. Hamiltonian parameters employed in the IBM-2 calculation of the wave functions along with their references. Nucleus ϵdνϵdπκχ νχπξ1ξ2ξ3cð0Þ νcð2Þ νcð4Þ νcð0Þ πcð2Þ πcð4Þ πωνν ωππ ωνπ wνyν 76Ge [122] 1.20 1.20 −0.21 1.00 −1.20 −0.05 0.10 −0.05 76Se [123] 0.96 0.96 −0.16 0.50 −0.90 −0.10 82Se [123] 1.00 1.00 −0.28 1.14 −0.90 −0.10 82Kr [124] 1.15 1.15 −0.19 0.93 −1.13 −0.10 −0.10 96Zra1.00 1.00 −0.20 −2.20 0.65 0.17 0.17 0.33 96Mo [125] 0.73 1.10 −0.09 −1.20 0.40 −0.10 0.10 −0.10 −0.50 0.10 100Mo [125] 0.55 1.00 −0.06 −1.20 0.40 −0.10 0.10 −0.10 −0.60 0.20 0.10 100Ru [126] 0.89 0.89 −0.18 −1.00 0.40 0.60 0.09 −0.13 110Pd [127] 0.78 0.60 −0.13 0.00 −0.30 0.20 0.04 0.00 −0.26 −0.29 −0.30 −0.26 −0.29 −0.03 110Cd [128] 0.92 0.92 −0.15 −1.10 −0.80 1.10 0.109 1.10 0.07 −0.17 0.16 116Cd [129] 0.85 0.85 −0.27 −0.58 0.00 −0.18 0.24 −0.18 −0.15 −0.06 116Sn [130] 1.32 −0.50 −0.22 −0.07 −0.06 0.04 124Snb1.10 −0.30 −0.16 −0.20 0.30 0.02 124Te [129] 0.82 0.82 −0.15 0.00 −1.20 −0.18 0.24 −0.18 0.10 128Te [129] 0.93 0.93 −0.17 0.50 −1.20 −0.18 0.24 −0.18 0.30 0.22 128Xe [131] 0.70 0.70 −0.17 0.33 −0.80 −0.18 0.24 −0.18 0.30 130Te [129] 1.05 1.05 −0.20 0.90 −1.20 −0.18 0.24 −0.18 0.30 0.22 130Xe [131] 0.76 0.76 −0.19 0.50 −0.80 −0.18 0.24 −0.18 0.30 0.22 136Xeb1.31 −0.04 0.01 −0.02 136Ba [131] 1.03 1.03 −0.23 1.00 −0.90 −0.18 0.24 −0.18 0.30 0.10 148Nd [132] 0.70 0.70 −0.10 −0.80 −1.20 −0.12 0.24 0.90 0.40 0.20 148Sm [132] 0.95 0.95 −0.12 0.00 −1.30 −0.12 0.24 0.90 0.05 150Nd [132] 0.47 0.47 −0.07 −1.00 −1.20 −0.12 0.24 0.90 0.40 0.20 150Sm [132] 0.70 0.70 −0.08 −0.80 −1.30 −0.12 0.24 0.90 0.05 154Sm [132] 0.43 0.43 −0.08 −1.10 −1.30 −0.12 0.24 0.90 0.05 154Gd [132] 0.55 0.55 −0.08 −1.00 −1.00 −0.12 0.24 0.90 −0.20 −0.10 160Gd [135] 0.42 0.42 −0.05 −0.80 −1.00 0.08 0.08 0.08 −0.20 −0.10 160Dy [135] 0.44 0.44 −0.06 −0.80 −0.90 0.08 0.08 0.08 −0.05 −0.15 198Pt [133] 0.58 0.58 −0.18 1.05 −0.80 −0.10 0.08 −0.10 0.00 0.02 0.00 198Hg [134] 0.55 0.55 −0.21 1.00 −0.40 0.08 0.37 0.25 0.16 232Tha0.26 0.26 −0.05 −0.80 −1.45 0.20 0.20 0.20 232Ua0.28 0.28 −0.05 −1.00 −1.30 0.12 0.12 0.12 0.20 0.10 238Ua0.22 0.22 −0.05 −0.40 −1.30 0.12 0.12 0.12 0.20 0.10 238Pua0.24 0.24 −0.05 −0.60 −0.05 0.12 0.12 0.12 0.02 0.05 −0.09 aParameters fitted to reproduce the spectroscopic data of the low lying energy states. bGS parameters fitted to reproduce the spectroscopic data of the low lying energy states. DEPPISCH, GRAF, IACHELLO, and KOTILA PHYS. REV. D 102, 095016 (2020) 095016-24