Lattice QCD calculation of $\pi^0\rightarrow e^+ e^-$ decay

We extend the application of lattice QCD to the two-photon-mediated, order $\alpha^2$ rare decay $\pi^0\rightarrow e^+ e^-$. By combining Minkowski- and Euclidean-space methods we are able to calculate the complex amplitude describing this decay directly from the underlying theories (QCD and QED) which predict this decay. The leading connected and disconnected diagrams are considered; a continuum limit is evaluated and the systematic errors are estimated. We find $\mathrm{Re} \mathcal{A} = 18.60(1.19)(1.04)\,$eV, $\mathrm{Im} \mathcal{A} = 32.59(1.50)(1.65)\,$eV, a more accurate value for the ratio $\frac{\mathrm{Re} \mathcal{A}}{\mathrm{Im} \mathcal{A}}=0.571(10)(4)$ and a result for the partial width $\Gamma(\pi^0\to\gamma\gamma) = 6.60(0.61)(0.67)\,$eV. Here the first errors are statistical and the second systematic. This calculation is the first step in determining the more challenging, two-photon-mediated decay amplitude that contributes to the rare decay $K\to\mu^+\mu^-$.

Introduction. -Recent advances in the methods of lattice QCD allow the calculation of increasingly complex processes which involve both QCD and QED. While initially motivated by the calculation of the light-by-light scattering process which contributes to the anomalous magnetic moment, g µ − 2, of the muon [1,2], these methods can be applied to other processes where an ab initio lattice QCD result would be of value. Of particular interest is the rare flavor-changing neutral current (FCNC) decay K L → µ + µ − . This decay results from both an accurately-calculated short-distance second-order weak process and a largely unknown electroweak process involving a two-photon intermediate state. A comparison of the second-order weak prediction with experiment is important for testing the standard model and exploring physics beyond it. However, to compare this standard model prediction with experiment, we require a firstprinciples calculation of the long-distance, two-photon contribution to this decay.
In contrast with light-by-light scattering entering g µ −2 which can be computed using Euclidean-space methods, K L → µ + µ − is inherently a Minkowski process with a complex amplitude resulting from physical time evolution. Here we develop a method that allows us to deal with such a process using lattice QCD and apply this method to calculate the complex π 0 → e + e − decay amplitude -a process which can be viewed as a simpler version of the K L → µ + µ − decay but suffering from the same problem of a potentially on-shell, two-photon intermediate state. It is the lattice QCD calculation of the π 0 → e + e − decay which is the subject of this paper.
To leading-order the π 0 → e + e − decay is mediated by a two-photon intermediate state as shown in Fig. 1. This decay is highly suppressed by the chiral symmetry of the electromagnetic interaction of the electron which results in a decay rate proportional to (m e /M π ) 2 where m e and M π are the masses of the electron and π 0 . This decay amplitude A is complex with an imaginary part determined by the optical theorem, giving the wellknown unitary bound for the π 0 → e + e − branching ratio. We will use lattice QCD to compute both the real and imaginary parts of A using a treatment of QED without power-law finite-volume errors. Thus, our result for the imaginary part provides an improved calculation of the partial width π 0 → γγ with finite volume errors that vanish exponentially with increasing volume, without the constraint to use photon momenta that obey quantization conditions related to the volume used in the lattice QCD calculation [3]. (Preliminary results from these methods appear in Refs. [4][5][6][7].) The π 0 → e + e − decay has been measured to near 1% accuracy in the KTeV experiment [8]. After the contribution of e + e − γ Dalitz decays have been removed, the conventional O(α) radiative corrections must be per-arXiv:2208.03834v2 [hep-lat] 23 Aug 2022 formed. These O(α) corrections [9, 10] require a twoloop calculation in which the zeroth-order π 0 → e + e − decay is not approximated as point-like. As is conventional, we use the calculation of Refs. [9,10] to remove these radiative effects from the experimental results, giving B π 0 → e + e − = 6.86(27) stat. (23) syst. × 10 −8 . This modified experimental result can be compared with our lattice QCD calculation which also does not include these corrections. These corrections depend on a low-energy constant χ (r) = 4.5, determined self-consistently from the experimental result with these corrections removed.
This corrected experimental value is larger than most theoretical results including that obtained here, B π 0 → e + e − = 6.22(5) stat (2) syst × 10 −8 . See Ref. [11] for a recent precise theoretical prediction and extensive references to earlier results.
Computational strategy.
-We start with the Minkowski-space expression for the decay amplitude: where w = u − v is the relative space-time position of the two electromagnetic currents, P = ( 0, M π ) is the four-momentum of initial pion and h = ± 1 2 is the helicity of the electron on which A does not depend. Note that we have integrated out their average position u+v 2 and removed the resulting delta function that imposes total energy and momentum conservation. The Minkowski metric tensor g µµ = diag(−1, 1, 1, 1).
A direct Euclidean-space calculation of the hadronic matrix element in Eq. 1 using lattice methods is feasible if we Wick-rotate the w 0 contour by replacing the real variable w 0 with the product e −iφ w 0 and increasing φ from 0 to π/2 while keeping w 0 real, thereby making w 0 a Euclidean time. At the same time, the p 0 contour must also be rotated so that the exponent ip 0 w 0 in Eq. 1 remains purely imaginary as |p 0 | → ∞. Because of the existence of the two-photon intermediate state whose energy may be lower than the energy of initial pion state, the rotated p 0 contour cannot simply follow the imaginary axis. Instead the p 0 contour must be distorted as in Fig. 2 to avoid the poles which cross the imaginary axis when the energy of the two-photon intermediate state is lower than the pion mass, 2| p| < M π .
This choice of contour guarantees that the p 0 integral will converge. However, once the variable w 0 has become imaginary, the factor e −ip·w will grow exponentially with |w 0 | for the real values of p 0 that appear for the contour of Fig. 2. Fortunately, the Euclidean-space hadronic matrix element H µν (w) will compensate for this behavior, so that both the p 0 and w 0 integrations are convergent. The dependence of the hadronic matrix element on w 0 can be determined by inserting a complete set of intermediate states between the two electromagnetic currents in the hadronic matrix element 0|T The lightest intermediate state is the two-pion state with E n = 2M π . Thus, this Euclidean-space hadronic matrix element will decay as exp(−3M π |w 0 |/2) for large |w 0 |. This large-|w 0 | fall-off is sufficient to overcome the exp (|w 0 |M π /2) growth due to the p 0 contour in Fig. 2. We can express the analytic expression for decay amplitude after contour deformation as follows: where g µµ = diag(i, 1, 1, 1) is introduced to connect the Minkowski conventions for the electromagnetic currents in L µν with the Euclidean conventions used in the hadronic matrix element. Note that the amplitude A is not altered by these changes of contour and remains complex with real and imaginary parts which can be computed from the Euclidean-space amplitude in Eq. (5).
The leptonic factor L µν is evaluated by performing the p 0 integral using Cauchy's theorem, resulting in a threedimensional integral over p. When integrated over | p | the singular factor 1 | p|−Mπ/2−i gives both real and imaginary parts. Recognizing that the integral is independent of the direction of the outgoing positron momentum k + , allows us average over this direction. Here we present the result for the spatial components, writing L ij (w) = L(w 0 , | w|) ijk w k /| w| 2 with Here β = 1 − 4m 2 e /M 2 π , F (x) = cos(x) − sin(x)/x and E e (k, p, θ) = k 2 + p 2 − 2pk cos(θ) + m 2 e . We then evaluate the leptonic factor L µν as a two-dimensional numerical integral, requiring that the integration error lies below 0.001%. The result is tabulated as a function of w 0 and | w|. We evaluate the four-dimensional integral over w in Eq. (3) as a sum over lattice points with the values of L µν (w) obtained from this table by linear interpolation.
Hadronic matrix element. -The hadronic matrix element H µν in Eq. (5) can be calculated from the threepoint function: which can be computed using lattice QCD. It is only in this evaluation of H µν that a finite volume is introduced. Since this Euclidean-space amplitude involves no massless particles, all finite-volume errors will decrease exponentially in the linear size of the volume. The factor Z V renormalizes the non-conserved, local current used in the lattice calculation of the right hand side of Eq. (8). The factor N π is obtain from the π 0 twopoint function and compensates for the normalization of the pion interpolating operator π 0 (t). We use a Coulombgauge-fixed wall source for this operator.
The two diagrams needed to evaluate this three-point function are shown in Fig. 3. The first is connected and can be constructed from two wall-source propagators and one point-source propagator. The two electromagnetic currents must be located sufficiently far from the pion interpolating operator that contamination from states more energetic than the pion can be neglected. We require the pion source at t to be separated from the closer current by a fixed time ∆t. Let T be time extent of the lattice and require that all time coordinates lie in the interval [0, T ). Label the time locations of current operators t > and t < where (t > − t < )%T < T /2. We then average the results from two choices for the time t of the pion wall source. For the first we chose t = (t < − ∆t)%T and for the second t = (t > + ∆t)%T .
The second diagram is quark-line-disconnected. Such disconnected diagrams typically involve large statistical noise and are difficult to calculate. In our calculation, we use for the quark loop shown in Fig. 3b the results for Tr D −1 (x, x)γ µ computed using all-to-all propagators computed with randomly-displaced, 3 4 grid sources from the RBC/UKQCD calculation of the disconnected contribution to the hadronic vacuum polarization component of g µ − 2 [12]. As shown in Tab. III, we determine this disconnected amplitude with a statistical error of 60%.
The hadronic matrix element is calculated on four ensembles, whose parameters are listed in Tab. I. All en-sembles use Möbius domain wall fermions (DWF), which achieve good chiral symmetry with a much smaller size in the 5th dimension than required by the conventional DWF action. All ensembles use the Iwasaki gauge action. The 24ID, 32ID and 32IDF ensembles use the dislocation suppressing determinant ratio action to reduce chiral symmetry breaking effects. For every configuration, we have 1024 or 2048 point-source propagators with randomly distributed sources and Coulomb-gauge-fixed wall-source propagators with sources on every time slice.  Table I. Table of gauge-field ensembles used here. These were generated by the RBC/UKQCD collaborations [13].
Results. -The calculated real and imaginary parts of the amplitude are listed in Tab. II and plotted in Fig. 4. In Tab. II the experimental value for the imaginary part is evaluated using the optical theorem and the experimental pion life time; the experimental real part is obtained by subtracting the imaginary part contribution from the experimental decay rate, with radiative corrections. In the table and plot only the contribution from the connected diagram is included. The contribution from the disconnected diagram is treated as a source of systematic error. The disconnected diagram is calculated for the 24ID ensemble and the result shown in Tab. III.

Source
Im We use the continuum limit extrapolated from the 48I and the 64I ensembles as shown in Fig 5 as   contributions of connected and disconnected diagrams for the 24ID ensemble listed in Tab. III. The error resulting from our pion mass being 4 MeV larger than physical is estimated from chiral perturbation theory while that arising from our choice of Z V is taken from Ref. [13].
The smaller error on this ratio results from statistical correlations and our method of estimating the systematic errors. We can combine our more accurate result for this ratio with the experimental result for the decay width Γ(π 0 → γγ) to obtain more accurate values for the real part of the decay amplitude and the branching ratio [14]: (10) B(π 0 → e + e − ) = 6.22(5) stat (2) syst × 10 −8 . (11) The error labeled "expt" arises from the error on the measured π 0 → γγ decay rate. This result for the branching ratio is 1.8 sigma below the experimental value [8] B(π 0 → e + e − ) expt = 6.86(27) stat. (23) syst. × 10 −8 , from which radiative corrections have been removed. Finally we can use our result for the imaginary part of A to determine a lattice QCD prediction for the decay width Γ(π 0 → γγ) = 6.60(0.61) stat (0.67) syst eV to be compared with the experimental result [15] Γ(π 0 → γγ) = 7.802(0.052) stat (0.105) syst eV. This new lattice QCD result is computed with physical quark masses, contains finite volume errors which are suppressed exponentially in the linear extent of the lattice and can be viewed as a refinement of earlier lattice results [16,17].
Conclusion and Outlook. -We have applied a combination of covariant Feynman perturbation theory and lattice QCD to calculate the complex amplitude describing the decay π 0 → e + e − . Our result for this decay branching ratio is accurate at the 1% level and follows the pattern of previous theoretical results lying below the experimental value. Our method for combining lattice QCD with photon and lepton propagators extends techniques developed to calculate the hadronic light-bylight scattering contribution to g µ − 2 and holds promise for the eventual calculation of the two-photon-exchange contribution to the rare decay K L → µ + µ − .