Light-Neutrino Exchange and Long-Distance Contributions to $0\nu2\beta$ Decays: An Exploratory Study on $\pi\pi\to ee$

We present an exploratory lattice QCD calculation of the neutrinoless double beta decay $\pi\pi\to ee$. Under the mechanism of light-neutrino exchange, the decay amplitude involves significant long-distance contributions. The calculation reported here, with pion masses $m_\pi=420$ and 140 MeV, demonstrates that the decay amplitude can be computed from first principles using lattice methods. At unphysical and physical pion masses, we obtain that amplitudes are $24\%$ and $9\%$ smaller than the predication from leading order chiral perturbation theory. Our findings provide the lattice QCD inputs and constraints for effective field theory. A follow-on calculation with fully controlled systematic errors will be possible with adequate computational resources.

Introduction -It is a fundamental question whether the neutrinos are Dirac or Majoranna-type fermions. Neutrinoless double beta (0ν2β) decay, if detected, would prove that neutrinos are Majorana fermions. Besides, it provides direct evidence that the fundamental law of lepton number conservation is violated in nature. According to light-neutrino exchange mechanism, the observation of 0ν2β decay would also give us information about the absolute neutrino mass, which oscillation experiments cannot predict.
The standard picture of 0ν2β involves the long-range light neutrino exchange -a minimal extension of Standard Model. On the other hand, current knowledge of second-order weak-interaction nuclear matrix elements needs to be improved, as various nuclear models lead to discrepancies on the order of 100% [17]. A promising approach [18,19] to improving the reliability of the 0ν2β predication is to constrain the few-body inputs to ab initio many-body calculations using lattice QCD [20][21][22][23].
In this work we perform the first lattice QCD calculation of the non-local matrix elements for the process of ππ → ee, where the light neutrinos are included as active degrees of freedom. We find that the decay amplitude receives dominant long-distance contributions from eνπ intermediate state. Although small, the excited-state contribution is identified with a clear signal in our calculation. At both unphysical and physical pion masses, we find that the lattice results are consistently smaller than the predication from leading order chiral perturbation theory [18].
Light-neutrino exchange in 0ν2β decay -We begin with the effective Lagrangian L eff for the single β decay which represents the standard Fermi charged-current weak interaction involving the left-handed fermionic fieldsū L , d L ,ē L and ν eL . Here G F is the Fermi constant and V ud is the CKM matrix element. One can introduce the neutrino mixing matrix to connect the neutrino flavor eigenstates to the mass eigenstates. For the electron flavor, we havē with U ek the mixing matrix element. The effective Hamiltonian for 2β decay can be constructed as where the hadronic factor H µν ( . Under the mechanism that 0ν2β decays are mediated by the exchange of light Majorana neutrinos, the leptonic factor can be written as [24] L µν (x) = −m ββ S 0 (x, 0)e L (x)γ µ γ ν e c L (0) e iqx q 2 a massless scalar propagator and m ββ = ∑ k m k U 2 ek the effective neutrino mass. The charge conjugate of a fermionic field ψ is given as ψ c = Cψ T = γ 4 γ 2ψ T .

arXiv:1809.10511v1 [hep-lat] 27 Sep 2018
For a general 0ν2β decay I(p I ) → F (p F )e(p 1 )e(p 2 ), its decay amplitude can be written as Here we use e 1,2 to specify the electron state carrying momentum p 1,2 . The leptonic matrix element is given by which is antisymmetric under the exchange of two electrons e 1 ↔ e 2 due to the Pauli exclusion principle. Here the spinors are defined as for i = 1, 2. Inserting the complete set of hadronic intermediate states, the decay amplitude can be written as Given the spatial momenta ⃗ p for the hadronic intermediate states specified by n, the neutrino's momenta are constrained by the conservation law ⃗ p ν,i = ⃗ p I − ⃗ p − ⃗ p i and the corresponding energies are denoted as E ν,i = ⃗ p ν,i . One can write the spinor product as a combination of . The coefficient of the second term is proportional to the difference in electron momenta and generically suppressed by a fac- MeV is the typical Fermi momentum of nucleons in a nucleus [18,24]. Keeping only the term of u L (p 1 )u c L (p 2 ), the decay amplitude is simplified as . Calculation of ππ → ee decay -In this work we calculate the ππ → ee decay amplitude with two pions at rest and two electrons carrying spatial momenta ⃗ p 1 = −⃗ p 2 , ⃗ p 1,2 = E ππ 2. While the condition of ⃗ p 1 − ⃗ p 2 k F ≪ 1 is no more valid, we target on the determination of the amplitude given in Eq. (9), which is more relevant for chiral effective field theory inputs to ab initio many-body calculation [18]. This setup has advantages as follows.
• Due to the non-zero momentum carried by the electron, the energies of any possible intermediate states eνn always lie above the initial-state energy E ππ ≈ 2m π . Therefore no exponentially growing contamination is associated with the intermediate states when one performs an integral over a Euclidean time. The effects of finite volume on a generic second-order weak amplitude [25] is not relevant here as well.
• We use the discrete lattice momenta (2π L) ⃗ m for the intermediate hadronic particles and the mo- p i is the momentum carried by the electron. As non-zero momenta are assigned for the neutrino propagator, one can keep the lowest mode of the propagator, which reduces the powerlaw finite-volume effects.
Note that no short-distance divergence appears as x approaches to 0 in L eff (x)L eff (0). This can be seen by the power counting in the integral In lattice QCD, a hard cutoff is introduced by the inverse of lattice spacing 1 a. Thus the unphysical short-distance contribution appears as an O(a 2 ) discretization effect. Using the Coulomb gauge fixed wall sources for the φ π interpolating operator, we construct the correlation function through C(t x , t y , t ππ ) = On the lattice, the scalar propagator withq i = 2 sin(q i 2) the lattice discretized momenta. V and T are the spatial volume and time extent of the lattice. We can calculate the zero mode (⃗ q = 0) of the propagator as 1 V T ∑ qt e iq t (tx −ty ) qt 2 +∑ iki modes (⃗ q ≠ ⃗ 0) of the propagator can be constructed as 1 Nr ∑ Nr r=1 φ r (x)φ * r (y) using stochastic method, with Here the stochastic sources ξ r (q) satisfy It is proposed by NPLQCD collaboration that the neutrino propagator can also be computed in an exact way by using double Fourier transformation [26]. Following Refs. [27][28][29][30][31][32][33] and integrating t x and t y over a fixed window [t a , t b ] with t a ≫ t ππ , we obtain with T box = t b − t a + 1 the time extent of the integration window. N ππ and E ππ are known from the correlation function ππ 2Eππ e −Eππt + e −Eππ(T −t) + constant by using the methods proposed in Ref. [34]. When T box is sufficiently large, the contamination from the exponential term vanishes as E n + E ν,i + E i > E ππ . The coefficient of the term proportional to T box provides a result for the decay amplitude A(ππ → ee).
Numerical results -We use two ensembles with m π = 420 and 140 MeV generated by the RBC and UKQCD collaborations [35]. The corresponding parameters are listed in Table I. We produce wall-source lightquark propagators on all time slices and make use of the time translation invariance to average the correlator over all T time translations. (To reduce the computational costs at m π = 140 MeV, we adopt the technique of all mode average [36,37] with T sloppy propagators used for correlator average and 1 precise propagator for correction.) We compute propagators for both periodic and anti-periodic boundary conditions in the temporal direction and use their average in the calculation, which effectively doubles the temporal extent of the lattice.
The Feynman diagrams corresponding to the process of ππ → ee are shown in Fig. 1. To show the time dependence of the C(t x , t y , t ππ ) explicitly, we define the unintegrated amplitude M(t) as a function of the variable t = t x − t y :  Table I. Ensembles used in this work. We list the pion mass mπ, the lattice spacing inverse a −1 , the space-time volume L 3 × T , the number, N conf , of configurations used and the number, Nr, of stochastic sources for the neutrino propagator. The time t x and t y are separated by at least 6 time units from the ππ sources (t x,y −t ππ ≥ 6) so that the φ † π φ † π interpolating operators can project onto the ground ππ state. At large t , the time dependence of M(t) is saturated by the ground intermediate state -eνπ where the matrix elements of ⟨0 J µL π⟩ V and ⟨π J µL ππ⟩ V are determined from the correlation functions ⟨J µL (t)φ † π (0)⟩ and ⟨φ π (t π )J µL (t J )φ † π φ † π (t ππ )⟩, respectively. The subscript ⟨⋯⟩ V indicates that the initial and final states are defined in the finite volume. The single-pion states π⟩ V satisfy the normalization condition ⟨π(⃗ p) π(⃗ p ′ )⟩ FV = (2E π )V δ ⃗ p, ⃗ p ′ , while the two-pion states ππ⟩ V can be connected to the states in the finite volume ππ⟩ ∞ through Lellouch-Lüscher relation [38,39] with the momenta k = E 2 ππ 4 − m 2 π and q = kL (2π). The time dependence of M(t) is shown in Fig. 2. At large t the data of M(t) is consistent with the contribution from the ground intermediate state. By subtracting the ground-state contribution, the remaining excitedstate contribution is shown by blue square points in the left panel of Fig. 2 and enlarged in the right panel. Although relatively small, the contribution from the excited intermediate states can be identified with a clear signal.
The integrated matrix element defined in Eq. (15) is shown in Fig. 3. We realize that the size of integration window T box ≈ 16 is not sufficiently large to discard the exponential term associated with the ground intermediate state. (This can be confirmed in Fig. 2 that at t ≈ 15 the values of M(t) are statistically larger than 0.) After removing this exponential term, we can fit the lattice   data to a linear function of T box and determine the values of A lat (ππ → ee). To convert A lat (ππ → ee) to the physical amplitude A(ππ → ee), a renormalization factor square Z 2 V A shall be multiplied, which relates the local lattice vector or axial-vector current (which we use) to the conserved or partially conserved ones. Besides, the Lellouch-Lüscher factor shall be multiplied to relate a finite-volume amplitude to the infinite-volume one. In our calculation, the two pions are in the ground state, i.e. at threshold. The large-L expansion of the Lellouch-Lüscher factor is given by with a ππ the scattering length and r ππ the effective range from the k-expansion of ππ scattering phase shift k cot δ(k) = a −1 ππ + r ππ The coefficients d i are given by The values of zeta function Z 00 (s, 0) have been provided by Ref. [40]. We evaluate a ππ using Lüscher's finite-size formula [40] and use it as an input to determine the finitevolume correction up to O(L −2 ). Another type of power-law finite-volume effect arises from the long-range property of the neutrino propagator. The finite-volume effects relevant for the eνπintermediate state can be evaluated as n the discrete momentum for the pion. The neutrino's energy is given by E ν = ⃗ p + ⃗ p e , with ⃗ p e the momentum carried by the electron. We define a function f (⃗ p) ≡ and split it as The term inside brackets does not contribute a power-law finite-volume effect. We thus simplify ∆ FV as The function κ(⃗ n e ) with ⃗ n e = ⃗ p e L (2π) can be computed numerically and we find κ(⃗ n e ) = 0.686(3) for m π = 420 MeV and 0.517(3) for m π = 140 MeV. Thus Eq. (22) indicates that the finite-volume correction appears as an O(L −2 ) effect. We expect that the size of f (−⃗ p e ) is significantly smaller than f ( ⃗ 0), as the total contribution to the decay amplitude from the intermediate hadronic states that carry non-zero lattice momenta only amounts for 3-4% when compared to the zero-momentum contribution. We therefore neglect this finite-volume effect in this work, and leave it for future studies.
In Table II, we show the ground-state, excited-state and total contributions to the decay amplitude as A (g) , A (e) and A (g) + A (e) , respectively. The results are presented in units of F 2 π T lept , where the decay constant F π is determined from the matrix element ⟨0 d γ µ γ 5 u π(p)⟩ = √ 2p µ Z A F π , with Z A the renormalization constant. Systematic effects associated with three choices of t a − t ππ = 6, 7, 8 are relatively smaller than the statistical errors, suggesting that a separation of 6 is a safe choice to neglect the excited ππ states.
0.055 (13) Table II. Results for ground-state (A (g) ), excited-state (A (e) ) and total (A (g) + A (e) ) contributions to the ππ → ee decay amplitude. All the results are listed in units of F 2 π T lept .
The decay amplitude of A(ππ → ee) is mainly contributed by the ground intermediate state via the process of ππ → πeν → ee. Although the size of excited-state contribution is only 3-4%, it is statistically significant (see Fig. 2) as the uncertainty of the amplitude has been reduced to below 1%. Without the signal-to-noise problem, the case of ππ → ee serves as an ideal laboratory to develop the necessary methods and tools for a calculation of 0ν2β decay with controlled uncertainties. Our exploratory study demonstrates the possibility of a first-principles calculation of the long-distance contribution to 0ν2β decay via light-neutrino exchange. At m π = 420 and 140 MeV, we find that the decay amplitude A(ππ → ee) are 24% and 9% smaller than the leading-order predication A LO (ππ → ee) = 2 F 2 π T lept in chiral perturbation theory [18]. Various systematic effects such as lattice artifacts and finite-volume effects require an accurate examination in the future work but are not expected to qualitatively alter the conclusions of this work. The techniques presented here can be directly applied to the study of other 0ν2β decays, such as nπ → pee and nn → ppee. From these decays, lattice QCD can provide more lowenergy QCD inputs for the effective field theory [18].