Sterile neutrinos with non-standard interactions in $\beta$- and $0\nu\beta\beta$-decay experiments

Charged currents are probed in low-energy precision $\beta$-decay experiments and at high-energy colliders, both of which aim to measure or constrain signals of beyond-the-Standard-Model physics. In light of future $\beta$-decay and LHC measurements that will further explore these non-standard interactions, we investigate what neutrinoless double-$\beta$ decay ($0\nu\beta\beta$) experiments can tell us if a nonzero signal were to be found. Using a recently developed effective-field-theory framework, we consider the effects that interactions with right-handed neutrinos have on $0\nu\beta\beta$ and discuss the range of neutrino masses that current and future $0\nu\beta\beta$ measurements can probe, assuming neutrinos are Majorana particles. For non-standard interactions at the level suggested by recently observed hints in $\beta$ decays, we show that next-generation $0\nu\beta\beta$ experiments can determine the Dirac or Majorana nature of neutrinos, for sterile neutrino masses larger than $\mathcal O(10)$ eV.


Introduction
Investigations of β-decay processes have played an important role in the development of the Standard Model (SM) of particle physics. In the modern era, precision analyses of pion, neutron, and nuclear β-decay are used to extract precise values of SM parameters (in particular V ud G F ) and to probe beyond-the-Standard-Model (BSM) physics. At the relatively low energies associated with β-decay experiments, which involve Q-values around an MeV, BSM physics can be described using effective-field-theory (EFT) techniques as was already proposed by Lee and Yang [1], for recent reviews see e.g. Refs. [2][3][4][5][6][7]. Within this framework the V − A structure of the SM charged weak interaction is described by a dimension-six operator, while deviations are captured by additional higher-dimensional operators which are suppressed by powers of v/Λ, where v 246 GeV is the electroweak scale and Λ the scale of new physics. At the leading dimension-six level, the SM current is then augmented by ten effective interactions each of which comes with an, in principle complex, Wilson coefficient that scales as v 2 /Λ 2 .
These charged-current interactions can be investigated both at high energies, e.g. in pp → eν at the LHC, as well as in low-energy precision measurements. Focussing on probes in the latter category, Falkowski, González-Alonso, and Naviliat-Cuncic recently performed a state-of-the-art analysis of neutron and nuclear β-decay data [8] (we will refer to this work by FGN). In addition to neutron decay and superallowed 0 + → 0 + transitions, FGN included data from mirror 1/2 + → 1/2 + decays. The authors then performed a simultaneous fit including all leading-order dimension-six Wilson coefficients (considering only the real parts), while taking into account theory uncertainties by marginalizing over nuisance parameters related to radiative corrections and nuclear matrix elements (NMEs). Armed with this framework, FGN then analyzed all relevant data in several scenarios. Assuming only SM interactions FGN obtained accurate values for V ud = 0.97370 (25) and g A = 1.27276 (45), illustrating that β-decay experiments can be used to extract precise determinations of fundamental SM parameters and hadronic matrix elements. Adding non-standard interactions involving only left-handed neutrinos to the global fit gives an impressive confirmation of the SM: non-standard scalar and tensor interactions are constrained at the per-mille level, corresponding to scales of BSM physics of several TeV (note that pseudoscalar interactions give suppressed contributions and were not considered). Assuming the BSM scale is well above the electroweak scale, LHC measurements can be used to derive similar, complementary, constraints [2,3].
This picture changes somewhat once non-standard interactions involving right-handed neutrinos are included as well. In this case, the global fit constrains vector, axial-vector, and scalar couplings to right-handed neutrinos at the few percent level, but prefers a nonzero right-handed tensor coupling at the 10% level (about 3.2σ away from the SM point). This discrepancy is driven by a single recent measurement of the "little-a" coefficient in neutron decay [9] and is not significant enough to warrant too much excitement. In particular it is unclear whether one can construct a viable BSM scenario, given the current LHC constraints from pp → eν [2,3]. Nevertheless, in light of the FGN result, the recent hints for BSM effects in CKM unitarity tests [10][11][12][13], and more generally in context of ongoing and future β-decay experiments and collider measurements, it is interesting to ask what a signal of BSM charged-currents would imply for other complementary experiments. For instance, Refs. [2,3] demonstrated a strong complementarity between β-decay and LHC observables.
Here we focus on couplings to right-handed neutrinos and the connection to probes of lepton number violation (LNV). The presence of right-handed neutrinos in nature would imply the existence of fields that are neutral under the full SM gauge group 1 . As a result, nothing forbids the existence of a Majorana mass for fields like this unless additional symmetries are invoked. Such a mass term would imply that the mass eigenstates of neutrinos are Majorana and the violation of Lepton Number (L) by two units. So, if β-decay or collider experiments uncover evidence for right-handed neutrino interactions, a reasonable follow-up question would be: in which scenarios can we expect to find correlated signals in experimental probes of LNV? Here, we investigate this question in the context of searches for neutrinoless double beta decay (0νββ), the process where two neutrons in a nucleus transform into two protons with the emission of two electrons but zero (anti-)neutrinos.
To compare constraints from β-and 0νββ-decay experiments on couplings to right-handed neutrinos we need to consider the mass scale of right-handed neutrinos. β decays are insensitive to neutrino masses whenever the neutrinos are light, m i Q, while a 0νββ signal is proportional to the Majorana masses of the neutrinos. We therefore focus on right-handed neutrinos with masses well below the MeV scale that can be considered approximately massless for β decay experiments (as well as at colliders). Additional constraints on sterile neutrinos could in principle come from cosmological and astrophysical probes [17,18] such as big bang nucleosynthesis [19,20] considerations or the cosmic microwave background [21]. Such constraints however depend on the thermal history of the universe and can be avoided in specific scenarios [22].
In what follows, we calculate the expected 0νββ decay rates of various isotopes as a function of the non-standard couplings and the right-handed neutrino mass. Barring significant cancelations between the contributions from the 'standard mechanism' and those from sterile neutrinos, we then discuss the range of neutrino masses for which a β-decay or collider signal could imply a measurable signal in 0νββ. We start by setting up the EFT framework in Sect. 2 and briefly describe the FGN analysis in Sect. 3.1. The 0νββ decay rates in the presence of light righthanded neutrinos with non-standard interactions are discussed in Sect. 3.2. We present our analysis of EFT couplings in Sect. 4 and consider a specific BSM model involving leptoquarks in Sect. 5. We conclude in Sect. 6. Appendices are devoted to the matching relations with the neutrino-extended Standard Model Effective Field Theory and the calculation of 0νββ decay rates.

The Lagrangian in the Standard Model Effective Field Theory
In this work, we focus on interactions that give rise to β decays, especially those that involve right-handed neutrinos. We therefore consider an effective field theory that consists of the SM supplemented by an SM gauge singlet, ν R , and include higher-dimensional operators up to dimension six. The Lagrangian of the resulting EFT, called the neutrino-extended Standard Model EFT (νSMEFT) can then be written as where we used Ψ c = CΨ T for a field Ψ in terms of the charge conjugation matrix C = −C −1 = −C T = −C † . We use the definition for chiral fields Ψ c L,R = (Ψ L,R ) c = CΨ L,R T = P R,L Ψ c , with P R,L = (1 ± γ 5 )/2. Furthermore, L = (ν L , e L ) T is the left-handed lepton doublet andH = iτ 2 H * with H the Higgs doublet ( Here v = 246 GeV is the Higgs vacuum expectation value (vev), h(x) is the Higgs field, and U (x) is an SU (2) matrix encoding the Goldstone modes. Generally, ν R is a column vector of n righthanded sterile neutrinos, making Y ν a 3 × n matrix of Yukawa couplings andM R a symmetric complex n × n mass matrix. We will work in the basis where the charged leptons e i L,R and quarks u i L,R and d i R are mass eigenstates (i = 1, 2, 3). After electroweak symmetry breaking this then where V is the CKM matrix. We discuss the mass and flavor bases for the neutrinos in the next subsection.
Finally, the relevant dimension-five operators can be written as, while the needed dimension-six operators involving left-and right-handed neutrinos are collected in Tables 3 and 4 of appendix A, respectively.

The Lagrangian below the electroweak scale
Below the electroweak scale the Lagrangian of Eq. (1) induces the following mass terms for the neutrinos, where This matrix can be diagonalized by an N × N unitary matrix, U , This allows us to write the kinetic and mass terms of the neutrinos in the simple form in terms of the Majorana mass eigenstates ν = N m + N c m = ν c . These eigenstates are related to the flavor basis by where P and P s are 3 × N and n × N projector matrices which project onto the active and sterile states, respectively. We can use the above to write the charged-current interactions, which are induced by the dimension-six operators in Eq. (1) below the electroweak scale, in the mass basis of the neutrinos. In the notation of Ref. [23], this gives, at the scale µ 2 GeV, The matching of the above coefficients to those in Eq. (1) is described in App. A. In the SM, only the C VLL coefficient is nonzero. In the limit of vanishing neutrinos masses, the couplings C    Eq. (9) provides the interactions that are relevant for both single-and double-β decays. As such, the above Lagrangian has been studied extensively in the literature in the context of β decays, although usually using a different convention. We provide the translation to the notation of Ref. [8] in Table 2.1. Since β decays are low-energy observables, they are commonly described in terms of the interactions involving nucleons instead of the quark-level interactions given above. The Lagrangian, often used in studies of β decays, is then written in the following form 3 : In our analysis, we will focus on the couplings to neutrinos with right-handed chirality, corresponding to the couplings, C − V,A,S,P,T in the above Lagrangian, or C ΓAR of Eq. (9). Given the matching discussed in App. A, the relevant couplings are all induced by interactions that involve sterile neutrinos 4 .
From the above Lagrangian one can already see the form of the contributions to single-β decay and 0νββ. The contributions to β decay take the form where Γ αβ (m i ) depends on the Wilson coefficients, C α,β , under consideration and is nearly independent of m i for small neutrino masses, m i Q. In this case, the contributions through the SM charged current reduce to Γ SM and are independent of neutrinos masses and mixings. For large m i the partial rates Γ αβ (m i ) depend more sensitively on the neutrino masses until they vanish for m i > Q, as the decay to ν i becomes energetically disallowed. This effect can be searched for experimentally by looking for kinks in the electron spectrum [25,26].
Instead, the 0νββ decay rate is an LNV observable and therefore requires an explicit insertion of the neutrino mass leading to, 2 . Similar to β decay, the amplitudes A αβ (m i ) depend on the Wilson coefficients under consideration and are nearly independent of m i for small neutrino masses, m i m π . For large neutrino masses, the amplitudes scale as A αβ (m i ) ∼ 1/m 2 i , with a more complicated dependence in between, see Sect. 3.2 and Ref. [23] for more details. 3 The C ± α couplings of Ref. [8] are defined in the flavor basis and thus related to those used here by the change of basis in Eq. (7). In our numerical analyses, we will often focus on the case with n = 1 and U4α = δ4α for simplicity. In this scenario the C − α couplings used here coincide with those of FGN. 4 Note that these interactions can also be induced by operators involving left-handed neutrinos once one considers operators of dimension seven, since ν c L is right handed. However, such operators violate lepton number by two units and are already stringently constrained by 0νββ, leading to Λ > 10 TeV [24]. This implies their effects are unlikely to be measured in β experiments, given the sensitivity to the C − i couplings. Table 1: Translation to the notation of [8].
The relation between the parameters in Eq. (10) and those in Table 2.1 is given by [5] where g V,A,S,P,T are the vector, axial, scalar, pseudoscalar, and tensor charges of the nucleon [5], which can be determined using lattice calculations. Here, g V = 1 up to (negligible) quadratic corrections in isospin-symmetry breaking [27] and we use the FLAG'19 averages [28] for the axial, scalar and tensor charges: g A = 1.251(33), g S = 1.022(100) and g T = 0.989(33) [29,30] (see also Ref. [31]). Although the pseudoscalar charge is enhanced by the pion pole, namely g P = 349(9) [31], the suppression of the pseudoscalar contributions to single β-decay is larger such that β-decay experiments do not significantly constrain C − P . The matching in Eq. 11 includes the short-distance (inner) radiative corrections ∆ V R and ∆ A R . Especially the former is important, because it is necessary to extract V ud from nuclear data. Four recent calculations of this quantity are available [32][33][34][35], which all agree within 1σ. In this analysis we use the Seng et al. evaluation, ∆ V R = 0.02467(22) [32], which has the smallest uncertainty. Given the assumption about the reality of X and˜ X , all C ± X are then also real.

Single β-decay
As mentioned the interactions in Eq. (10) give tree-level corrections to the SM charged current which can be probed in β decays. We will follow the analysis of Ref. [8] to constrain the couplings in Eq. (10) 5 . Here we briefly summarize this analysis by FGN and discuss some of the details, related to theory uncertainties, of the fit. Since the observables considered by FGN include the β decays of the neutron as well as nuclei, the theoretical expressions depend of the NMEs in addition to the hadronic matrix elements appearing in Eq. (11). At leading order, nuclear effects are encapsulated in the so-called Fermi (F) and Gamow-Teller (GT) matrix elements, M F,GT . Leading-order expressions for beta decay observables in terms of the Lee-Yang Wilson coefficients C ± X can be found in Refs. [36,37]. However, given the experimental precision, subleading effects such as weak-magnetism and long-distance electromagnetic corrections have to be included in the SM terms. These small contributions can be calculated with high accuracy for the transitions that are included in this work [3,[38][39][40].
Apart from the corrections that can be accurately determined, there are several places where theoretical uncertainties are sizable. In particular, the "polluted" mixing parameters,ρ i , of [8] involve ratios of GT to F matrix elements which have not been calculated to sufficient precision. They are therefore treated as nuisance parameters and are marginalized over in the fit. Furthermore, theoretical uncertainties related to the nuclei-dependent parts of the radiative corrections is taken into account by including two additional nuisance parameters, η 2,3 [8].
In Ref. [8], the best-fit results are obtained by marginalizing over nuisance parameters and other Wilson coefficients. In this study, depending on the scenario under consideration, we do not always marginalize over all the Wilson coefficients. In some cases we keep certain Wilson coefficients fixed, in order to be able to compare directly with the 0νββ limits. As a result, the β limits we present are not always identical to those presented in Ref. [8]. We do always marginalize over all the nuisance parameters, even when only turning on one or two Wilson coefficients at a time.
Finally, although FGN neglected all neutrino masses, a complete analysis would in principle include the effects that the neutrino masses, m i , have on single-β observables. Instead of taking into account the full m i dependence of all β-decay observables in the fit, we will restrict ourselves to a range of m i for which the FGN analysis is still applicable. We estimate an upper bound on this range by considering the ways in which β observables will receive corrections due to nonzero m i . First, the contributions induced by the SM charged current are proportional to N i=1 |U ei | 2 = 1 for negligible neutrino masses, with corrections of the form 6 Here, the Q values of the isotopes that appear in the fit are in the MeV range, Q 0.7 MeV. Combining this with the upper limit on the mixing parameter U e4 , roughly |U e4 | 2 < 10 −3 for m 4 = 100 keV from the β-decay spectrum of 35 S [18,25,41], the corrections to the β decay rates due to a fourth neutrino are expected to be below the O(10 −4 ) level, for m 4 100 keV. In addition, the contributions of dimension-six operators will receive corrections ∼ Since the higher-dimensional operators will be suppressed compared to the SM, v 2 C − α 0.1, we expect these corrections to contribute at most at the O(10 −4 ) level. The observables in the FGN analysis, such as phase-space corrected lifetimes (Ft values) and β-decay asymmetries, have O(10 −3 ) uncertainties and are thus not sensitive to neutrino-mass corrections for m 4 100 keV.

Neutrinoless double beta decay
We now turn to the calculation of 0νββ decay rates of various nuclear isotopes. We closely follow the framework of Ref. [23]. Our starting point is the following effective Lagrangian of Eq. (9) This Lagrangian is written in the neutrino mass basis where ν denotes a 3 + n column vector of neutrino mass eigenstates. The Lagrangian is defined at a scale µ = 2 GeV and consists of operators that can be induced by operators of the dimension-six νSMEFT Lagrangian involving one ν R field. The only exception is the C VLL term which also includes the SM charged weak interaction of the active neutrinos.
We consider a simplified scenario where n = 1 implying a single sterile neutrino with arbitrary mass m 4 . While a pure 3 + 1 model (without dimension-five operators) would lead to two massless neutrinos, we assume some mass mechanism for the active neutrinos for instance through additional decoupled sterile neutrinos. As discussed above, we focus on m 4 masses well below the MeV scale when comparing to single-β experiments, where (in most cases) such neutrinos can be treated as massless 7 , while we consider larger masses when investigating the reach of 0νββ experiments. Under these assumptions it is possible to write down a relatively simple expression for the 0νββ decay rate in terms of electron phase-space factors, G 0i , and three so-called subamplitudes A L,R,M . The G 0i phase space factors are defined in Ref. [24] and have been calculated in the literature [45][46][47]. They include Coulomb corrections between the outgoing electrons and the nucleus and the (partial) screening effect of the nuclear charge by the surrounding electron cloud. Numerical values are given in Table 2.
The subamplitudes depend on the underlying LNV mechanism (in our case, the Majorana mass terms and the non-standard neutrino couplings) as well as hadronic and nuclear matrix elements. Each subamplitude is written as a sum over neutrino mass eigenstates Strictly speaking this description is only valid if m i < 1 GeV or so. Heavier neutrinos must be integrated out leading to effective dimension-nine interactions which give additional contributions in Eq. (14). In what follows we capture these effects by demanding that the m i -dependent nuclear [47] 76 and hadronic matrix elements have the correct m i behavior to reproduce the amplitudes that would result from integrating out the neutrinos at the quark level, in the m i 1 GeV region. As we focus on neutrino masses smaller than a GeV or so, we relegate a discussion of the procedure to App. B and refer to Ref. [23] for more details.
The first subamplitude is given by to be non-vanishing, C VLL requires an insertion of M R , see the last term in Eq. (22), so that the relevant terms C GeV. Such contributions lead to O(1) corrections to the amplitudes, but in most cases depend on QCD matrix elements that are not well known. For our numerical analysis we have used lattice QCD determinations [49] when possible, while the remaining matrix elements are estimated using naive dimensional analysis (NDA), see App. B and Ref. [23] for more details. Similarly while the final subamplitude contains only interference terms Here the A (ν) R,M (m i ) subamplitudes are again due to the exchange of hard neutrinos, which we discuss in App. B.
As mentioned, the M I (m i ) are linear combinations of different NMEs which are discussed in App. B. These elements in principle depend on the mass of the neutrino that is being exchanged. However, for neutrino masses below a few MeV, this mass dependence can be safely neglected. The mass-independent NMEs, M I (0), have been calculated with various nuclear methods. A collection of results is given in App. B, Table 5 for several isotopes. The differences between nuclear calculations lead to an additional O(1) uncertainty on the 0νββ decay results [50].
For heavier neutrino masses, the mass dependence of the NMEs and hadronic low-energy constants cannot be neglected. We follow the approach of Ref. [23] which interpolates between regions where m i k F ∼ m π where k F is the typical Fermi momentum in nuclei and m i > Λ χ ∼ 1 GeV where neutrinos can be integrated out at the quark level leading to local LNV operators. In the small-mass regime, NMEs and LECs are essentially mass independent and the amplitude scales linearly with m i . In the large-mass regime, the neutrino mass-dependence is captured by the Wilson coefficients of dimension-nine operators, ∼ (ūd) 2ē e c , that are generated after integrating out the neutrinos. The Wilson coefficients then have the form C (6) α C (6) β /m i so that the amplitude scales as 1/m i , while the NMEs and LECs are again mass independent. The mass dependence in the intermediate regime arises from two sources: 1) the neutrino propagator scales as m i /(q 2 + m 2 i ) instead of the Coulomb-like m i /|q|, 2) LECs associated to the exchange of hard virtual neutrinos pick up an intrinsic m i dependence. Here we will apply the formulae that were constructed in Ref. [23] by matching to results in the low-and high-mass regimes and interpolating in between. These formulae thus have the largest uncertainty in the m i ∼ Λ χ regime where the chiral and perturbative-QCD descriptions cannot be trusted, see App. B for explicit expressions.

Results
We now turn to the calculation of 0νββ constraints on the various effective interactions. We begin our analysis by considering a single C − X coupling, while setting the remaining Wilson coefficients to their SM values. We compare the constraints from 0νββ decay against single-β limits [8]. The Figure 1: 90% C.L. limits from the 0νββ decay of 136 Xe (red) compared to a fit of β decays (blue) in the scenario where only one of the v 2 C − i has been turned on. The width of the 0νββ bands indicate the variation induced by using shell-model [51] or QRPA [52] NMEs.
former are only valid for Majorana neutrinos. The latter are strictly speaking only valid in the limit of massless neutrinos but we will consider neutrino masses m ν R < 100 keV, where this is a reasonable approximation, see the discussion in Sect. 3.1.
To derive the contributions of the scenarios considered below to 0νββ we make several assumptions. We will consider a case with one additional neutrino, n = 1, and assume that the non-standard interactions only couple to the sterile state, i = 4 (this is equivalent to assuming there is no mixing between sterile and active neutrinos, U 4i ∝ δ 4i ). In addition, we will not consider the 'usual' mechanism for inducing 0νββ proportional to m ββ = m i U 2 ei . Both simplifications neglect additional contributions to 0νββ; while the neglected mixing effects can only contribute through interference terms and will tend to be suppressed (see Sect. 3.2), the ∼ m ββ terms can in principle be sizable. The reason we nevertheless make these simplifications is that these effects depend on the details of the neutrino masses and mixings, determined by M L,D,R in Eq. (4), which can only be assessed in a full model of neutrino masses. Barring cancelations, these additional contributions will tend to increase the predicted 0νββ rate, thereby strengthening the limits presented here. Our approach is thus conservative.

Single coupling analysis
Our results are shown in Figs. 1 and 2 for the 5 different couplings. The plots depict the m 4 − v 2 C − i plane, where the blue (red) shaded areas show the parameter space allowed by the single-β (double-β) constraints. To obtain the latter we assume there are no significant cancelations between the ν R contributions and those proportional to m ββ . The width of the 0νββ bands show the variation induced by using shell-model [51] or quasi-particle random phase approximation (QRPA) [52] NMEs, illustrating the impact of the nuclear theory uncertainties. The panels do not show the effect of changes in the LECs discussed in App. B.2, significant uncertainties of which would lead to additional errors at the O(1) level. For reference we also show the LHC limits and, for C − P , the constraint from Γ(π → eν)/Γ(π → µν) in green [2,30]. The collider limits were derived from pp → e + MET + X based on 20(36) fm −1 data recorded at GeV, with a more complicated mass dependence in between. The exact location of the peak sensitivity depends on the coupling under consideration because of non-trivial mass dependence of nuclear and QCD matrix elements. To get an idea of the sensitivity to BSM scales that are probed we can assume a typical scaling C − i = 1/Λ 2 . The constraints on C − V then imply a BSM scale between Λ > 0.7 TeV for m 4 = 10 eV to a maximum sensitivity of Λ > 47 TeV for m 4 = 0.5 GeV. These sensitivities should be compared to Λ > 1.5 TeV in the case of single-β decay. The dips that appear around m 4 ∼ GeV for C − A and C − S arise from the fact that the 0νββ decay rate vanishes. In these cases, for our particular choice of hadronic and nuclear matrix elements, the long-distance contributions cancel the short-distance terms. At this point, higher-order corrections should be included and our calculations are not reliable. Considering that this happens only for a tiny range of neutrino masses we do not pursue this further. Finally, we denote the prospective limit, T 0νββ 1/2 > 10 28 yr from EXO-200 [53] by the dashed red lines. As this limit is a factor ∼ 100 stronger than the current constraint, it leads to an improvement of ∼ √ 10 on the C − i constraints. The collider limits shown in Figs. 1 and 2 tend to be stronger than the single-β constraints by a factor of a few to an order of magnitude. However, the LHC limits from pp → eν only apply as long as the BSM physics responsible for the dimension-six interactions can still be treated as a contact interaction at high energies, i.e. as long as, Λ √ s ∼ TeV. Instead, the limits from single-and double-β decay are subject to a milder assumption, Λ GeV. In any case, if future LHC analyses find hints for interactions with light right-handed neutrinos, future 0νββ experiments should see a signal if the neutrinos are Majorana and have masses above roughly 1 keV.
Additional constraints on C − X arise from a careful analysis of the electron spectrum in neutrinoful double-β decay [54,55]. In particular, the angular distribution of outgoing electrons in 100 Mo double-β decay [56], .03 for m ν R < 0.1 MeV comparable to the single-β decay bounds. Future experiments are expected to improve this to the 0.01 level [54]. We have not indicated these bounds in the plot.
As seen from Fig. 2, single-β decay experiments are not sensitive to C − P because of a Q 2 /m 2 π suppression of the amplitude where Q denotes the Q-value of the experiment. However, the ratio Γ(π → eν)/Γ(π → µν) is rather sensitive to this coupling, leading to a stringent constraint that is overtaken by the 0νββ limit for m 4 30 keV.

Multi-coupling analysis
We now turn to a scenario where two couplings are turned on simultaneously, again setting the remaining couplings to their SM values. We focus on a comparison of single-β and 0νββ decay experiments. We show contours in the v 2 C − V -v 2 C − A plane in the left panel of Fig. 3, corresponding to a case in which right-handed neutrinos couple to quarks through (axial) vector currents. The single-β bounds are depicted in blue and indicate good agreement with the SM expectation, which is represented by the origin. The blue contour resides within a square with boundaries −0.03 < v 2 C − A,V < 0.03. In the case of Majorana neutrinos, the (future) 0νββ constraints are planes in the left and right panels, respectively. The blue (red) regions are those allowed by single-(double-)β decay constraints. In both cases we turned on only two Wilson coefficients at a time, setting the rest to their SM values.
shown by the (dashed) red contours, we again assumed there are no cancelations between the ∼ m ββ and non-standard contributions, while we used the NMEs of Ref. [52] and neglected the O(1) hadronic and nuclear uncertainties. The bounds span an oval that can roughly be described 0.05 × 120 eV m 4 . The current 0νββ constraints are thus competitive with single-β probes for m 4 > 200 eV or so. Next-generation experiments will push this to m 4 > 20 eV.
The right panel of Fig. 3 instead shows the v 2 C − A -v 2 C − T plane where the (not too significant) discrepancy found in Ref. [8] is visible. If this anomaly would be confirmed, the complementary nature of 0νββ could prove useful in identifying the nature and mass of right-handed neutrinos. With the current experimental 0νββ limits a nonzero value of C − A,T , consistent with the β-decay data, require neutrinos to be either Dirac or have masses below 300 eV. Next-generation 0νββ experiments will bring this down to 30 eV. Complementary constraints on neutrino masses in the eV range come from oscillation experiments [57,58], while searches for kinks in the electron spectrum of single-β decay experiments [18,25] are sensitive to masses up to the MeV range.
Finally, we use the same machinery to perform a fit of the β-decay data, where the single-β data is marginalized over all 5 couplings, while we only consider the C − A,T contributions to 0νββ. The resulting contour is shown in Fig. 4, which more clearly shows a preference for nonzero C − T . In this case, the preferred fit values of the Wilson coefficients are only consistent with Majorana neutrinos if the sterile neutrino has a mass below 30 eV. Next-generation experiments will push this 3 eV, covering essentially the entire possible mass range.

A leptoquark model
So far our analysis was performed purely in an EFT framework and we constrained effective couplings of the massive sterile neutrino. While this is convenient to demonstrate the complementarity between single-β and 0νββ experiments, it might also lead to oversimplifications.
To get a sense of more realistic BSM scenarios that induce multiple dimension-six operators at the same time, we consider a toy model with several leptoquarks (LQs) that have interactions with right-handed neutrinos. Integrating out the leptoquarks leads to dimension-six νSMEFT operators [6]. Following the notation of Ref. [59] we consider a BSM model with two scalar LQs, S 1 andR 2 . S 1 has mass M 1 and transforms under SU (3) c as an anti-triplet, under SU (2) L as a singlet, and carries nonzero hypercharge: S 1 (3, 1, 1/3), whileR 2 has mass M 2 and transforms under SU (3) c as a triplet, under SU (2) L as a doublet, and also carries nonzero hypercharge:R 2 (3, 2, 1/6). The relevant couplings to fermions can be written as where y LL,RL and y RR,LR are 3×3 and 3×n Yukawa matrices, respectively, while i, j are SU (2) L indices. Integrating out the two LQs at tree level gives rise to dim-6 operators and at low energies the effective interactions between neutrinos and the first-generation quarks and charged-leptons become where we neglected neutral-current operators that are not relevant for single-and double-β decays. Furthermore, α = 1, . . . 3 and a = 1, . . . n indicate the active and sterile neutrinos and the Wilson coefficients are given by In the mass basis the matching coefficients become A simplification is made possible by the fact that, as pointed out in Ref. [5], the Wilson coefficient C (6) VLL cannot be independently determined by the nuclear data alone. In particularly, in Eq. (11), it is not possible to distinguish the pure CKM element V ud from the new physics contamination parameterized by V ud (1 + L + R ). 9 We therefore absorb C (6) VLL into the CKM element V ud , and only consider two Wilson coefficients C (6) SRR and C (6) TRR in this simple LQ model. We show the resulting allowed regions in the M 1 −M 2 plane in Fig. 5, where we set y LR 11 y RL * 1e = y RR 11 y LL * 1e = 1 for simplicity. Here the red (dashed) contours show the (future) 0νββ limits for m 4 = 10 eV, while the blue regions denote the ∆χ 2 = 1, 2, 3 contours from single-β decays. These 0νββ limits were obtained assuming there are no significant cancelations between the ∼ m ββ and non-standard contributions. In addition, we used the NMEs from Ref. [52] and do not show the theoretical uncertainties due to the nuclear and hadronic matrix elements. The figure clearly shows the preference of the single-β fit for a nonzero tensor coupling, C − T , represented by the closed contours at ∆χ 2 = 1, 2, while for ∆χ 2 = 3 the fit allows for C − T → 0 and for the LQ masses to decouple. By assuming that neutrinos are Majorana one can again compare to the current 0νββ limits. From the figure, one can see that these constraints are already cutting into the region preferred by single-β decay at ∆χ 2 = 1, for m ν R = 10 eV. The prospective EXO-200 measurements would allow one the completely exclude the ∆χ 2 = 1 region and large parts of the ∆χ 2 = 2 parameter space.
It should be mentioned that the considered scenario does not provide a complete description of neutrino masses or attempt to reconcile with direct collider observables which typically exclude 9 The data only constrain four Wilson coefficients C + V,A,S,T , which depend on five parameters: V ud and L,R,S,T , leaving one flat direction. Figure 5: The red (dashed) line shows the current (future) limit from 0νββ for m ν R = 10 eV. The blue regions denote the ∆χ 2 = 1, 2, 3 contours from the single-β decay fit.
LQs with masses below 1−2 TeV [60,61]. A realistic model would therefore require more involved model building beyond the relatively simple extension of the SM considered here. We nevertheless consider the above scenario as a toy model as it clearly shows the complementary information on the nature and size of neutrino masses 0νββ could provide once a signal of nonzero BSM couplings would be found in probes of single β.

Conclusions
Sterile neutrinos play a role in numerous promising extensions of the SM and are strongly motivated by the necessity to generate neutrino masses. If sterile neutrinos exist, it is very possible that neutrinos are Majorana particles as no SM symmetries forbid sterile Majorana masses. In broad classes of SM extensions, sterile neutrinos only appear sterile at low energies, but interact through the exchange of heavy beyond-the-Standard-Model particles. Such SM extensions include left-right symmetric models, GUTs, leptoquark, and Z models. In this work we studied how such interactions can be probed in β-and 0νββ-decay experiments and in particular how these experiments complement each other in determining the mass and nature, Dirac or Majorana, of neutrinos.
To do so we employed the analysis of Ref. [8] which performed a comprehensive analysis of non-standard neutrino interactions in β-decay experiments. They demonstrated that precision β-decay experiments probe non-standard couplings to sterile neutrinos at the O(10 −1 − 10 −2 ) level, corresponding to BSM scales of a few TeV. They also found a slight hint for a nonzero tensor coupling to sterile neutrinos. We addressed what these results imply for the search of sterile neutrinos with 0νββ decay experiments. On the 0νββ side, we calculated the decay rates as a function of sterile neutrino masses and non-standard couplings at the dimension-six level of the neutrino-extended Standard Model Effective Field Theory. We focussed on rather light sterile neutrinos as only these can be produced in neutron and nuclear β-decay processes. It is worthwhile to mention that the KATRIN experiment, which measures the tritium β-decay spectrum, may also provide competitive constraints on the non-standard couplings of the sterile neutrino within our mass range. Given the current limits, the non-standard couplings could induce relative distortions of the spectrum at the per mille level [62]. As was shown for a benchmark calculation with m ν R ∼ 5 keV, these effects could allow for improved constraints assuming that the whole electron energy spectrum of KATRIN is accessible and theoretical and systematical uncertainties can be reduced below the per mille level.
For sufficiently small masses, β-decay constraints are neutrino-mass independent whereas 0νββ decay rates scale quadratically with the Majorana mass. We showed that, barring cancelations between ∼ m ββ and non-standard contributions, 0νββ measurements start to become competitive for right-handed neutrino masses larger than about 100 eV with a mild dependence on the particular interaction under consideration. Under the same assumptions, we also discussed the range of neutrino masses for which one would expect a measurable signal in 0νββ, assuming couplings of right-handed neutrinos were to be confirmed in either single β decays or at the LHC. For right-handed neutrino couplings of the size indicated by the hint in the recent β-decay results, we show that next-generation 0νββ experiments would be able to determine the Dirac or Majorana nature of neutrinos for m ν R 3 eV. Finally, we discussed the implications of the hints in single-β and the 0νββ limits in the context of a specific (toy) model of BSM physics involving two leptoquarks. Such a model with TeV leptoquark masses can account for the discrepancy and imply a measurable signal in next-generation 0νββ experiments unless neutrinos are Dirac particles.
Our study demonstrates the complementary nature of different neutrino probes, and the various aspects of the theory space they individually are sensitive to. In case of a future discovery of a deviation from the Standard Model predictions in one of the experiments, looking at other probes in a global approach provides valuable information to better understand the nature of neutrinos.
Hud P U, C where the first contribution to C VLL is the SM contribution. The Wilson coefficient of the operators involving right-handed chiralities are, B Details of the 0νββ calculation Here we briefly discuss the ingredients involved in the calculation of 0νββ decay rates that were not fully explained in Sect. 3.2. We start with the definition of the nuclear matrix elements (NMEs), after which we discuss the subamplitudes, A

B.1 Nuclear matrix elements
The NMEs used in Sect. 3.2 are defined as follows, Table 5 can be related to the previously discussed NMEs. For the magnetic GT matrix element we have, and the for short-distance NMEs

B.2 Hard neutrino exchange
The contributions due to hard neutrinos can be written as, in terms of the effective ππ, πN , and N N interactions, C ππ i , C πN i , and C N N i . A

(ν)
R (m i ) can be obtained by replacing C α i L → C α i R . In the case of the dimension-six couplings under consideration here, together with the approximation in which we drop the interference terms with the SM weak current as discussed in Sect. 3.2, several of the above terms vanish. It turns out that the derivative pion terms as well as the pion-nucleon couplings go to zero with the above approximations, so that we have C ππ i L = C πN i L = C πN i V = 0. The remaining C ππ α couplings are given by, Here, all LECs scale as g ππ i = O(F 2 π ) and the right-handed coupling c νππ i R can be obtained from c νππ i L by interchanging the L, R labels on the Wilson coefficients, L ↔ R, while leaving those on the LECs unchanged.   Table 5: Comparison of NMEs computed in the quasi-particle random phase approximation [52], shell model [51], and interacting boson model [66,67] for several isotopes of experimental interest. See e.g. [68] for a more recent calculation in the interacting boson model which, however, uses slightly differently defined matrix elements.
The relevant N N couplings can be written as, The right-handed couplings c NN,νNN i R can again be obtained from c NN,νNN i L by interchanging the L, R labels on the Wilson coefficients, L ↔ R, while leaving those on the LECs unchanged. Note that some of the terms in Eqs. (28) and (30) can be neglected given our assumptions, but can play a role in the right-handed couplings after the above replacement. Using naive-dimensional analysis would lead one to suspect that the g N N α scale as 1/Λ 2 χ , where Λ χ ∼ 1 GeV is the breakdown scale of chiral perturbation theory. However, it turns out that these LECs need to be enhanced by Λ 2 χ /F 2 π in order to obtain renormalized amplitudes, so that g N N α ∼ 1/F 2 π , see Ref. [23] for more details.
Similar to the m i dependence of the NMEs, we use the low-and high-mass limits of the LECs to construct interpolation formulae for these hadronic matrix elements. The interpolation we use is of the form g α (m i ) naive = g α (0) where m 0 2 GeV is a scale at which the procedure of integrating out the heavy neutrino becomes reliable. Theḡ α (m) are effective LECs, scaling as 1/m 2 , which are needed in the m i ≥ m 0 region to ensure the amplitude reproduces the result from integrating out the neutrino at the quark level. The above expression reduces to g α (0) for m i m π and m 2 0 m 2 iḡ α (m 0 ) for m i → ∞. The required input to these interpolation formulae are thus g α (0) andḡ α (m i ), where the latter depends on the matrix elements of the dimension-nine operators, ∼ (ūd) 2ē e c , that arise after integrating out neutrinos at the quark level. For simplicity we neglect the QCD evolution of these matrix elements, which can in principle be captured by Eq. (31).
Performing the matching to the amplitudes that result from integrating out neutrinos with m i 2 GeV gives rise to the following expressions for the ππ LECs Similar expressions can in principle be derived for the N N LECs [23], currently, however, none of the needed LECs are known in the N N sector. The techniques of Refs. [69,70] where the short-distance LECs for the LL case were calculated can in principle be applied to non-standard interactions as well, but this has not been done so far. In our numerical analysis we therefore estimate the hard-neutrino contributions by considering the ππ terms, while neglecting the N N pieces. We take the following values for the needed LECs, g ππ LR (0) = g ππ TT (0) = −g ππ S1 (0) = −g ππ S2 (0) = F 2 π , g N N α (0) = 0 , g ππ Here the values for g ππ α (0) are consistent with NDA estimates, while we use the lattice QCD results of Ref. [49] for g ππ 2,3,4,5 evaluated at µ = 2 GeV.