Long-distance Contributions to Neutrinoless Double Beta Decay $\pi^- \to\pi^+ e e$

Neutrinoless double beta decay, if detected, would prove that neutrinos are Majorana fermions and provide the direct evidence for lepton number violation. If such decay would exist in nature, then $\pi^-\pi^-\to ee$ and $\pi^-\to\pi^+ ee$ (or equivalently $\pi^-e^+\to\pi^+ e^-$) are the two simplest processes accessible via first-principle lattice QCD calculations. In this work, we calculate the long-distance contributions to the $\pi^-\to\pi^+ee$ transition amplitude using four ensembles at the physical pion mass with various volumes and lattice spacings. We adopt the infinite-volume reconstruction method to control the finite-volume effects arising from the (almost) massless neutrino. Providing the lattice QCD inputs for chiral perturbation theory, we obtain the low energy constant $g_\nu^{\pi\pi}(m_\rho)=-10.89(28)_\text{stat}(74)_\text{sys}$, which is close to $g_\nu^{\pi\pi}(m_\rho)=-11.96(31)_\text{stat}$ determined from the crossed-channel $\pi^-\pi^-\to ee$ decay.


I. INTRODUCTION
Observation of neutrinoless double beta (0ν2β) decays would prove neutrinos as Majorana fermions and lepton number violation in nature. The light-neutrino exchange is the most widely discussed mechanism to explain 0ν2β decays. Under this mechanism, the decay amplitude is proportional to the effective neutrino mass m ββ and thus the detection of 0ν2β decay would provide the information about the absolute neutrino mass, while the neutrino oscillation experiments are only sensitive to the mass differences among neutrinos.
On the theoretical side, 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% [15]. The interpretation of 0ν2β experiments relies on reliable calculations of the nuclear matrix elements, with robust uncertainty estimation. While the heavy nuclei system is beyond the capability of the current lattice QCD calculation, computations of the double beta decay for a light nuclei system shall be feasible [16,17]. The lattice results are required as inputs to determine the relevant low energy constants in the effective field theory [18][19][20][21], with which the nuclear matrix elements for heavy nuclei system can be calculated.
Without the signal-to-noise-ratio problem, the decay channels π − π − → ee and π − → π + ee serve as an ideal laboratory to perform a lattice QCD study of the 0ν2β decay and to test the prediction from effective field theory. Our exploratory study [2] has demonstrated the possibility of a first-principles calculation of the π − π − → ee decay, where we obtained the decay amplitude with A LO the leading-order prediction from chiral perturbation theory (χPT). By putting the amplitude into the χPT formula [18] A(ππ → ee) we obtain the low energy constant g ππ ν (µ) at µ = m ρ = 775 MeV where the uncertainties are statistical only. The two values of g ππ ν differ by ∼ 30%. This can be accounted by the systematic effects in the lattice calculation such as finite-volume effects and lattice artifacts, as well as higher-order truncation effects from χPT.
The lattice QCD calculations of π − → π + ee decay have been first carried out by CalLat Collaboration [22] for the short-distance contribution and NPLQCD Collaboration [23] for the long-distance contribution. While the vanishing phase space does not allow the π − → π + ee decay happen in nature (This problem does not exist for the K − → π + ee decay, which is proposed by Ref. [24]), the hadronic matrix element is well defined within the Standard Model and is equivalent to the one from π − e + → π + e − scattering, where π ± and e ± carry zero spatial momentum. As the crossed-channel analog to the π − π − → ee decay, the process of π − → π + ee can be combined together with π − π − → ee and serves as a cross-check for the prediction from χPT. Since in π − → π + ee decay the initial and final state only involves a single stable hadron, the study of the finite-volume effects is simplified. For example, we can adopt a newly developed technique called infinite-volume reconstruction [1] to determine the decay amplitude, where the finite-volume effects are exponentially suppressed even a massless neutrino propagator is included in the lattice calculation. Using four ensembles with different volumes and lattice spacings, we obtain the decay amplitude as where the first uncertainty is statistical and the second one is an estimation for both finitevolume effects and lattice artifacts. Using the χPT formula for π − → π + ee decay [18] A(π − → π + ee) we obtain the low energy constant Although the functional forms of χPT formulae (2) and (5) are quite different, the results for g ππ ν given in (3) and (6) are close to each other, demonstrating the success of χPT prediction.
The decay amplitude of a general 0ν2β process I(p I ) → F (p F )e(p 1 )e(p 2 ) can be written as A = F, e 1 , e 2 |H eff |I .
Here we use e 1,2 to specify the electron state carrying momentum p 1,2 . The second-order weak effective Hamiltonian is defined as Here G F is the Fermi constant and V ud is the CKM matrix element. The left-handed fermion fields are defined as ψ L = P L ψ,ψ L =ψP R (for ψ = u, d, e, ν e ) with projectors The effective Hamiltonian can be written as a product of hadronic and leptonic factors where the hadronic factor H µν (x) is defined as . Under the mechanism that 0ν2β decays are mediated by the exchange of light Majorana neutrinos, the leptonic factor L µν (x) is given by q 2 denotes a massless scalar propagator and the effective neutrino mass m ββ = | i U 2 ei m i | combines the neutrino masses m i and the elements U ei of the Pontecorvo-Maki-Nakagawa-Sato (PMNS) matrix. The charge conjugate of a fermionic field ψ is given as ψ c = Cψ T = γ 4 γ 2ψ T .
A. Decay amplitude of π − → π + ee For specific process π − → π + ee, with two electrons carrying vanishing momenta, the decay amplitude in Minkowski space-time becomes where leptonic part is factorized in T lept = 4G 2 F V 2 ud m ββēL (p 1 )e c L (p 2 ). The superscript M denotes the Minkowski space-time. The factor of 2 comes from interchange of electrons.

The hadronic function is defined by
This calculation is similar to the calculation of QED corrections to self energy using Feynman gauge. We can adopt the infinite-volume reconstruction (IVR) method proposed in Ref.
It should be noted that the hadronic function receives contribution from the vacuum state, which is lighter than the single pion state. In the Euclidean space-time, the hadronic function would grow exponentially as the time separation between the two current operators increases. To reproduce the amplitude in Eq. (11), one needs to treat the vacuum state properly as we will describe later.
In the following sections, we will first introduce our approach to calculate the Euclidean space-time hadronic function H(x) on lattice (For simplicity, we have left out the superscript of E for Euclidean space-time) and connect it with the Minkowski space-time integral, Eq.
(11). Then we use two different methods, QED L and IVR, to calculate the integral in Eq. (11). The results from the two methods are compared and discussed later.

B. Calculation of the hadronic function
In order to calculate the Euclidean space-time hadronic function on lattice, we define the following four point correlation function with wall-source pion interpolating operators φ † π − and φ π + . Here the time slices t i and t f are chosen as with sufficiently large ∆T for the ground-state saturation. Since the wall-source operators have a good overlap with the π ground state, we find the ground-state saturation for ∆T 1  Table I. Ensembles used in this work are generated by the RBC/UKQCD collaborations [25,26].
We list the pion mass m π , the space-time volume L 3 × T , the lattice spacing a, the number, N conf , of configurations used, the values of m π L and the time separation, ∆T , used for the π ground-state saturation.
fm at the physical pion mass. For the ensembles used in this work, the values of ∆T are chosen conservatively and listed in Table I.
The Euclidean space-time hadronic function is given by: where V = L 3 is the spatial-volume factor. In H(x − y) the vacuum-intermediate-state contribution from π + |J µL (x)|0 0|J µL (y)|π − leads to an exponentially growing factor e mπ(tx−ty) in Euclidean correlator when t x − t y increases. Here, we define the subtracted Euclidean space-time hadronic function: with H 0 (x) defined as The contractions of correlation function C(t f , x, y, t i ) are given by the type1 and type2 diagrams in Fig. 1, while the contractions of C 0 (t f , x, y, t i ) are given by the vacuum diagram, where the two hadronic parts are only connected through a neutrino propagator. It should be noted that either x t < y t or x t > y t are possible. After the subtraction, we remove the unphysical, exponentially-growing contributions. In Minkowski space-time, the amplitude contributed by the vacuum diagram, A M 0 , are known analytically and takes the simple form as where the decay constant F π is defined as Comparing A M 0 with the χPT formula [18] we can find that A M 0 is just the leading order term in χPT. In the Euclidean space-time, although H 0 function cannot reproduce the physical vacuum contribution A M 0 , by removing it and using the hadronic function H (x) as inputs, we can obtain the subtracted amplitude, A = A M − A M 0 , which includes the higher-order χPT contributions. Through out the paper, we will calculate the dimensionless, normalized which can be used to determine low energy constant g ππ ν via

C. Lattice setup
We use four ensembles at the physical pion mass generated by the RBC and UKQCD Collaborations [25,26]. The corresponding parameters are listed in Table I. For the ensembles 24D and 32D, lattice spacings are the same but the lattice volumes are different. It allows us to study the finite-volume effects. The ensemble 48I, 32D-fine and 24D have different lattice spacings but similar volumes. These ensembles provide us the information to examine the lattice artifacts. Note that the 48I uses Iwasaki gauge action in the simulation while the other three ensembles use Iwasaki+DSDR action.
We produce wall-source light-quark propagators on all time slices and make use of the time translation invariance to average the correlator over all T time translations. We have used AMA [27] and low modes deflation with compressed eigen-vectors [28]. These techniques have greatly reduced the computational cost for generating propagators. The correlators for the type1 and type2 diagrams in Fig. 1 are given by where S the light-quark propagator. We can write the above correlator in a general form of Such form allows us to obtain a spatial volume average of C(x) by using the double Fourier transformation. We have whereH i (t, p) (i = 1, 2) is a spatial Fourier transformation of H i (t, x). Using the spatial volume average, we can obtain a precise lattice data for both connected (type1) and disconnected (type2) diagrams.

D. The infinite volume reconstruction and QED L method
As pointed out earlier, the calculation of π − → π + ee transition is similar with the one of QED self energy contributions to π + -π 0 mass difference. So we adopt both methods of QED L and IVR in our analysis, although no QED effects are involved here. Note that a particularity of double beta decay is that the J µL interpolating operator involves both vector (V) and axial-vector (A) currents. The VA+AV contributions vanish due to the parity symmetry. After combining the VV and AA contributions, we obtain a large cancellation in π − → π + ee transition amplitude and a significant reduction in its uncertainty. As a result, the finite-volume (FV) effects are enhanced compared to the statistical errors.

QED L
We start the calculation of the normalized amplitude A defined in Eq. (21) using the QED L method, first introduced in Ref. [29] A where c 1 = 2.83729. For the AA part, the FV corrections arise from the excited states, e.g.

Infinite volume reconstruction
The detailed description of the IVR method have been given in Ref.
[1], where the time integral is split into the range of |t| > t s and |t| < t s In the fourth line of Eq. In this calculation, the lattice volumes listed in Table I are relatively small. Besides, FV effects are enhanced due to the cancellation between VV and AA contributions. As a consequence, even exponentially suppressed, the size of δ IVR (L) are statistically significant.
This can be confirmed by Fig. 2, where we compare the amplitude A IVR for ensembles 24D and 32D. Although VV and AA parts of the amplitudes are nearly consistent for the two ensembles, a significant discrepancy is found after the VV and AA parts are combined together.
where A π IVR and A π IVR (L) are determined using the hadronic functions H π (x) and H (L) π (x), respectively. Depending on the functional forms of F π (q 2 ), we have two estimates for δ π IVR .
• In the framework of scalar QED, we set F π (q 2 ) = 1, where the internal electromagnetic structure of pion is neglected. As shown in Fig. 3, at large t and x, e.g. t/a = 15 and | x|/a > 12, the lattice results of the VV hadronic function H V V lat (x) agree relatively well with the H (L) π (x) from scalar QED. We denote the δ π IVR using F π (q 2 ) = 1 as δ π,(1) IVR .
• If we adopt the expression of F π (q 2 ) = 1 + (r 2 π /6)q 2 and use the PDG value of the pion charge radius r π = 0.659(4) fm as input [32], a better consistency is found between the lattice data of H V V lat (x) and H (L) π (x) at t/a = 15. In this case, δ π IVR is denoted as δ π,(2) IVR .
From Fig. 3 we confirm that the long distance behavior of H V V lat (x) can be well described by the π intermediate state. Since FV effects mainly come from long distance physics, δ π IVR (L) can provide a good estimate for δ IVR (L). It is worthwhile to point out that both δ IVR (L) and δ π IVR (L) are exponentially suppressed as L increases, thus the residual FV effects in δ IVR (L) − δ π IVR (L) are also exponentially suppressed and thus well under control.  The IVR amplitudes A IVR as a function of t s are shown in Fig. 4 together with a fit to a constant. All the ensembles shown in Fig. 4 visibly agree with the corresponding fit in the window of 3 fm t s 4.5 fm and lead to reasonable values of χ 2 per degree of freedom.
The fitting results are shown in Table II.

A. Finite-volume effects
For ensembles 24D (m π L = 3.3) and 32D (m π L = 4.5), the results A IVR disagree by ∼ 10%, which is much larger than their statistical errors. We evaluate the FV corrections δ π,(1) IVR and δ π,(2) IVR at t s 3.75 fm by adopting Eq. (32) and using the input of F π (q 2 ) = 1 and F π (q 2 ) = 1 + (r 2 π /6)q 2 , respectively. As can be seen in Table II, after adding the corrections, the large discrepancy between 24D and 32D results vanishes. For the ensembles with smallest volume, e.g. 24D and 32D-fine, the results for A IVR + δ π,(1) IVR and A IVR + δ π,(2) IVR still differ by  In Fig. 5 we compare the results for ensemble 24D and 32D from QED L and IVR methods.
With LO and partially NLO FV corrections, the results of still have significant dependence on lattice volumes; while for the IVR results, such dependence is very mild. We therefore use the IVR results in the final analysis. For all the four ensembles, the amplitude A IVR + δ π,(1) IVR as a function of lattice spacing square is shown in Fig. 6. We use the results from 24D (a −1 = 1.015 GeV) and 32D-fine (a −1 = 1.378 GeV) ensembles to perform a continuum extrapolation and obtain A cont IVR = 0.1045(34) at the continuum limit. We have another ensemble 48I with a −1 = 1.73 GeV to further examine the lattice artifacts. As the 48I ensemble is simulated with Iwasaki gauge action, it contains the different lattice artifacts compared to the 24D and 32D-fine ensembles, where Iwasaki+DSDR action is used. Although the 48I result cannot be used in the continuum extrapolation directly, it helps us to estimate the size of the lattice artifacts.
the FV corrections δ π,(1) IVR or δ π,(2) IVR contributed by the pion intermediate state, the 24D and 32D results become consistent. The residual FV effects are estimated by the size of δ π,(1) IVR − δ π,(2) IVR . After the continuum extrapolation, the final result of A is given in Eq. (33), with a ∼ 3% statistical error. For the systematic effects including both FV effects and lattice artifacts, we estimate them at the level of ∼10%.
By putting the amplitude A into the χPT formula, we determine the low energy constant g ππ ν in Eq. (34). This result is close to the g ππ ν from π − π − → ee decay, which is given in Eq. (3), suggesting that χPT works well in the pion sector. It has been found in Ref. [20] that a leading-order, short-range contribution needs to be introduced in the nn → ppee decay, which breaks down Weinberg's power-counting scheme. Moving the calculation from the pion to the nucleon sector is the next step of our 0ν2β project. It is interesting to examine the impact of this short-range contribution in our future study.