Relativistic model-free prediction for neutrinoless double beta decay at leading order

Starting from a manifestly Lorentz-invariant chiral Lagrangian, we present a model-free prediction for the transition amplitude of the process $nn\rightarrow pp e^-e^-$ induced by light Majorana neutrinos, which is a key process of the neutrinoless double beta decay ($0\nu\beta\beta$) in heavy nuclei employed in large-scale searches. Contrary to the nonrelativistic case, we show that the transition amplitude can be renormalized at leading order without any uncertain contact operators. The predicted amplitude defines a stringent benchmark for the previous estimation with model-dependent inputs, and greatly reduces the uncertainty of $0\nu\beta\beta$ transition operator in the calculations of nuclear matrix elements. Generalizations of the present framework could also help to address the uncertainties in $0\nu\beta\beta$ decay induced by other mechanisms. In addition, the present work motivates a relativistic {\it ab initio} calculation of $0\nu\beta\beta$ decay in light and medium-mass nuclei.


Introduction
Neutrinoless double beta decay (0νββ) is a secondorder weak process, in which a nucleus decays to its neighboring nucleus by turning two neutrons to two protons, emitting two electrons but no corresponding antineutrinos [1].Its observation would provide direct evidence of lepton number violation beyond the standard model, prove the Majorana nature of neutrinos [2,3], constrain the neutrino mass scale and hierarchy [4], and shed light on the matter-antimatter asymmetry in the Universe [5].Therefore, it becomes one of the top priorities in the field of nuclear and particle physics, and invigorates worldwide experiments either activate or planned [6,7,8,9,10,11,12,13] (for a recent review, see Ref. [14]).
Nuclear matrix element, which encodes the impact of the nuclear structure on the decay half-life, is crucial to interpreting the experimental limits and even more potential future discoveries.Within the standard picture of 0νββ, which involves the long-range light neutrino exchange [15], a minimal extension of the standard model, current knowledge of the nuclear matrix element is not satisfactory [16], as various nuclear models lead to discrepancies as large as a factor around 3.
Chiral effective field theory (EFT) [17,18] plays an important role in addressing such uncertainties.It can provide the nuclear Hamiltonian and weak currents in a consistent and systematically improvable manner [19,20,21].Based on chiral EFT, the issue of g A quenching in single β decays [22] was resolved as a combination of two-nucleon weak currents and strong nuclear correlations [23,24].The chiral-EFT-based 0νββ transition operators under various sources of lepton number violation [25,26,27,28,29,30,31], as well as their impacts on the nuclear matrix elements have been extensively studied in the literature [32,33,34,35,36,37].In addition, ab initio calculations of nuclear matrix elements using a chiral Hamiltonian are available recently [38,39,40,41].
In the context of the standard light Majorana neutrinoexchange mechanism [15], Cirigliano et al. showed that a chiral EFT description of 0νββ decay, already at leading order (LO), requires a contact operator to ensure renormalizability [28,29].The size of this contact term should be determined in principle by matching to first-principle gauge field theory calculations, which however are not yet available currently (see Refs. [42,43,44,45,46,47,48,49,50] for related progress).This unknown contact term, thus, leads to an additional source of uncertainty for the nuclear matrix elements besides the nuclear-structure ones.
Cirigliano et al. [30,31] first modeled the nn → ppe − e − transition amplitude with an integral representation, i.e., a momentum integral of the neutrino propagator times the generalized forward Compton scattering amplitude.This approach introduces, in the intermediate-momentum region, model-dependent inputs from elastic intermediate states in analogy to the Cottingham formula [51].In the absence of experimental datum and costly calculations based on lattice gauge theory, the obtained model-dependent amplitude could be used as a pseudo datum to determine the contact term in nuclear structure calculations [41].However, a model-free prediction of the transition amplitude is still missing.
In this Letter, we present a model-free prediction of the nn → ppe − e − amplitude by developing a relativis-tic framework based on chiral EFT.The term "modelfree" here means the obtained results are free of modeldependent inputs beyond the framework of the chiral EFT.We show that the transition amplitude can be renormalized at LO without the vague contact term in the relativistic framework starting from a manifestly Lorentz-invariant chiral Lagrangian.Matching to our results will allow one to stringently determine the contact term needed in nonrelativistic nuclear-structure calculations, and to greatly reduce the corresponding uncertainty in the calculations of nuclear matrix elements for 0νββ decay.It also paves the way to relativistic ab initio calculations of 0νββ decay rates in light and medium-mass nuclei.

Relativistic framework
We start from a manifestly Lorentz-invariant effective Lagrangian relevant at the leading order of chiral EFT [20], where f π = 92.2MeV is the pion decay constant, g A = 1.27 is the nucleon axial coupling, M denotes the nucleon mass, and C α (α = S, V, A, AV, T ) are the low-energy constants (LECs).This Lagrangian consists of the pion field u = exp[i τ • π/(2f π )] and the nucleon field Ψ = (p, n) T , which are coupled to the weak current l µ via the axial vector u µ = iu † (∂ µ − il µ )u − iu∂ µ u † and the chirally covariant derivative The weak current reads l µ = −2 √ 2G F V ud τ + e L γ µ ν eL + h.c., with the Fermi constant G F and the V ud element of the Cabibbo-Kobayashi-Maskawa (CKM) matrix [52,53].
For the standard mechanism of 0νββ decay, the lepton number violation at low energy is dominated by the electron-neutrino Majorana mass where C = iγ 2 γ 0 denotes the charge conjugation matrix, and m ββ the effective neutrino mass.Contrary to the previous studies based on the heavy baryon approach [28,29], which relies on a nonrelativistic expansion of the Lagrangian, here we apply the manifestly Lorentz-invariant Lagrangian to the problem of 0νββ decay.The LO contribution to the scattering amplitude can be obtained by solving the relativistic scattering equation where E is the total energy, and p ′ and p are the nucleon outgoing and incoming momenta in the center of mass frame, respectively.
Equation ( 3) is consistent with the three-dimensional reduction of the Bethe-Salpeter equation [54] and it satisfies relativistic elastic unitarity.The LO potential V = ϕ 0 ⊗ ϕ 0 Vϕ 0 ⊗ ϕ 0 is defined by the LO two-nucleon irreducible diagrams sandwiched between the leading term of the Dirac spinor ϕ(p, s) expanding in powers of small momenta p.The subleading terms of the Dirac spinor and the retardation effects on pion and neutrino exchanges are to be included perturbatively at higher orders.As a result, the LO strong and neutrino potentials take the form as those in the Weinberg's approach, where q = p ′ − p, and Note that the present derivation is similar to the socalled modified Weinberg approach [55], which was applied to nucleon-nucleon scattering problem.It has proven to be useful to improve the renormalizability of nucleon-nucleon scattering [55] and few-body systems [56,57].In the heavy baryon approach, the nonrelativistic expansion of the Lagrangian leads, instead of Eq. ( 3), to the Lippmann-Schwinger equation, which contains a nonrelativistic limit (1/M → 0) of the two-nucleon propagator (6) Obviously, the relativistic propagator scales as O(Λ −3 ) at the ultraviolet (UV) region |k| ∼ Λ, while the nonrelativistic one scales as O(Λ −2 ).As a result of this milder UV behavior, we will show that the 0νββ amplitude can be renormalized without promoting a contact term to the LO neutrino potential.

Renormalization of the 0νββ amplitude
We focus on the scattering process nn → ppe − e − in the 1 S 0 wave, as it is the only channel that requires a contact term to ensure the renormalizability in the heavy-baryon approach [28,29].Apart from the 1 S 0 channel, there are other channels with nonzero angular momenta that have fairly large contributions to the nuclear matrix elements in nuclear structure calculations [58,59], but the LO amplitude can be renormalized in these channels without the contact terms.
Without loss of generality for our arguments, we consider the kinematics n(p i )n(−p i ) → p(p f )p(−p f )e(p e1 = 0)e(p e2 = 0).The LO amplitude can be schematically written as where ) is a phase space factor1 , G 0 the two-nucleon free propagator, and T s the T -matrix resumming the strong potential V s .The four terms in Eq. ( 7) correspond to the four diagrams depicted in the first row of Fig. 1, and here we denote them as A A , A B , A B , and A C from left to right.
To study the renormalization of the transition amplitude, we now discuss the degree of UV divergence for the diagrams A ∼ O(Λ D ) (log Λ for D = 0).The A A tree diagram has no divergence.For A B , A B , and A C , it was shown that the degree of divergence is dominated by the loop integrals involving the insertion of neutrino potential, once T s is made finite [28].The divergence structures are shown in the second row of Fig. 1.Counting the powers of loop momenta, one finds the degree of divergence with L being the number of loops and g the UV scaling of the two-nucleon propagator, and −2 comes from the |q| −2 dependence of V ν .From Eq. ( 6), we know the nonrelativistic propagator has g = −2 and the relativistic one has g = −3.Therefore, in the nonrelativistic framework, A B and A B are convergent as O(Λ −1 ), but A C is logarithmic divergent.In the present relativistic framework, however, A B , A B , and A C are all convergent as O(Λ −2 ), so the corresponding LO amplitude is renormalizable.Such analysis should hold true not only for chiral EFT, but also for pionless EFT, since the dominating UV divergence structure does not involve pion exchanges.
The renormalizability of the LO amplitude A LO ν can be demonstrated explicitly in any specific regularization scheme.Here, we regulate the strong potential with a separable gaussian function, After projecting to the 1 S 0 channel, the LEC in its shortrange part C 1 S0 = C 1 − 3C 2 is determined by reproducing the scattering length a np = −23.74fm.We have checked that the 1 S 0 phase shifts are indeed cutoff independent as Λ → ∞ .Then, the amplitude is evaluated with Eq. ( 7) in the momentum space. is plotted against the cutoff Λ in the upper panel of Fig. 2. A logarithmic divergence of the amplitudes can be clearly seen for the results given by the nonrelativistic pionless and chiral EFTs, consisting with the analysis above.This leads to the introduction of a contact term to ensure the renormalizability, as proposed in Ref. [28].In contrast, the amplitudes obtained in the relativistic framework, both chiral and pionless EFTs, converge against the increasing cutoff.This means that the amplitudes can be predicted at LO without any undetermined LECs at all kinematics.
In the lower panel of Fig. 2, the leading-order contributions to the amplitude (see the first row of Fig. 1) are depicted separately, as functions of the cutoff.The two-loop diagram A C dominates the cutoff dependence of the full LO amplitude in the present relativistic framework, while it is logarithmic divergent in the nonrelativistic case.
Note that recent studies on the renormalization of chi-ral EFT in a scheme with finite cutoffs have been applied to nucleon-nucleon scattering [60,61], and extending such a scheme to 0νββ decay process could also be interesting in the future, and the present results will provide a benchmark for such studies.One can define the LO 0νββ operator for calculating nuclear matrix elements.In the momentum space, it is defined as V LO ν = −A A with A A the tree-level amplitude [62], Here, V ν (p ′ , p) is the LO neutrino potential [Eq. ( 5)].After expanding the operator in terms of 1/M , one finds the leading relativistic correction has the opposite sign from the nonrelativistic part.The relativistic correction grows with increasing momentum cutoff and is significant at large cutoffs Λ M .In Fig. 3, we benchmark our calculated amplitude against the one estimated previously at the kinematics |p i | = 25 MeV and |p f | = 30 MeV [30,31].Such a comparison is not trivial at all, because the two amplitudes are obtained in distinct approaches.The amplitude obtained in this work is properly renormalized by the relativity and there is no need to introduce any unknown counter terms, while in the previous estimation [30,31], one has to renormalize the amplitude by introducing an unknown contact term, whose size was constrained by additional model-dependent inputs within a certain range.
It is remarkable that the renormalized amplitude obtained with the relativistic chiral EFT, is quite consistent with the previous estimation [30,31], and their difference is only about 10%.This demonstrates the validity of both approaches for 0νββ transitions, but the present relativistic framework avoids the uncertainties introduced by the model-dependent inputs.In addition, the theoretical uncertainties in the present framework can be systematically improved by moving to higher orders of the chiral expansion.The theoretical uncertainty originating from the truncation of the chiral expansion was rarely discussed in previous works on the nn → ppe − e − amplitude.Following Ref. [63], a very rough estimation could be given by m π /Λ b with Λ b = 600 MeV the breakdown scale of chiral expansion, i.e., around 23%.Compared to the relativistic chiral EFT, its pionless counterpart predicts an amplitude larger by about 40%.We would expect that the predictions of the amplitudes with chiral and pionless EFTs should come closer at higher orders, as is the case for N N scattering [64,65], and it would be interesting to examine this in the future.[30].With the fitted contact term, the nonrelativistic amplitudes at other kinematic points can also be renormalized.The uncertainty of the synthetic datum reflects its systematic error from the model dependences, which include the neglect of inelastic intermediate states, the parametrization of momentum dependence of nucleon-nucleon scattering amplitude, and the phenomenological form factors [30,31].These uncertainties are in turn propagated to the contact term required in the nonrelativistic chiral EFT by fitting to the synthetic datum, and also the resultant amplitudes at all other kinematic points.In contrast, the present relativistic framework provides a model-free prediction of the amplitudes, and they are generally consistent with the nonrelativistic estimations.The relativistic results are slightly larger than the nonrelativistic ones by about 10%-20% for the on-shell amplitude, while for the off-shell amplitude, the enhancement becomes more sizable especially at high final momentum, e.g., about 40% at |p f | = 150 MeV.This should be understandable because the contact term in the nonrelativistic framework is determined at extremely lowenergy kinematics.
In the present relativistic framework, the contact term is expected to appear at next-to-next-to-leading order, so its contribution is significantly suppressed as compared to the nonrelativistic case.Therefore, the present results are model-free at LO, in the sense that the present approach is free of model-dependent inputs beyond the framework of chiral EFT.There exist some variants of power counting schemes in both nonrelativistic and relativistic EFTs [21], but it not preclude the chiral EFT from providing model-independent results for physical observables.Such scheme dependence could in principle be controlled by the EFT paradigm within the truncation errors.
The present work provides the opportunity for a direct benchmark between the chiral EFT predictions and the lattice QCD calculations for 0νββ decay, which is currently only possible at unphysically heavy quark masses due to the computational cost [66,67,47].Extending the present framework to meet the lattice QCD conditions is straightforward, since it only takes the axial coupling constant g A and the 1 S 0 scattering length a np as inputs.A direct comparison against the upcoming lattice QCD simulations for the nn → ppe − e − amplitudes, instead of using them as inputs, will provide a stringent test on the validity of the present framework.
The present results can help to address the uncertainty of nuclear matrix elements in nuclear-structure calculations for isotopes of experimental interest for 0νββ searches.On the one hand, the present amplitudes can serve as alternative synthetic data to determine the contact term in nonrelativistic nuclear-structure calculations.Note that the contact-term effects could be amplified in realistic 0νββ transitions [28].On the other hand, a consistent relativistic ab initio calculation of nuclear 0νββ decay rates starting from relativistic chiral forces is highly demanded in view of the recent significant progress in both relativistic ab initio methods [68,69,57] and nuclear forces [70,71,72,73].
It is interesting to extend the present framework to investigate other processes related to 0νββ decay, e.g., the two-neutrino double-beta decay and the pion double charge exchange reactions.For the cases of these two processes, the strong interactions in the intermediate state should be taken into account at LO, in contrast to the case of 0νββ decay.As a result, the renormalization-group analysis for the two-nucleon amplitudes for these processes could be more involved, due to the tensor forces generated by the one-pion exchange in the intermediate states.

Summary and Outlooks
We present a model-free prediction of the nn → ppe − e − amplitude by developing a relativistic framework based on chiral EFT.Contrary to the nonrelativistic case, we show that the amplitude can be renormalized at LO without any uncertain contact operators.The calculated amplitude is slightly larger than the previous model-dependent estimation at the kinematics |p i | = 25 MeV and |p f | = 30 MeV [30,31] by about 10%, providing a highly nontrivial validation for the previous result.The enhancement of the amplitude could be more sizable for either off-shell amplitude or high-momentum kinematics.Therefore, the present amplitude defines a stringent benchmark for the previously estimated amplitude, and can help to address the uncertainty of the nuclear matrix elements in nuclearstructure calculations for isotopes of experimental interest for 0νββ searches.In addition, the present work motivates a relativistic ab initio calculation of 0νββ decay in light and medium-mass nuclei.

Figure 1 :
Figure 1: Leading-order contributions to the amplitude of nn → ppe − e − (first row) and the corresponding ultraviolet divergence structures (second row).The double and plain lines denote nucleon and lepton fields, respectively.The squares denote an insertion of neutrino potential Vν .The circles denote the nucleon axial and vector currents coupled to Vν .The gray ellipses represent the T matrix generated by iteration of the strong potential Vs.

Figure 2 :
Figure 2: (Color online).The upper panel depicts the amplitudes A LO ν for |p i | = 1 MeV and |p f | = 38 MeV, as functions of the cutoff Λ.The solid and open circles (squares) denote the results of the relativistic (nonrelativistic) chiral and pionless EFTs, respectively.The lower panel depicts the leading-order contributions to the amplitude (see the first row of Fig. 1) obtained with the relativistic chiral EFT.Taking the kinematics |p i | = 1 MeV and |p f | = 38 MeV as an example, the LO amplitude A LO ν

Figure 3 :
Figure 3: (Color online).The amplitude A LO ν for |p i | = 25 MeV and |p f | = 30 MeV given by the relativistic chiral and pionless EFTs, in comparison with the result estimated previously in Ref. [30].

Figure 4 :
Figure 4: (Color online).The renormalized on-shell (upper) and off-shell (lower) amplitudes A LO ν at various kinematic points.The blue band represents the nonrelativistic results, in which the contact term is fitted to the synthetic datum (shown as diamond) for the amplitude at kinematics |p i | = 25 MeV and |p f | = 30 MeV, i.e., A LO ν = −0.0195(5)MeV −2 [30].The width of the band reflects the uncertainty propagated from the synthetic datum.

Figure 4
Figure 4 depicts the renormalized on-shell and off-shell amplitudes at various kinematic points.The results given