Lattice QCD study of the rare kaon decay $K^+\to\pi^+\nu\bar{\nu}$ at a near-physical pion mass

The rare kaon decay $K^+\to\pi^+\nu\bar{\nu}$ is an ideal process in which to search for signs of new physics and is the primary goal of the NA62 experiment at CERN. In this paper we report on a lattice QCD calculation of the long-distance contribution to the $K^+\to\pi^+\nu\bar{\nu}$ decay amplitude at the near-physical pion mass $m_\pi=170$ MeV. The calculations are however, performed on a coarse lattice and hence with a lighter charm quark mass ($m_c^{\bar{\mathrm{MS}}}(\mbox{3 GeV})=750$ MeV) than the physical one. The main aims of this study are two-fold. Firstly we study the momentum dependence of the amplitude and conclude that it is very mild so that a computation at physical masses even at a single kinematic point would provide a good estimate of the long-distance contribution to the decay rate. Secondly we compute the contribution to the branching ratio from the two-pion intermediate state whose energy is below the kaon mass and find that it is less than 1% after its exponentially growing unphysical contribution has been removed and that the corresponding non-exponential finite-volume effects are negligibly small.


I. INTRODUCTION
The rare kaon decays K → πνν have attracted increasing interest during the past few decades. As flavor changing neutral current (FCNC) processes, these decays are highly suppressed in the standard model (SM) and thus provide ideal probes for the observation of new physics. Moreover, these decays are short-distance (SD) dominated and are therefore theoretically clean. In this paper we denote contributions from distances greater than O(1/m c ), where m c is the mass of the charm quark, as long-distance (LD) contributions. Such LD contributions can safely be neglected in the CP-violating K L → π 0 νν decays and are expected to be small, perhaps of O(5 -10%), in K + → π + νν decays. The study reported in this paper is the next step in our project to compute the LD contributions non-perturbatively using lattice QCD and to obtain the branching ratio for K + → π + νν decays with a QCD precision of O(1%) (uncertainties in the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements need to be reduced in parallel with the reduction of QCD errors). The usual approach simply uses perturbation theory to integrate out the charm quark. An estimate of the missing LD effects in this approach is based on chiral perturbation theory and suggests a (6 ± 3)% enhancement of the branching ratio Br(K + → π + νν) [1], which is comparable to the 8% total SM error [2]. We will test this result using lattice QCD.
The KOTO experiment at J-PARC, is designed to search for K L → π 0 νν decays [10]. It has observed one candidate event based on the first 100 hours of physics running in 2013 and set an upper limit of 5.1 × 10 −8 for the branching ratio at 90% confidence level [11]. For a recent description of the status of this experiment see Ref. [12]. Since the LD contributions to K L → π 0 νν decays are negligible, we do not discuss these decays further in this paper. This paper is the latest in a series of lattice QCD studies of the rare kaon decays K → π + − and K → πνν, in which we have developed the theoretical framework and performed exploratory numerical calculations [13][14][15][16][17][18][19][20][21][22][23][24]. Here we focus on the LD contributions to the K + → π + νν decay amplitude. Since we have two neutrinos in the final state and the strangeness changes by one unit, this is a doubly weak process and the long-distance contributions arise either from the exchange of two W-bosons (W -W diagrams) or the exchange of one Z-boson and one W -boson (Z-exchange diagrams) [17]. When evaluating the LD contribution of both types of diagram the W and Z exchanges are replaced by four-quark operators. However, the integral over the locations of these four-quark operators is logarithmically divergent when their locations coincide, requiring the introduction of an additional counter term. In the full SM these divergences are cut off by the mass of the W -boson. To obtain the physical amplitude, we must consistently combine the long distance contributions determined from lattice QCD with the short distance contributions which can be reliably determined using perturbation theory [25][26][27][28][29][30]. An exploratory calculation of the K + → π + νν decay amplitude has been performed in Refs. [23,24] with meson masses m π 420 MeV, m K 563 MeV and a charm-quark mass m MS c (2 GeV) 860 MeV. Because of the small phase space available with these masses, the allowed momenta for the final-state particles are constrained to lie in a narrow region and provide little information on the momentum dependence of the decay amplitude. For this reason, we computed the amplitude in Refs. [23,24] at a single choice of momenta. The momentum dependence of the decay amplitude was therefore unresolved and is the focus of the study reported here.
In this paper we calculate the amplitude at the near-physical pion mass m π 170 MeV and kaon mass m K 493 MeV and study the momentum dependence of the K → πνν decay amplitude. In addition since now, in contrast to the earlier studies in Refs. [23,24], the mass of the kaon is above the two-pion threshold, we are able to study the contribution of the ππ intermediate state to the decay amplitude together with the associated finite-volume effects.
The conclusions are presented in Sec. V but we anticipate them briefly here to say that we find that the momentum dependence of the LD contributions is very mild and that the contribution of the lowest-energy finite-volume two-pion intermediate state to the branching ratio is at the sub-percent level.
The plan for the remainder of this paper is as follows. In Section II we present the details of our computation. This is followed by the results for the momentum dependence in Sec. III.
The contribution of the lowest-energy, finite-volume two-pion intermediate state, computed at the maximum value of s = (p K − p π ) 2 , is discussed in Sec. IV. Finally in Sec. V we present our conclusions.

II. DETAILS OF THE LATTICE CALCULATION
In this study we use the 2 + 1 flavor domain wall fermion configurations with lattice size The decay amplitude for the rare kaon decay K + (p K ) → π + (p π )ν(p ν )ν(pν) can be written as the product of a scalar amplitude F (s, ∆) and a spinor productū(p ν ) / p K (1 − γ 5 )v(pν) [17]: Here we use Euclidean four-momenta with imaginary time components for on-shell particles, e.g. p π = (iE π , p π ) for the pion. The Lorentz invariant s ≡ −(p K − p π ) 2 is the square of the invariant mass of the νν pair and takes values in the range s ∈ [0, s max ] where We denote by ∆ max = m 2 K − m 2 π the absolute maximum of ∆; this can be reached for s = 0. One can also write ∆ as ∆ = ∆ max (s) cos θ, where θ indicates the angle between p π and p ν in the rest frame of the νν pair. For a particular K + → π + νν event the quantity ∆ cannot be determined experimentally and so a suitable phase-space integral must be performed to obtain either the full rate, or the differential distribution with respect to s (the variable s depends on only the momenta of the charged mesons and is therefore measurable).
Note that the square of the spinor product |ū( , so that at the edge of the allowed momentum region, where ∆ = ∆ max (s), the decay amplitude vanishes. We are interested therefore in computing the amplitude at momenta away from this kinematic edge. In our calculation, we compute the amplitude at the following four pairs of (s, ∆): The scalar amplitude F (s, ∆) is computed for both the W -W and Z-exchange diagrams. We denote the contribution from the W -W diagrams by F W W (s, ∆); as indicated this depends on both s and ∆. We denote the contribution from the Z-exchange diagrams by F Z + (s) which, again as indicated, depends on only the single variable s since the neutrino and antineutrino are emitted from a single point. The quantity F Z + (s) is analogous to the corresponding form factor in K 3 decays, but with the charged weak current replaced by the neutral one.
For massless neutrinos, the second form factor F Z − (s) does not contribute to the amplitude. We calculate the Z-exchange diagram at one additional kinematic point (s, ∆) = (s max , 0).
Although at this point F Z + (s) cannot be determined directly, we can obtain the amplitude For a complete explanation on how we compute the scalar amplitudes F W W (s, ∆), F Z ± (s) and F Z 0 (s), we refer the readers to Refs. [17,23,24].

III. MOMENTUM DEPENDENCE OF THE AMPLITUDE
For the Z-exchange diagrams, we calculate the scalar amplitudes F Z +,−,0 (s) at s = 0, s max /3 and s max /2 and present the results in Table I. For s = 0, we have F Z + (s) = F Z 0 (s). For s = s max /3 and s max /2, we find that the values of F Z 0 (s) are a little smaller than those of F Z + (s). However, within the statistic uncertainties, F Z 0 (s) and F Z + (s) are consistent, suggesting that F Z 0 (s) is a good approximation for F Z + (s). In addition we evaluate the scalar amplitude F Z 0 (s) at s = s max . At this kinematic point the pion and kaon are both at rest and so we can only determine a single form factor, F Z 0 (s max ). The main motivation for the computation of F Z 0 (s max ) is to provide an estimate of the contribution of the disconnected diagrams. At the other kinematic points we use twisted boundary conditions to obtain results for a range of momenta, all with s ≥ 0, whereas at s = s max we use periodic boundary conditions. With twisted boundary conditions for the valence down quark, the neutrinos carry momenta which are not integer multiples of 2π/L, and such momenta cannot be transferred by the gluons which satisfy periodic boundary conditions (see Fig. 1). From the results in Table I we conclude that the disconnected diagrams make only a small contribution to F Z 0 (s max ) and since we expect F Z + (s) F Z 0 (s) and have found the momentum dependence of F + (s) to be very mild (see Fig. 2 and the discussion below), this suggests that the contribution of the disconnected diagrams to the K + → π + νν decay amplitude is much less than 1%.
The momentum dependence of the form factors is plotted in Fig. 2. In the upper panel, we show the form factor f + (s) for the K + → π 0 + ν (K 3 ) decay. In the middle panel, we Ref. [24]): Here we use a superscript (Z) to indicate that this is the contribution from the Z-exchange diagrams. In order to obtain the full contribution of the Z-exchange diagrams to F Z + (s), we need to subtract the short-distance counter-term from C lat 1 Q 1 +C lat 2 Q 2 , the combination given in the last column of Table I, following the procedure described in detail in Refs. [17,24]. This short-distance contribution to F Z + (s), which is difficult to compute on such a coarse lattice, is proportional to f + (s) and hence appears as a constant term in P (Z) c (s). Since the focus of this paper is to explore the s-dependence of P    For the W -W diagrams, we calculate the scalar amplitude F W W for the momentum pairs (s, ∆) given in Eq. (2). This scalar amplitude is divided into two parts corresponding to the different contractions labelled as Type 1 and Type 2 in Fig. 3. The W -W diagrams contain charged-lepton propagators and therefore in Table II we present the contributions to F W W from each of three lepton flavors, = e, µ, τ , separately.
Combining the contributions from Type 1 and Type 2 diagrams and from the three lepton flavors, we obtain the results for F W W . In the upper panel of Fig. 4, we show the s and ∆ dependence of F W W (s, ∆). In the lower panel, we present the results for P  is obtained from F W W using We parametrize the momentum dependence of P (W W ) c (s, ∆) using the linear function and from the correlated fit to the lattice data, we find b As estimated in Ref. [24], using the linear parametrization of P c (s, ∆), the branching ratio of K + → π + νν is proportional to where b ∆ = b   The observation that the momentum dependence is so mild provides a useful guide for our future lattice computations of the K + → π + νν decay amplitude. To perform the calculation at physical kinematics is very challenging: on the one hand one needs a large volume in order to accommodate a pion with a mass of about 140 MeV and on the other hand one needs a fine lattice to avoid lattice artifacts from the physical mass of the charm quark. We would require very significant computational resources in order to perform the calculations for a wide range of values of (s, ∆). The study reported here suggests that the s-and ∆-dependence of the scalar amplitude has only a minimal effect on the branching ratio. Thus even computing the K + → π + νν decay amplitude at a single kinematic point (s, ∆), we can obtain a good estimate of the LD contribution to the decay rate.

IV. CONTRIBUTION OF THE TWO-PION INTERMEDIATE STATE
Two of the less familiar obstacles that must be overcome when computing the K + → π + νν decay amplitude using lattice QCD are the potentially large finite-volume distortions of the two-pion intermediate state energy spectrum at low energies and the unphysical terms which grow exponentially as the Euclidean time separation between the initial-and final-state interpolating operators is increased. Such terms which arise from intermediate two-pion states with energy E ππ below M K . In this section we discuss the corrections that we make for these effects when s = s max . Thus, for K + → (π + π 0 ) I=2 → π + νν we need to study contributions from the |π + π 0 intermediate state with I = 2. We first determine its energy.
In the top panel of Fig. 5 we show the ratio as a function of t. Here C ππ (t) is a two-point correlation function for the (ππ) I=2 interpolating operator in the center-of-mass frame and C π (t) is a correlation function for the single π interpolating operator at rest (see Ref. [24] for the choice of interpolating operators). The subtraction in C ππ (t + 1) − C ππ (t) is introduced to remove the constant term arising from the around-the-world effect and the corresponding difference C 2 π (t + 1) − C 2 π (t) is used in the denominator to insure that the ratio R(t) in Eq. (8) is proportional to exp(−∆E). Thus, by studying the t-dependence of R(t) as explained in Ref. [32], we determine the energy difference ∆E = E ππ − 2m π = 0.001687(41) (note that the pion mass in lattice units is 0.12535 (27)). Using Lüscher's finite-volume formula, we find the (ππ) I=2 scattering length to be given by m π a I=2 ππ = −0.0659 (16), in good agreement with the prediction from leading- In the next step, we calculate the 0 ππ|H W |K and π|J µ |ππ 0 matrix elements, where the subscript 0 indicates the finite-volume two-pion ground state, and construct the two ratios: where the normalisation factors are N π = 0 |π(0)|π( 0 ) , N K = 0 |K(0)|K( 0 ) and N ππ = 0 |ππ(0)|ππ 0 . We present R ππ→π (t) and R K→ππ (t) as functions of t in the central and lower panels of Fig. 5 respectively. Performing fits for R ππ→π (t) and R K→ππ (t) in the plateau regions, as indicated by the gray bands in Fig. 5, we determine the matrix elements to be 0 ππ|H W |K = −i 9.653(25) × 10 −5 and π|J 4 |ππ 0 = i 3.1910(60). The value of π|J 4 |ππ 0 is close to the factorization approximation π|J 4 |ππ 0 ≈ π|π 0|J 4 |π = i 3.2149(62). Using these matrix elements as inputs, we find that the ππ contribution to F Z 0 (s max ) is F Z,ππ 0 (s max ) = −1.536(5)×10 −3 , which is about 7.5% of the total scalar amplitude F Z 0 (s max ). Such effects are not small and result in a sizeable exponentially growing contamination. In Fig. 6, we show the lattice results with and without correcting for the exponentially growing contamination from the ππ intermediate state. The corresponding systematic effects are significant when compared to the statistical ones. We thus include these corrections in our analysis when producing the results shown in Table I.
In addition to removing the exponentially growing contamination, we also evaluate and correct for the corresponding finite-volume effects due to the ππ intermediate state. As shown in Ref. [24,33], the finite-volume correction which must be subtracted from the finitevolume amplitude to obtain the infinite-volume result is where φ is a kinematic function of E ππ and δ(E ππ ) is the I = 2 s-wave ππ phase shift. Note that in Eq. (10) all quantities are to be evaluated at E ππ = M K and the two-pion matrix elements must be evaluated in a spatial volume that has been chosen so that the finitevolume state |ππ has the energy of the kaon. We use the approximation δ ≈ k a ππ where k = E 2 ππ /4 − m 2 π and dδ/dE ππ ≈ a ππ Eππ 4k to estimate ∆ FV and find ∆ F V = 4.28(7) × 10 −4 , corresponding to a correction to P c of about 4.8 × 10 −4 . Here we approximate the energyconserving matrix elements needed in Eq. (10) by those determined above when E ππ = 0.25239(60). Since the expected value of the full charm-quark contribution to P c is about 0.4 [2], this finite-volume correction has a negligible, per-mille effect on the branching ratio (see Eq. (97) and the following discussion in Ref. [24]).

V. SUMMARY AND CONCLUSIONS
In this paper we report on a study of the long-distance contribution to the amplitude of the rare kaon decay K + → π + νν at the near physical pion mass of m π = 170 MeV.
This study is an intermediate step in our long-term project to evaluate the rate for this decay at physical kinematics and with sub-percent precision. Having previously developed the necessary theoretical framework for the evaluation of the long-distance contributions in Ref. [17] and performed a complete computation at unphysical quark masses and for a single choice of meson and neutrino momenta in Refs. [23,24], in this paper we investigate two areas to guide us and to provide experience for our future calculations at physical masses.
These are: 1. The momentum dependence of the amplitude. Since varying the momenta of the mesons and neutrinos requires very significant computational resources, it was instructive to gain some understanding of the behaviour of the amplitude as a function of the momenta. Given the limited kinematic range available, we expect the momentum dependence to be very mild and in Sec. III we explicitly confirm this and conclude that a computation at physical masses even at a single kinematic point (s, ∆) would provide a good estimate of the long-distance contribution to the decay rate. Of course, eventually as computing resources allow, this conclusion can be tested also at physical quark masses.
An additional conclusion from this study is that the contribution from disconnected diagrams appears to be of the order of 2% of that coming from the connected part. This was the case at s = s max , which is the only point at which disconnected diagrams could be evaluated, since it was at this point that twisted boundary conditions for the quarks were not used. M K , K + → (π + π 0 ) I=2 → π + νν, from the Z-exchange diagrams. Two consequence of this in the evaluation of long-distance contributions in a finite-volume Euclidean-space calculation are: (i) the presence of unphysical terms which grow exponentially with the range of the timeseparation of the two weak operators and which must be subtracted; (ii) finite-volume effects which are not exponentially small in the volume and which must be removed.
The procedures to remove the exponentially growing terms and the power-like finite-volume corrections are well understood. In Sec. IV we perform the subtraction of the exponentially growing terms and find that the contribution from the two-pion intermediate state is a little less than 1% of the full decay amplitude. We calculate the corresponding FV effects and find them to be negligibly small.
As the final step in this K + → π + νν project we have generated data on a 64 3 × 128 lattice with Möbius domain wall fermions and the Iwasaki gauge action at an inverse lattice spacing of a −1 = 2.359(7) GeV, pion mass m π = 135.9(3) MeV and kaon mass m K = 496.9(7) MeV [34]. We use a series of charm quark masses surrounding the physical value. The use of this ensemble will remove the significant systematic effect arising from the lighter than-physical charm quark mass of the present study, resulting in a physical rare kaon decay amplitude with the main systematic effects under control. The resources required for the 64 3 × 128 calculation are relatively demanding making the determination of the amplitude at several kinematic points computationally expensive. The observation of the very mild momentum dependence of the decay amplitude observed in this study implies that we can obtain the rate with good precision, even if we perform a computation at a single kinematic point of (s, ∆).