Neutrinoless Double Beta Decay from Lattice QCD: The Long-Distance π − → π + e − e − Amplitude
MMIT-CTP/5196
Neutrinoless Double Beta Decay from Lattice QCD: TheLong-Distance π − → π + e − e − Amplitude
W. Detmold and D.J. Murphy (NPLQCD Collaboration) Center for Theoretical Physics, Massachusetts Institute of Technology, Boston, MA 02139, USA (Dated: April 17, 2020)Neutrinoless double beta decay (0 νββ ) is a hypothetical nuclear decay mode withimportant implications. In particular, observation of this decay would demonstratethat the neutrino is a Majorana particle and that lepton number conservation isviolated in nature. Relating experimental constraints on 0 νββ decay rates to theneutrino masses requires theoretical input in the form of non-perturbative nuclearmatrix elements which remain difficult to calculate reliably. This work marks a firststep toward providing a general lattice QCD framework for computing long-distance0 νββ matrix elements in the case where the decay is mediated by a light Majorananeutrino. The relevant formalism is developed and then tested by computing thesimplest such matrix element describing an unphysical π − → π + e − e − transition on aseries of domain wall fermion ensembles. The resulting lattice data is then fit to next-to-leading-order chiral perturbation theory, allowing a fully-controlled extractionof the low energy constant governing the transition rate, g ππν ( µ = 770 MeV) = − . stat (51) sys . Finally, future prospects for calculations of more complicatedprocesses, such as the phenomenologically important n n → p + p + e − e − decay, arediscussed. PACS numbers: 11.15.Ha, 12.38.Gc
I. INTRODUCTION
Neutrinoless double beta decay (0 νββ ), depicted in Figure 1, is a hypothetical nucleardecay process which, if observed, would provide a wealth of information about the propertiesof neutrinos. In particular, it is the only known experimentally viable method for resolvingthe long-standing question of whether neutrinos are Majorana or Dirac particles. In addition,it would also provide a first example of a lepton-number violating process, which may helpto explain baryogenesis in the early universe, as well as provide additional constraints on theparameters describing the neutrino sector in the Standard Model of particle physics. While0 νββ has not been observed, it is the subject of a large and active experimental search effort,with bounds on the half-lives of relevant nuclei at the level of T ν / (cid:38) − yrs [1].Next-generation experiments currently underway are aiming to probe half-lives that are anadditional one to two orders of magnitude larger in the near future [2].Relating a future experimental measurement of a 0 νββ decay rate T ν / for a particularnucleus to the effective Majorana neutrino mass m ββ = | (cid:80) k U ek m k | — where { m k } are the a r X i v : . [ h e p - l a t ] A p r u ud dd uν e e − e − ν e d ud du uW − W − n n p + p + (a) 2 νββ u ud dd ue − e − d ud du uW − ν e W − n n p + p + (b) 0 νββ FIG. 1. Quark-level Standard Model processes responsible for neutrinoful (left) and neutrinoless(right) double beta decay. neutrino eigenstate masses and U ek are elements of the Pontecorvo-Maki-Nakagawa-Sakata(PMNS) neutrino mixing matrix — requires theoretical input in the form of a nuclear matrixelement, M ν , describing the non-perturbative, hadronic part of the decay. These quantitiesare related by (cid:0) T ν / (cid:1) − ∝ | m ββ | G ν (cid:12)(cid:12) M ν (cid:12)(cid:12) , (1)where G ν is a known kinematic factor. Reliably calculating M ν for nuclear systems relevantto experimental searches has proven to be a difficult challenge. A variety of phenomenologicalnuclear models have been used to perform these calculations [3, 4], with predictions fora given nucleus from different models typically varying by 100% or more [5], and withno principled method for assigning systematic uncertainties. Improving this situation willbe crucial for interpreting experimental results from 0 νββ searches as constraints on theparameters of particular models of neutrinoless double beta decay moving forward.In principle, lattice QCD and the electroweak theory jointly provide an entirely ab-initio method for determining M ν . However, in practice, computing matrix elements of the largenuclei relevant to 0 νββ searches is well beyond the computational and algorithmic limits oflattice QCD for the forseeable future. More realistically, one could hope to compute QCDmatrix elements of sub-processes such as the n n → p + p + e − e − decay, and then relate theseto matrix elements of many-body systems within an effective field theory framework [6].Another possibility is to compute matrix elements of small nuclei which could then be usedto probe the systematics of nuclear model calculations by directly comparing lattice andmodel predictions.First calculations of the long-distance contributions to the neutrinoful double beta de-cay process n n → p + p + e − e − ν e ν e , and of the leading order short-distance contributions toneutrinoless double beta decay arising from new physics beyond the electroweak scale, werereported in Refs. [7] and [8], respectively. More recently, first calculations of the simplestlong-distance 0 νββ amplitude describing an unphysical π − → π + e − e − transition have ap-peared in the literature [9, 10], as well as a calculation of the related π − π − → e − e − decayamplitude [11].This work presents a complete calculation of the long-distance π − → π + e − e − amplitudeusing a series of domain wall fermion ensembles. The paper is organized as follows: SectionII and Appendix A develop the necessary formalism, including a novel treatment of the lightMajorana neutrino on the lattice using a regulated form of the continuum, infinite volumescalar propagator. Section III describes the lattice ensembles and numerical implementationsof the two- and four-point correlation functions needed to extract the π − → π + e − e − matrixelement, as well as a series of fits to next-to-leading-order chiral perturbation theory ( χ PT)used to extrapolate the lattice data to the physical mass, infinite volume, and continuumlimit, as well as determine the relevant χ PT low energy constant g ππν ( µ ) and assign a fullstatistical and systematic error budget. Finally, Sections IV and V discuss the results of thiscalculation in the context of other calculations in the literature, and lay out the prospectsfor future work. II. METHODOLOGY
It is assumed throughout this work that neutrinoless double beta decay is mediated bythe long-distance, light Majorana neutrino exchange mechanism. At low energies, and afterintegrating out the W boson, the underlying Standard Model interaction responsible forbeta decay is described by the effective electroweak Hamiltonian H W = 2 √ G F V ud ( u L γ µ d L ) ( e L γ µ ν eL ) , (2)where G F is the Fermi constant and V ud is the Cabibbo-Kobayashi-Maskawa (CKM) matrixelement describing the strength of the d → u transition arising from the flavor-changing weakinteraction. 0 νββ is induced at second order in electroweak perturbation theory, leading tothe bilocal matrix element [12] (cid:90) d x d y (cid:10) f ee (cid:12)(cid:12) T { H W ( x ) H W ( y ) } (cid:12)(cid:12) i (cid:11) = 4 m ββ G F V ud (cid:90) d x d y H αβ ( x, y ) L αβ ( x, y ) , (3)which can be factorized into tensors L αβ ≡ e L ( p ) γ α S ν ( x, y ) γ β e CL ( p ) e − ip · x e − ip · y ≡ Γ lept .αβ S ( x, y ) e − ip · x e − ip · y , (4)describing the leptonic part of the decay and H αβ ≡ (cid:10) f (cid:12)(cid:12) T { u L ( x ) γ α d L ( x ) u L ( y ) γ β d L ( y ) } (cid:12)(cid:12) i (cid:11) ≡ (cid:10) f (cid:12)(cid:12) T { j α ( x ) j β ( y ) } (cid:12)(cid:12) i (cid:11) , (5)describing the hadronic part of the decay, respectively. In addition, S ( x, y ) is the neutrinopropagator, e CL ≡ Ce (cid:62) L denotes charge conjugation, and T {· · · } denotes the time-orderingoperation. Since current constraints from oscillation experiments [13] suggest that m ββ isvery small compared to typical scales relevant to QCD or nuclear physics, it is also assumedthroughout this work that the massless scalar propagator S ( x, y ) = (cid:90) d q (2 π ) q e iq · ( x − y ) (6)is sufficient to describe the neutrino up to corrections which are much smaller than thepercent-scale statistical and systematic errors of the lattice calculations .To develop methodology, it is instructive to begin by considering the simplest 0 νββ process from the perspective of lattice field theory: an unphysical π − → π + e − e − transitionfor pions at rest. While this decay does not occur in nature, it is a well-defined amplitude Previous, exploratory work in Ref. [9] examined the neutrino mass dependence of the π − → π + e − e − am-plitude and found that the computed signals were indeed indisinguishable within statistical uncertaintieswhen m ββ (cid:28) m π . in quantum field theory, and serves as a natural starting point for lattice calculations sincesystems of pions are free of the well-known signal-to-noise issue plaguing calculations ofnucleon and nuclear systems [14, 15]. In addition, since this transition has only singlehadron initial and final states, the volume dependence of the hadronic matrix element isexpected to be mild and exponentially suppressed. The π − → π + e − e − amplitude has beencomputed at next-to-leading-order in chiral perturbation theory [6, 16–20], allowing for asimple case study of matching lattice results to χ PT in the context of 0 νββ amplitudes.The desired 0 νββ matrix element can be extracted in lattice QCD using methods whichhave been successfully applied to other second-order electroweak processes, including theneutrinoful double beta decay (2 νββ ) amplitude for the process n n → p + p + e − e − ν e ν e [7],as well as various kaon decays [21–23]. The key observation underlying these calculations isthat after integrating the Euclidean-space four-point function C π → πee ( t − , t x , t y , t + ) = (cid:88) (cid:126)x,(cid:126)y (cid:90) d q (2 π ) q e iq · ( x − y ) Γ lept .αβ (cid:10) O π + ( t + ) T { j α ( x ) j β ( y ) } O † π − ( t − ) (cid:11) (7)over t x and t y , the required matrix element M ν = ∞ (cid:88) n =0 (cid:88) (cid:126)x,(cid:126)y (cid:90) d q (2 π ) Γ lept .αβ (cid:10) πee (cid:12)(cid:12) j α ( (cid:126)x ) (cid:12)(cid:12) n (cid:11)(cid:10) n (cid:12)(cid:12) j β ( (cid:126)y ) (cid:12)(cid:12) π (cid:11) E n | (cid:126)q | ( | (cid:126)q | + E n − m π ) e i(cid:126)q · ( (cid:126)x − (cid:126)y ) (8)appears as the slope of the linear contribution in the T → ∞ regime: C π → πee ( T ) = T − ∆ (cid:88) t x =∆ T − ∆ (cid:88) t y =∆ C π → πee (0 , t x , t y , T )= ∞ (cid:88) n =0 (cid:88) (cid:126)x,(cid:126)y (cid:90) d q (2 π ) Γ lept .αβ (cid:10) πee (cid:12)(cid:12) j α ( (cid:126)x ) (cid:12)(cid:12) n (cid:11)(cid:10) n (cid:12)(cid:12) j β ( (cid:126)y ) (cid:12)(cid:12) π (cid:11) E n | (cid:126)q | ( | (cid:126)q | + E n − m π ) e i(cid:126)q · ( (cid:126)x − (cid:126)y ) × (cid:18) ( T − e − ( | (cid:126)q | + E n − m π )( T − − | (cid:126)q | + E n − m π (cid:19) , (9)where T is the size of the integration window, and n indexes all possible intermediate states.In practice, T is taken as large as possible to suppress the additional exponential and constantcontributions appearing in Eq. (9), and the cutoff ∆ (cid:28) T is chosen sufficiently large to avoidpotential coupling to excited initial or final states which may enter if the current insertionsare near the pion sources and sinks. At large T the matrix element (8) can be extractedfrom a simple linear fit to the T dependence of Eq. (9).The procedure described above can be spoiled by the appearance of long-distance inter-mediate states which would introduce exponentially growing, rather than exponentionallysuppressed, contamination into Eq. (9). This is certainly the case for the π − → π + e − e − transition, for which a pion-to-vacuum transition (cid:10) (cid:12)(cid:12) j µ (cid:12)(cid:12) π ( p ) (cid:11) = − ip µ f π (10)is allowed. A standard procedure for dealing with this contamination is to compute all suchtransition amplitudes on the lattice, allowing their contributions to be removed from thefour-point function (7) prior to performing the temporal integration, and thus removing theexponential divergence [23]. The contributions to the matrix element (8) from these low-lying states can then be reintroduced ex post facto . This particular aspect of the calculationis more difficult for Majorana exchange processes than for purely hadronic decays, since,in general, this subtraction would require the relevant first-order matrix elements to becomputed for the full range of momenta needed to saturate the integral over the neutrinomomentum (cid:126)q . Fortunately, the only relevant long-distance intermediate state for the π − → π + e − e − decay is the vacuum, for which the integration over (cid:126)q and the hadronic matrixelement decouple. For the phenomenologically important n n → p + p + e − e − decay, thisissue is likely avoided altogether, since the lightest long-distance intermediate state is thedeuteron [7], and the power-law fall-off of the neutrino propagator at large separations isexpected to overwhelm the exponentionally growing hadronic contribution arising from thesmall energy splitting between the dinucleon and deuteron states.A second complication in evaluating the four-point function defined by Eq. (7) on thelattice is that the continuum neutrino propagator (6) is divergent in the limit x → y . In thecontext of a lattice calculation, this divergence must be explicitly regulated, and a varietyof choices of this regulator have been explored in the literature. The approach taken inthe exploratory long-distance π − → π + e − e − calculation preceding this work [9], as well asin some lattice QCD+QED calculations which implement photon-exchange processes [24]and suffer from a similar divergence, is to use a lattice-regularized propagator with a non-zero bare mass, which can ultimately be extrapolated to zero. Another possibility exploredextensively in the lattice QCD+QED literature is to work directly with a massless latticepropagator after removing the divergent zero-mode contribution — one such example is thefirst-principles determination of the neutron-proton mass difference reported in Ref. [25] —which is known to introduce power-law finite volume effects [26, 27]. Yet another regular-ization scheme is the infinite volume reconstruction method introduced by Feng and Jin[28] and applied to neutrinoless double beta decay in Ref. [10]. In the present work analternative to these methods is explored: the neutrino propagator is implemented with aGaussian-regulated form of the continuum, infinite volume, massless scalar propagator, S Λ ( x, y ) = (cid:90) d q (2 π ) q e i · q ( x − y ) e − q / Λ , (11)which reduces to Eq. (6) in the limit Λ → ∞ . This approach has a number of advantages:Eq. (11) is computationally cheap and easily implemented, and does not introduce power-law finite volume effects. In addition, there is a natural choice of the regulator cutoff onthe lattice — Λ = π/a , where a is the lattice spacing — which ensures that the regulatoris removed in the continuum limit a → III. CALCULATION
The calculations detailed in this work make use of a series of N f = 2 + 1 domain wallfermion gauge field ensembles generated by the RBC/UKQCD collaboration and summarizedin Table I. These ensembles use the Iwasaki gauge action [29] and the domain wall fermionaction with the Shamir Kernel [30, 31] for the quarks. Each ensemble incorporates the seaeffects of two isospin-symmetric light quark flavors with bare mass am l and a single heavyquark flavor with bare mass am h . While the bare mass of the heavy flavor has been tunedto closely reproduce the physical strange quark mass, the bare masses of the light quarks aresomewhat heavier than the physical up and down quark masses, leading to simulated pionmasses in the range 300 MeV (cid:46) m π (cid:46)
430 MeV. The range of simulated masses, as wellas the two independent lattice spacings and physical volumes, allow for the 0 νββ matrixelement M ν to be matched to its predicted pion mass dependence from χ PT, as well as forthe results to be extrapolated to the infinite volume and zero lattice spacing limits. Detailsof the ensemble generation and fits to the low-energy spectrum are described in Refs. [32] and[33] for the 24I and 32I ensembles, respectively. The scale-setting analysis used to extractthe lattice cutoffs in physical units is described in Ref. [34].
Ensemble am l am s β L × T × L s m π L m π (MeV) a − (GeV)24I 0.01 0.04 2.13 24 × ×
16 5.81(1) 432.2(1.4) 1.784(5)24I 0.005 4.57(1) 339.6(1.2)32I 0.008 5.53(1) 410.8(1.5) 2.382(8)32I 0.006 0.03 2.25 32 × ×
16 4.84(1) 359.7(1.2)32I 0.004 4.06(1) 302.0(1.1)TABLE I. Summary of the ensembles and input parameters used in this analysis. Here, β is thegauge coupling, L × T × L s is the lattice volume decomposed into the length of the spatial ( L ),temporal ( T ), and fifth ( L s ) dimensions, and am l and am h are the bare, input light and heavyquark masses. Details of the ensemble generation and scale setting can be found in Refs. [32–34]. The remainder of this section describes the results of the calculations that were performed,as well as the fits that were used to extract physical quantities of interest. Section III Adescribes fits to two-point correlation functions used to extract the pion masses, decayconstants, and normalization factors of each simulation. Section III B describes fits to thefour-point function used to extract M ν . Finally, in Section III C, chiral perturbation theoryis used to extrapolate the lattice results for M ν to the physical pion mass, continuum, andinfinite volume limit, as well as to extract the relevant low energy constant g ππν ( µ ). A. Spectrum
Extracting M ν from the lattice four-point function (7) requires four inputs: the pionmass, the pion decay constant, the renormalizaton factor for the local V − A electroweakcurrent, and the pion-to-vacuum transition matrix element N s s O = (cid:10) (cid:12)(cid:12) O s s π (cid:12)(cid:12) π (cid:11) , (12)where O s s π are the pion interpolating operators with the same source ( s ) and sink ( s )smearing as used to compute the four-point function. These quantities can be determinedentirely from appropriate two-point functions, which, in this analysis, are constructed fromCoulomb-gauge fixed wall source lattice propagators computed using a deflated, mixed-precision conjugate gradient solver [35] with 1000 low-mode deflation vectors and a stoppingtolerance of r = 10 − . One such propagator is computed for each time slice, and the corre-lation functions are computed using both a local sink (L) and a zero-momentum projectedwall sink (W), and ultimately time-translation averaged over the entire lattice to improvethe signal. These techniques, as well as the details of the specific correlation functions andfitting procedures described below, have been developed and used previously in Refs. [34, 36],and will only be briefly discussed here.In this analysis six types of two-point functions are computed: the pseudoscalar-pseudoscalar correlator (cid:104) P P (cid:105) with the interpolating operator P ( x ) = q ( x ) γ q ( x ) and alocal or wall sink, the axial-pseudoscalar correlator (cid:104) AP (cid:105) with A µ ( x ) = q ( x ) γ µ γ q ( x ) and alocal or wall sink, and the correlators C A ( t ) ≡ (cid:68) (cid:12)(cid:12)(cid:12) (cid:88) (cid:126)x ∂ µ A µ ( (cid:126)x, t ) (cid:12)(cid:12)(cid:12) π (cid:69) (13)and C A ( t ) ≡ (cid:68) (cid:12)(cid:12)(cid:12) (cid:88) (cid:126)x ∂ µ A µ ( (cid:126)x, t ) (cid:12)(cid:12)(cid:12) π (cid:69) , (14)where A µ is the non-local, five-dimensional conserved axial current [36] and A µ is the local,four-dimensional axial current as defined above. The first four correlators can be used todetermine m π , f π , and the overlap factors N s s O by fitting the lattice results to the expectedtime-dependence of the ground states, (cid:10) (cid:12)(cid:12) O s s ( t )( O s s ) † (0) (cid:12)(cid:12) (cid:11) t (cid:29) (cid:39) N s s O N s s O † m π (cid:16) e − m π t ± e − m π ( T − t ) (cid:17) , (15)where the sign is + (-) for (cid:104) P P (cid:105) ( (cid:104) AP (cid:105) ), and extracting f π from the relation f π = 1 m π V |N W LA ||N
W LP ||N
W WP | . (16)The final two correlators involving the divergences of the axial currents are used to extractthe axial current renormalization coefficient Z A by fitting a constant to the ratio12 (cid:20) C A ( t −
1) + C A ( t )2 C A ( t − ) + 2 C A ( t ) C A ( t + ) + C A ( t − ) (cid:21) t (cid:29) (cid:39) Z A Z A . (17)In the remainder of this work it is assumed that Z A ≈ Z V ≈ Z A , so that therenormalization factor for the V − A electroweak current may also be approximated by Z A .These approximations are valid up to small O ( m res ) and O ( m ) corrections, respectively,where m res is the domain wall residual mass. Additional detail can be found in Refs. [34, 36].The fits are performed simultaneously to the four (cid:104) P P (cid:105) and (cid:104) AP (cid:105) two-point functions,as well as to the ratio defined in Eq. (17), by minimizing the fully correlated χ χ = N (cid:88) i =1 N (cid:88) j =1 (cid:16) y i − f ( x i ; (cid:126)β ) (cid:17) Σ − ij (cid:16) y j − f ( x j ; (cid:126)β ) (cid:17) , (18)where y i are the data, Σ ij = (cid:10) ( y i − (cid:104) y (cid:105) ) ( y j − (cid:104) y (cid:105) ) (cid:11) (19)is the covariance matrix computed from the data, and f is the assumed fit form dependingon independent variables x i and parameters (cid:126)β . The extracted parameters and corresponding χ / dof for each ensemble are summarized in Table II, and are found to be consistent with theprevious determinations in Refs. [32, 33]. In addition, plots of the data and the correspondingfits can be found in Appendix B. Ensemble am l |N W WP | am π af π Z A χ / dof24I 0.01 1 . × . × . × . × . × N W WP for the pseudoscalar interpolatingoperator P ( x ) = q ( x ) γ q ( x ) with a Coulomb gauge-fixed wall source and zero-momentum projectedwall sink (WW), the axial current renormalization factor Z A , and the correlated χ / dof for eachlattice ensemble. The errors are purely statistical and are computed using the jackknife resamplingtechnique. B. Long-Distance π − → π + e − e − Amplitude
Applying Wick’s theorem to the hadronic matrix element, Eq. (5), results in two classesof diagrams and four total contractions, depicted in Figure 2. In practice, computing these d uνxα, it − y β, jt + (a) Neutrino block(b) Type 1 contraction (c) Type 2 contraction (cid:104) S † u ( t − → x ) γ α (1 − γ ) S d ( t − → x ) (cid:105) · Tr (cid:104) S † u ( t + → y ) γ β (1 − γ ) S d ( t + → y ) (cid:105) (20)2 = Tr (cid:104) S † u ( t + → x ) γ α (1 − γ ) S d ( t − → x ) S † u ( t − → y ) γ β (1 − γ ) S d ( t + → y ) (cid:105) (21)FIG. 2. Top: diagrammatic representation of the neutrino block construction (Eq. (22)). Thelabels ( α, i ) and ( β, j ) reflect the open spin and color indices at the source and sink, respectively.Bottom: two classes of hadronic contractions for the π − → π + e − e − decay. Crossed circles denoteinsertions of the electroweak current. contractions by brute force is prohibitively expensive due to the double summation over thespacetime locations of the current insertions. In Ref. [7], this problem was solved for neutri-noful double beta decay amplitudes by computing quark propagators in the presence of anadditional background axial field, which can be shown to implicitly induce this summation.Unfortunately, background field techniques do not easily generalize to include the neutrinopropagator, requiring the development of other techniques for 0 νββ decays.A general method for computing 0 νββ contractions, including the diagrams in Figure2, with the full integration over the locations of both current insertions was introduced inRef. [9], and is briefly reviewed here. This method works by exploiting the convolutiontheorem and the translational invariance of the neutrino propagator to reduce the cost ofthe summation from O ( V ), where V is the lattice volume, to O ( V log V ) using the fastFourier transform (FFT). The key idea, depicted in the top panel of Figure 2, is to usethe FFT to integrate the leptonic tensor, Eq. (4), against the quark lines passing throughone of the two weak current insertions. More explicitly, for each fixed time ordering of theoperators a 12 ×
12 spin-color matrix-valued field B α ( x ; t − , t + ) = (cid:90) d y L αβ ( x − y ) (cid:104) S † u ( t − → y ) γ β (1 − γ ) S d ( t + → y ) (cid:105) = F − (cid:104) F ( L αβ ) · F (cid:0) S † u γ β (1 − γ ) S d (cid:1) (cid:105) (22)is computed using the FFT ( F ) and its inverse. The full contractions — for example,Eqs. (20) and (21) for the π − → π + e − e − transition — are then assembled by integratingthis “neutrino block” B α ( x ) against quark propagators to x and contracting the remainingopen indices in the appropriate combinations, for a total cost scaling as O ( V log V ). Scalingbenchmarks for this algorithm on CPUs and GPUs which demonstrate its efficiency werereported in Ref. [9].In this work, the four-point function, Eq. (7), is computed using the algorithm describedabove for all π − source and π + sink separations between 12 and 24 lattice units, and forall time orderings of the weak current insertions which are a distance of at least 6 latticeunits from the source and sink . In addition, on each gauge field configuration the type1 contraction is averaged over all time translations for each fixed time-ordering of the op-erators, while the more expensive type 2 contraction is time-averaged over four randomlychosen translations. The neutrino propagator is regulated using the UV cutoff imposed bythe lattice itself, Λ = π/a . The procedure for extracting M ν from this data is as follows:first, the fitted values of m π , f π , Z A , and N W Wπ from Section III A are used to remove theexponentionally divergent contribution from the vacuum intermediate state by subtractingEq. (38) from the data. Then, the subtracted four point function is normalized accordingto Eq. (35) to remove the dependence on the source-sink separation. Next, this normalizedfour-point function is integrated in the remaining time dependence of the current insertionsfor each source-sink separation, and a linear fit is performed to this signal at large separation T . The slope of this fit determines M ν up to the contribution from the vacuum interme-diate state ( M ν (0) ). Finally, this missing contribution is reintroduced using Eq. (42). Figure3 illustrates the vacuum subtraction procedure using data computed on the 24I am l = 0 . For the special case of the type 1 contraction (20) this cost can be reduced by a further factor of 144 sincethis contraction factorizes into two independent spin-color traces. Thus, it suffices to compute a scalarneutrino block rather than the full 12 ×
12 spin-color matrix This minimum separation corresponds to a physical distance of 0.7 (0.5) fm on the 24I (32I) ensembles.
10 5 0 5 10 ( t y − t x ) /a C π → π ee ( , t x , t y , ) n = 0 n ≥ n ≥ T/a π → π ee ( T ) n = 0 n ≥ n ≥ FIG. 3. Left (right): Example signals for the raw (integrated) four-point function, shown inblack (red) and defined in Eq. (7) (Eq. (9)), before and after subtracting the contribution fromthe vacuum intermediate state on the 24I am l = 0 .
01 ensemble. The gray band is the vacuumintermediate state contribution computed using Eq. (38) and the results of Section III A. to the choice of fit range within this window — the slopes extracted from linear fits tothe integrated four-point functions are observed to be somewhat sensitive to the choice offit range, while still maintaining acceptable χ / dof (cid:39)
1. To account for the systematicuncertainty in the choice of fit window, and to avoid potentially introducing a bias into theanalysis, the procedure for averaging over fits introduced in Refs. [37, 38] is adopted. Allpossible linear fits to the window [ T min /a, T max /a ] with 9 ≤ T min /a ≤
19, 21 ≤ T max /a ≤ χ (18).The parameters x = x i ± δx stat .i extracted from these fits, labeled by the index i , are thenaveraged according to ˆ x = (cid:80) i w i x i (cid:80) i w i , w i = p i ( δx stat .i ) , (23)where p i is the p -value corresponding to the χ /dof obtained in fit i , and w i is its associatedweight. This choice of weight is constructed to penalize both poor fits with small p i and fitswith large δx stat .i which poorly constrain the parameters. Uncertainties are assigned by alsoaveraging the statistical uncertainties from each fit (cid:0) δ ˆ x stat . (cid:1) = (cid:80) i w i ( δx stat .i ) (cid:80) i w i , (24)and computing the weighted average deviation from Eq. (23)( δ ˆ x sys . ) = (cid:80) i w i ( x i − ˆ x i ) (cid:80) i w i (25)as an estimate of the systematic uncertainty. Results for M ν obtained from this procedure,including both statistical and systematic uncertainties, are reported in Table III, along withthe dimensionless matrix element S ππ = M ν M ν (0) , (26)1in anticipation of the fits to χ PT discussed in the following section. Example fits to thedata for the window [ T min /a, T max /a ] = [13 ,
25] are shown in Figure 4.
Ensemble am l a M ν (0) a ( M ν − M ν (0) ) a M ν S ππ χ / dof24I 0.01 0 . . . . . . . . . . . . . . . . . . . . . . . . . νββ matrix element M ν from the vacuum inter-mediate state ( M ν (0) ), the correlated difference M ν − M ν (0) obtained from the slope of a linear fit tothe integrated four-point function C π → πee ( T ) in the limit T (cid:29)
1, the full matrix element M ν , andthe dimensionless matrix element S ππ = M ν /M ν (0) . The quoted uncertainties are statistical andsystematic, respectively, as computed by the fit averaging procedure described in the text. M ν (0) is assigned a systematic error of zero since it is computed from Eq. (42), and does not depend onthe linear fits performed in Section III B. T/a π → π ee ( T ) am l = 0 . am l = 0 . (a) 24I ensembles T/a π → π ee ( T ) am l = 0 . am l = 0 . am l = 0 . (b) 32I ensembles FIG. 4. Lattice signals and example fits to the window [ T min /a, T max /a ] = [13 ,
25] for the integratedfour-point function, Eq. (9). The data has been processed by first normalizing the raw four-pointfunction, Eq. (7), using Eq. (35), and then removing the exponentionally divergent contributionfrom the vacuum intermediate state using Eq. (38).
C. Chiral/Continuum Extrapolation
The final step in the calculation is to extrapolate the lattice data to the combined limitsof physical pion mass, zero lattice spacing, and infinite volume. These extrapolations are2performed simultaneously using the ansatz S ππ = 1 + m π π f π (cid:32) (cid:18) µ m π (cid:19) + 6 + 56 g ππν ( µ ) (cid:33)(cid:124) (cid:123)(cid:122) (cid:125) NLO χ PT + c NLO
F V e − m π L ( m π L ) / (cid:124) (cid:123)(cid:122) (cid:125) FV + c a a (cid:124)(cid:123)(cid:122)(cid:125) Continuum , (27)which includes the next-to-leading-order (NLO) pion mass dependence computed in χ PT[18], as well as models of the leading order discretization effects and finite volume effects.The term linear in a is motivated by the observation that the leading discretization artifactsenter at O ( a ) for domain wall fermions. The finite volume term is motivated by the leadingorder asymptotic expansions of the NLO χ PT finite volume corrections for f π [32] and thepion vector form factor [39], which enter as the n = 0 and n = 1 first-order hadronic matrixelements in Eq. (8), respectively. In both cases χ PT predicts∆
NLO
F V ∝ e − m π L ( m π L ) / , (28)up to higher order contributions suppressed by additional powers of m π L . In principle,the finite volume corrections could be computed self-consistently within the framework of χ PT, but this calculation has not been performed for the π − → π + e − e − amplitude in theliterature. We also consider a second, more general fit ansatz S ππ = 1+ m π π f π (cid:32) (cid:18) µ m π (cid:19) +6+ 56 g ππν ( µ ) (cid:33) + (cid:32) c NLO
F V ( m π L ) / + c NNLO
F V ( m π L ) / (cid:33) e − m π L + c a a , (29)which includes a model of the next-to-next-to-leading-order (NNLO) finite volume correc-tions. This generalized ansatz is only used for the purpose of studying fit systematicsassociated with the volume dependence.Table IV summarizes a variety of fits to the ans¨atze Eq. (27) and Eq. (29) using differentsubsets of the data, and the superjackknife resampling technique to propagate uncertaintiesfrom independent ensembles into a global fit [40]. The uncertainties in the inverse latticespacings used to convert to physical units (Table I), as well as the uncertainties in thephysical m PDG π − = 139 . f PDG π = 130 . g ππν is fixed at the conventional value µ = 770 MeV.Fits (A1)-(A4) are performed using the ansatz Eq. (27) and different cuts on the ensemblesincluded in the fit. It is observed that fit (A1) including all data has a poor χ / dof, arisingfrom tension between the data with the lightest pion mass and the data with the heaviestpion mass, but that the χ / dof improves significantly if either of these ensembles is pruned.While some improvement is observed in fit (A2), which prunes the ensemble with the heaviestpion mass and largest value of m π L , this is still a relatively poor fit with χ / dof = 3 .
8. Amuch more substantial improvement is observed in fit (A3), which instead prunes the lightestensemble with the smallest value of m π L , suggesting that residual finite volume effects drivethe observed tension rather than the truncation of the chiral expansion to NLO. Thus, fit(A3) is chosen as the preferred fit determining the central values and statistical errors of g ππν and the matrix elements S ππ and M ν at the physical point, and is depicted in Figure 5.3 Label m min π (MeV) m max π (MeV) g ππν ( µ ) c NLO FV c NNLO FV c a (fm − ) S phys .ππ M ν phys . (GeV ) χ / dofA1 302.0(1.1) 432.2(1.4) -10.71(11)(4) 0.3(3.2)(1.2) ≡ ≡ ≡ ≡ ≡ ≡ TABLE IV. Summary of fits of the lattice data reported in Table III to the ans¨atze Eq. (27) andEq. (29). Fits (A1)-(A4) use the ansatz Eq. (27) but vary the subset of the data included in the fit.Fit (B1) uses the same data and ansatz as fit (A3) but discards the infinite volume extrapolationby fixing c NLO
F V ≡
0. Finally, fit (C1) uses all of the available data and the more general ansatzEq. (29). The matrix elements S phys .ππ and M ν phys . are obtained by using the fit to extrapolate tothe physical m π and f π , as well as to zero lattice spacing and infinite volume. The χ PT LEC g ππν is determined at the renormalization scale µ = 770 MeV. The remaining fits (A4), (B1), and (C1) are variations on fit (A3) used to assign system-atic errors. Fit (A4) also prunes the heaviest mass ensemble from fit (A3) and is used toestimate the systematic error associated with truncating the chiral expansion to NLO. Thisresults in a fit with N dof = 0, so no χ / dof can be assigned. Fit (B1) uses the same dataas fit (A3) but removes the finite volume correction by fixing c NLO
F V ≡
0. Finally, fit (C1) isperformed to all of the data using the more general ansatz Eq. (29), and includes a second,higher-order finite volume correction term ∝ ( m π L ) − / . Including this additional term alsoresults in a good fit with χ / dof <
1, providing further evidence that the tension observedin fit (A1) is driven by residual finite volume effects.
D. Results and Error Budget
Based on the arguments presented in the previous section, fit (A3) in Table IV is chosenas the preferred fit to define the central values and statistical uncertainties of the mainresults of this work. In addition, the following systematic errors are estimated:1.
Sensitivity of the linear fits determining the lattice results for M ν to the choice offit range : In Section III B, a systematic uncertainty associated with the variation inthe extracted 0 νββ matrix elements as the fit window is varied is computed using aprocedure for averaging over possible fits introduced in Refs. [37, 38]. This systematichas been propagated through the chiral extrapolations performed in Section III C.2. Residual finite volume effects : Since the finite volume term included in the chiralansatz, Eq. (27), is a model rather than a quantity which has been computed self-consistently in χ PT, it is possible that the final results still contain residual finitevolume errors, and the fit variations studied in Section III C suggest that this is indeedthe dominant systematic uncertainty. Two procedures for estimating this systematichave been considered: the first, implemented as fit (B1), uses the same data as fit(A3) but drops the finite volume term altogether. The second, implemented as fit(C1), includes all of the available data and adds an additional term ∝ ( m π L ) − / modeling the neglected NNLO and higher order finite volume corrections. The larger4 m π (GeV ) S ππ χ PT24I Pred.32I a (fm ) S ππ Continuum24I Pred.32I ( m π L ) − S ππ FV24I Pred.32I a (fm ) S ππ am l = 0 . am l = 0 . am l = 0 . am l = 0 . FIG. 5. The chiral (top left), continuum (top/bottom right), and infinite volume (bottom left)extrapolations corresponding to the preferred fit (A3) in Table IV. In all but the bottom rightplot the fit has been used to shift the lattice data to the physical point ( m π = m PDG π − , f π = f PDG π , a = 0 , L = ∞ ), excluding the quantity specified on the horizontal axis. For the continuumextrapolation we also plot the raw data without correcting in m π , f π , and the lattice volume(bottom right) to illustrate that most of the uncertainty in the top right figure is associated withapplying this correction. The vertical dashed line in the upper left plot corresponds to the physical m π − = 139 . of the differences in central values between fits (A3) and (B1) or (C1) is used as aconservative estimate of this systematic.3. Truncation of the chiral expansion:
It is possible that higher-order terms in the chiralexpansion are needed to accurately describe the lattice simulations over the full rangeof pion masses reported in this work . One way to estimate the potential influenceof higher order terms is to successively prune the heaviest data from the chiral / It was found in Ref. [34], for example, that next-to-next-to-leading-order corrections to the quark massdependence of f π were needed to obtain a good fit describing a range of lattice data extending from thephysical point to the heaviest m π ≈
430 MeV 24I ensemble. g ππν (770 MeV) = − . stat (4) fit (50) FV (9) χ PT ,S ππ = 1 . stat (6) fit (61) FV (10) χ PT ,M ν = 0 . stat (2) fit (10) FV (2) χ PT GeV . (30) IV. DISCUSSION
The final results, including all sources of error — g ππν (770 MeV) = − . stat (51) sys and S ππ = 1 . stat (62) sys — are in good agreement with an independent lattice QCDstudy of the long-distance π − → π + e − e − amplitude by Tuo, Feng, and Jin [10], which de-termined g ππν ( m ρ ) = − . stat (74) sys and S ππ = 1 . stat (74) sys . This calculationalso used a variant of the domain wall fermion discretization for the quarks, but was per-formed on a different set of ensembles with near-physical pion masses and coarser a ≈ . L [27] and infinite volume reconstruction[28] techniques for this purpose. Since the calculation was performed at the physical pionmass, g ππν ( µ ) could be extracted directly by inverting S ππ = 1 + m π π f π (cid:18) (cid:18) µ m π (cid:19) + 6 + 56 g ππν ( µ ) (cid:19) , (31)rather than by performing a chiral fit as in Section III C of this work. The same authorsalso calculated g ππν ( m ρ ) = − . π − π − → e − e − decay amplitude inRef. [11], which is in ≈ σ tension with the determinations from π − → π + e − e − . This lattercalculation does not attempt to quantify any sources of systematic error, which, presumably,would help to explain the disagreement. Finally, in Ref. [18] Cirigliano et al. estimate g ππν ( m ρ ) (cid:39) − . χ PT [41, 42], which is also in reasonableagreement with the results presented here.One advantage of the approach taken in this work is that performing simulations at arange of different pion masses allows for a controlled study of how well NLO χ PT describeslattice data. Since connecting first-principles lattice QCD calculations to predictions forthe matrix elements of large nuclei used in 0 νββ searches will almost certainly involve ananalogous matching to an effective field theory — allowing for an extrapolation from thefew-body systems accessible on the lattice to the many-body systems relevant to experi-ment — this study is important to bridge from theory to phenomenology and experiment.Furthermore, lattice calculations of nuclear systems are currently performed at significantlyheavier than physical pion masses to ameliorate the signal-to-noise problem, making it cru-cial to understand how reliably such calculations can be matched to existing effective fieldtheory formalisms.6The chiral fits performed in Section III C exhibit a degree of tension between the latticedata — which spans the range of pion masses 300 MeV (cid:46) m π (cid:46)
430 MeV and volumes4 (cid:46) m π L (cid:46) χ PT amplitude, as well as models of the leading finite volume and discretiza-tion artifacts. Dropping the ensemble with the smallest m π L or adding an additional termparametrizing the neglected, higher-order finite volume corrections is observed to dramati-cally reduce this tension, resulting in good fits with χ / dof <
1, and suggesting that finitevolume artifacts are the dominant systematic uncertainty. In light of this observation, fu-ture lattice QCD calculations of long-distance neutrinoless double beta decay amplitudesmay benefit from self-consistently addressing the volume dependence within the effectivefield theory framework used to match to the lattice data, or from performing simulationswith sufficiently large volumes that finite volume artifacts are further suppressed. Importantformal work in this direction has been performed in Ref. [43].
V. CONCLUSIONS
In this work, a novel and general lattice QCD framework for computing four-point func-tions describing long-distance neutrinoless double beta decay amplitudes mediated by a lightMajorana neutrino has been presented and used to compute the amplitude for the unphys-ical π − → π + e − e − transition with pions at rest. Data for a series of lattice ensembles withpion masses in the range 300 MeV (cid:46) m π (cid:46)
430 MeV was fit to the next-to-leading-orderchiral perturbation theory amplitude for this decay, and the fit was used to predict the cor-responding matrix element in the physical mass, infinite volume, and continuum limit withsub-percent total uncertainty, as well as to determine the χ PT low energy constant g ππν with ≈
5% uncertainty. The results were found to be consistent with other estimates of thesequantities in the literature. Future work will apply these methods to the phenomenologicallyimportant n n → p + p + e − e − decay. ACKNOWLEDGMENTS
The authors wish to thank Z. Davoudi for pointing out to us in early 2015 the similarityof the 0 νββ matrix element to other matrix elements describing rare kaon decays, and foran earlier lattice formulation of the Euclidean correlation function for this problem. Theauthors also wish to thank X. Feng, L. Jin, E. Mereghetti, H. Monge-Camacho, A. Nicholson,A. Pochinsky, M. Savage, P. Shanahan, M. Wagman, A. Walker-Loud and B. Wang for usefuldiscussions. The calculations presented in this work were performed using the IBM BlueGene/Q computers of the RIKEN-BNL Research Center and Brookhaven National Lab andthe Stampede2 supercomputer of the Texas Advanced Computing Center (TACC) at theUniversity of Texas at Austin. In addition, computations for this work were carried out inpart on facilities of the USQCD Collaboration, which are funded by the Office of Science ofthe U.S. Department of Energy. WD and DJM are supported in part by the U.S. Departmentof Energy, Office of Science, Office of Nuclear Physics under grant Contract Number DE-SC0011090. WD is also supported within the framework of the TMD Topical Collaborationof the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, and by the7SciDAC4 award DE-SC0018121. [1] A. Gando et al. (KamLAND-Zen), Phys. Rev. Lett. , 082503 (2016), [Addendum: Phys.Rev. Lett.117,no.10,109903(2016)], arXiv:1605.02889 [hep-ex].[2] M. J. Dolinski, A. W. P. Poon, and W. Rodejohann, Ann. Rev. Nucl. Part. Sci. , 219 (2019),arXiv:1902.04097 [nucl-ex].[3] J. Engel and J. Menndez, Rept. Prog. Phys. , 046301 (2017), arXiv:1610.06548 [nucl-th].[4] J. D. Vergados, H. Ejiri, and F. Simkovic, Rept. Prog. Phys. , 106301 (2012),arXiv:1205.0649 [hep-ph].[5] A. Giuliani and A. Poves, Adv. High Energy Phys. , 857016 (2012).[6] V. Cirigliano, W. Dekens, J. de Vries, M. Graesser, and E. Mereghetti, JHEP , 097,arXiv:1806.02780 [hep-ph].[7] B. C. Tiburzi, M. L. Wagman, F. Winter, E. Chang, Z. Davoudi, W. Detmold, K. Orginos,M. J. Savage, and P. E. Shanahan, Phys. Rev. D96 , 054505 (2017), arXiv:1702.02929 [hep-lat].[8] A. Nicholson et al. , Phys. Rev. Lett. , 172501 (2018).[9] W. Detmold and D. Murphy,
Proceedings, 36th International Symposium on Lattice Field The-ory (Lattice 2018): East Lansing, MI, United States, July 22-28, 2018 , PoS
LATTICE2018 ,262 (2019), arXiv:1811.05554 [hep-lat].[10] X.-Y. Tuo, X. Feng, and L.-C. Jin, Phys. Rev.
D100 , 094511 (2019), arXiv:1909.13525 [hep-lat].[11] X. Feng, L.-C. Jin, X.-Y. Tuo, and S.-C. Xia, Phys. Rev. Lett. , 022001 (2019),arXiv:1809.10511 [hep-lat].[12] S. M. Bilenky and C. Giunti, International Journal of Modern Physics A , 1530001 (2015),https://doi.org/10.1142/S0217751X1530001X.[13] M. Tanabashi et al. (Particle Data Group), Phys. Rev. D , 030001 (2018).[14] G. Parisi, Phys. Rept. , 203 (1984).[15] G. Lepage (1989) pp. 97–120.[16] V. Cirigliano, W. Dekens, M. Graesser, and E. Mereghetti, Phys. Lett. B , 460 (2017),arXiv:1701.01443 [hep-ph].[17] V. Cirigliano, W. Dekens, J. de Vries, M. Graesser, and E. Mereghetti, JHEP , 082,arXiv:1708.09390 [hep-ph].[18] V. Cirigliano, W. Dekens, E. Mereghetti, and A. Walker-Loud, Phys. Rev. C , 065501(2018), [Erratum: Phys.Rev.C 100, 019903 (2019)], arXiv:1710.01729 [hep-ph].[19] V. Cirigliano, W. Dekens, J. De Vries, M. L. Graesser, E. Mereghetti, S. Pastore, andU. Van Kolck, Phys. Rev. Lett. , 202001 (2018), arXiv:1802.10097 [hep-ph].[20] V. Cirigliano, W. Dekens, J. De Vries, M. Graesser, E. Mereghetti, S. Pastore, M. Piarulli,U. Van Kolck, and R. Wiringa, Phys. Rev. C , 055504 (2019), arXiv:1907.11254 [nucl-th].[21] Z. Bai, N. H. Christ, T. Izubuchi, C. T. Sachrajda, A. Soni, and J. Yu, Phys. Rev. Lett. ,112003 (2014).[22] Z. Bai, Proceedings, 34th International Symposium on Lattice Field Theory (Lattice 2016):Southampton, UK, July 24-30, 2016 , PoS
LATTICE2016 , 309 (2017), arXiv:1611.06601[hep-lat].[23] Z. Bai, N. H. Christ, X. Feng, A. Lawson, A. Portelli, and C. T. Sachrajda, Phys. Rev.
D98 ,074509 (2018), arXiv:1806.11520 [hep-lat]. [24] M. G. Endres, A. Shindler, B. C. Tiburzi, and A. Walker-Loud, Phys. Rev. Lett. , 072002(2016), arXiv:1507.08916 [hep-lat].[25] S. Borsanyi et al. , Science , 1452 (2015), arXiv:1406.4088 [hep-lat].[26] A. Duncan, E. Eichten, and H. Thacker, Phys. Rev. Lett. , 3894 (1996), arXiv:hep-lat/9602005.[27] S. Uno and M. Hayakawa, Progress of Theoretical Physics , 413 (2008),https://academic.oup.com/ptp/article-pdf/120/3/413/5203205/120-3-413.pdf.[28] X. Feng and L. Jin, Phys. Rev. D100 , 094509 (2019), arXiv:1812.09817 [hep-lat].[29] Y. Iwasaki and T. Yoshie, Phys. Lett. B , 449 (1984).[30] D. B. Kaplan, Phys. Lett.
B288 , 342 (1992), arXiv:hep-lat/9206013 [hep-lat].[31] Y. Shamir, Nucl. Phys.
B406 , 90 (1993), arXiv:hep-lat/9303005 [hep-lat].[32] C. Allton et al. (RBC-UKQCD), Phys. Rev.
D78 , 114509 (2008), arXiv:0804.0473 [hep-lat].[33] Y. Aoki et al. (RBC, UKQCD), Phys. Rev.
D83 , 074508 (2011), arXiv:1011.0892 [hep-lat].[34] P. A. Boyle et al. , Phys. Rev.
D93 , 054502 (2016), arXiv:1511.01950 [hep-lat].[35] A. Stathopoulos and K. Orginos, SIAM J. Sci. Comput. , 439 (2010), arXiv:0707.0131[hep-lat].[36] T. Blum et al. (RBC, UKQCD), Phys. Rev. D93 , 074505 (2016), arXiv:1411.7017 [hep-lat].[37] E. Rinaldi, S. Syritsyn, M. L. Wagman, M. I. Buchoff, C. Schroeder, and J. Wasem, Phys.Rev.
D99 , 074510 (2019), arXiv:1901.07519 [hep-lat].[38] S. R. Beane et al. , (2020), arXiv:2003.12130 [hep-lat].[39] C. Alexandrou et al. (ETM), Phys. Rev.
D97 , 014508 (2018), arXiv:1710.10401 [hep-lat].[40] J. D. Bratt et al. (LHPC), Phys. Rev.
D82 , 094502 (2010), arXiv:1001.3620 [hep-lat].[41] B. Ananthanarayan and B. Moussallam, JHEP , 047, arXiv:hep-ph/0405206 [hep-ph].[42] R. Baur and R. Urech, Nucl. Phys. B499 , 319 (1997), arXiv:hep-ph/9612328 [hep-ph].[43] R. A. Briceo, Z. Davoudi, M. T. Hansen, M. R. Schindler, and A. Baroni, Phys. Rev. D ,014509 (2020), arXiv:1911.04036 [hep-lat].
A. FORMALISM
In this appendix, a derivation of the formalism describing how to extract the relevant0 νββ matrix element from a lattice calculation with a regulated, infinite volume, continuumneutrino propagator, Eq. (11), is outlined, beginning from the Euclidean four-point functiondefined in Eq. (7). While this derivation focuses on the π − → π + e − e − transition amplitude,the formalism generalizes straightforwardly to other 0 νββ processes.The first step in this calculation is to isolate the time dependence of the four-pointfunction arising from the leptonic and hadronic contributions, respectively. The Euclideantime dependence of the neutrino propagator can be extracted by performing the integrationover the temporal component of the neutrino’s four-momentum, which gives ∞ (cid:90) −∞ dq π | (cid:126)q | + q e iq ( t x − t y ) e − q / Λ = 14 | (cid:126)q | e | (cid:126)q | (cid:18) e −| (cid:126)q || t x − t y | Erfc (cid:20) | (cid:126)q | Λ − Λ2 | t x − t y | (cid:21) + e | (cid:126)q || t x − t y | Erfc (cid:20) | (cid:126)q | Λ + Λ2 | t x − t y | (cid:21)(cid:19) , (32)9where Erfc( x ) = 2 √ π ∞ (cid:90) x dt e − t (33)is the complementary error function. The time dependence of the hadronic matrix elementcan be extracted by inserting complete sums over eigenstates of the QCD Hamiltonian,Γ lept .αβ (cid:10) O π + ( t + ) T { j α ( (cid:126)x, t x ) j β ( (cid:126)y, t y ) } O † π − ( t − ) (cid:11) = ∞ (cid:88) l =0 ∞ (cid:88) m =0 ∞ (cid:88) n =0 Γ lept .αβ (cid:104) | O π + ( t + ) | l (cid:105)(cid:104) l | j α ( (cid:126)x, t x ) | m (cid:105)(cid:104) m | j β ( (cid:126)y, t y ) | n (cid:105)(cid:104) n | O † π − ( t − ) | (cid:105) E l E m E n (cid:39) |N π | m π e − m π | t + − t − | ∞ (cid:88) m =0 Γ lept .αβ (cid:104) πee | j α ( (cid:126)x ) | m (cid:105)(cid:104) m | j β ( (cid:126)y ) | π (cid:105) E m e − ( E m − m π ) | t x − t y | , (34)and assuming the current insertion time slices, t x and t y , are sufficiently separated from thesource and sink time slices, t − and t + , that the sums over l and n are saturated by theirrespective ground states. The dependence on the source-sink separation can be removed bydefining a normalized four-point function C π → πee ( t x , t y ) ≡ m π |N π | e m π | t + − t − | C π → πee ( t − , t x , t y , t + ) , (35)with m π and N π = (cid:10) (cid:12)(cid:12) O π (cid:12)(cid:12) π (cid:11) determined from the corresponding two-point function. Com-bining these results and relabling m → n results in C π → πee ( t x , t y ) = ∞ (cid:88) n =0 (cid:88) (cid:126)x,(cid:126)y (cid:90) d q (2 π ) Γ lept .αβ (cid:104) πee | j α ( (cid:126)x ) | n (cid:105)(cid:104) n | j β ( (cid:126)y ) | π (cid:105) E n | (cid:126)q | e i(cid:126)q · ( (cid:126)x − (cid:126)y ) e − ( E m − m π ) | t x − t y | × (cid:18) e −| (cid:126)q || t x − t y | Erfc (cid:20) | (cid:126)q | Λ − Λ2 | t x − t y | (cid:21) + e | (cid:126)q || t x − t y | Erfc (cid:20) | (cid:126)q | Λ + Λ2 | t x − t y | (cid:21)(cid:19) . (36)For the vacuum ( n = 0) intermediate state, the matrix elementΓ lept .αβ (cid:104) πee | j α (0) | (cid:105)(cid:104) | j β (0) | π (cid:105) E = m π f π Z A e L e CL (37)decouples from the remaining integration over the neutrino’s three-momentum, and theintegration can be performed explicitly. The result C (0) π → πee ( t x , t y ) = m π f π Z A (cid:88) (cid:126)x,(cid:126)y π | x − y | (cid:16) − e − Λ24 | x − y | (cid:17) e ( m π − m e ) | t x − t y | e L e CL (38)can be used to remove the contribution of the vacuum intermediate state to the four-pointfunction Eq. (7), and is manifestly finite in the limit x → y for finite cutoff Λ, with C (0) π → πee ∝ Λ , and exponentially divergent in the limit | t x − t y | → ∞ for m e < m π , as expected.To derive the analogue of Eq. (8), which is needed to compute the contribution of the vac-uum intermediate state to M ν , requires performing the time-ordered integration of Eq. (36)0in the operator insertion times. The finite sums over lattice times can be approximated asintegrals C π → πee ( T ) ≈ T (cid:90) dt x T (cid:90) t x dt y C π → πee ( t x , t y ) , (39)and the asymptotic behavior in the limit T → ∞ can be isolated using the expansionErfc( x ) = e − x x √ π ∞ (cid:88) n =0 ( − n (2 n − x ) n . (40)Keeping only the terms proportional to T results in M ν = ∞ (cid:88) n =0 (cid:88) (cid:126)x,(cid:126)y (cid:90) d q (2 π ) Γ lept .αβ (cid:104) πee | j α ( (cid:126)x ) | n (cid:105)(cid:104) n | j β ( (cid:126)y ) | π (cid:105) E n | (cid:126)q | (cid:2) | (cid:126)q | − ( E n − m π ) (cid:3) e i(cid:126)q · ( (cid:126)x − (cid:126)y ) × (cid:18) | (cid:126)q | e −| (cid:126)q | En − mπ )2Λ2 Erfc (cid:20) E n − m π Λ (cid:21) − ( E n − m π ) Erfc (cid:20) | (cid:126)q | Λ (cid:21)(cid:19) . (41)From this expression it is easily verified that Eq. (41) reduces to Eq. (8) in the limit Λ → ∞ ,while also rendering the matrix element finite for finite Λ.Using this expression, the contribution of any particular long-distance intermediate stateto M ν can be calculated provided one has calculated, or has otherwise modeled using experi-ment or phenomenology, the corresponding first-order hadronic matrix element as a functionof the three-momentum transfer (cid:126)q . For the vacuum intermediate state the hadronic matrixelement (37) and the integration over the momentum again decouple, and the contributionto M ν may be parametrized as M ν (0) = m π f π Z A f ( m π , L, Λ) e L e CL , (42)with f ( m π , L, Λ) =
P V ∞ (cid:90) dq (cid:40) (cid:88) (cid:126)x,(cid:126)y π | (cid:126)x − (cid:126)y | sin (cid:0) q | (cid:126)x − (cid:126)y | (cid:1) q − ( m π − m e ) × (cid:18) qe − q mπ − me )2Λ2 Erfc (cid:20) m e − m π Λ (cid:21) + ( m π − m e ) Erfc (cid:104) q Λ (cid:105)(cid:19) (cid:41) , (43)where P V denotes the Cauchy principal value . B. TWO-POINT FUNCTIONS
This section presents Figures 6-10, summarizing the fits to two-point functions performedin Section III A. Figures 6-9 show the effective pion masses am eff π = cosh − (cid:20) C π ( t −