Recent progress in applying lattice QCD to kaon physics

Standard lattice calculations in kaon physics are based on the evaluation of matrix elements of local operators between two single-hadron states or a single-hadron state and the vacuum. Recent progress in lattice QCD has gone beyond these standard observables. I will review the status and prospects of lattice kaon physics with an emphasis on non-leptonic $K\to\pi\pi$ decay and long-distance processes including $K^0$-$\overline{K^0}$ mixing and rare kaon decays.


Introduction
Since the discovery of kaons, the kaon physics plays a key role in the building of the Standard Model. The main mission for lattice QCD in kaon physics is to evaluate the low-energy hadronic effects to test the Standard Model parameters or to constrain on new physics. Lattice QCD has been successful for the calculations of the observables such as the pion and kaon decay constants f K ± π ± , the K → π ν semileptonic form factor f + (0) and the neutral kaon mixing parameter B K . We refer these observables as "standard". Their relevant hadronic matrix elements have only one local operator insertion. The initial and final states involve at most one stable hadron. Besides, the spatial momenta carried by initial/final-state particles are much smaller than the ultraviolet lattice cutoff 1 a, with a the lattice spacing. These standard observables can be computed with high statistical precision and controlled systematic errors using lattice QCD simulations.
Many interesting observables in kaon physics, however, are not "standard". One example is the calculation of K → ππ decay where the final state involves multiple hadrons. Another example is the evaluation of the long-distance contributions to flavor changing processes such as the calculation of the real and imaginary parts of K 0 -K 0 mixing amplitudes, which are related to the K L -K S mass difference ∆M K and the indirect CP violating parameter . Rare kaon decays including K → πνν and K → π + − also belong to this category. As these transitions proceed via the second-order weak interaction, the calculations would involve the construction of 4-point correlation function and the treatment of nonlocal matrix elements with two effective operator insertions. To tackle such quantities, one needs to develop new techniques.
In this report, I will first summarize the lattice QCD calculation of standard observables. They include f K ± π ± , f + (0), and inclusive τ → s decay. All these quantities are related to the determination V us f + (0) = 0.2165 (4), Lattice inputs of f + (0) and f K ± f π ± together with the experimental data give a precise determination of the CKM matrix elements V us = 0.2231 (7), V us V ud = 0.2313 (7).
In the Standard Model, the CKM matrix is unitary. Most stringent test of CKM unitarity is given by the first row condition V u 2 ≡ V ud 2 + V us 2 + V ub 2 = 1.
Using the results of V us and V ud given in Eq. (3), one finds that V u 2 = 0.9798(82), which has a 2.5 σ deviation from CKM unitarity. Currently the most precise determination of V ud = 0.97420 (21) is from superallowed nuclear β decay [4,5]. Using V us from K 3 decay and V ud from nuclear β decay sharpens the unitarity test with a much smaller uncertainty. However, the deviation is still around 2.4 σ, as shown in the second line of Eq. (5). If using V us V ud from leptonic decays and V ud from nuclear β decay, then the result confirms the CKM unitarity; see the third line of Eq. (5). The above tests of CKM unitarity are put together here for a comparison 0.9798(82), K 3 + leptonic decays, 0.9988(5), K 3 + nuclear β decay, 0.9998 (5), leptonic + nuclear β decay.
To clarify the 2.x σ deviation in the unitarity test, it is important to reduce the uncertainty from the lattice QCD determination of f + (0). One of the recent updates for f + (0) is from Fermilab Lattice-MILC collaboration. HISQ fermions on 2+1+1 flavor MILC configurations are used in the calculation and preliminary results are shown in Fig. 1. Compared to their report last year [6], more lattice ensembles are used in the analysis. Employing 4 ensembles at the physical pion mass and 2 ultra-fine lattice spacings allows them to reduce the statistical error to 0.14%. At 0.12 fm and m l m s = 0.1, they use three different volumes. Three volumes together with one-loop chiral perturbation theory (ChPT) [7] allow for a good estimate of the finite-volume effects. After chiral and continuum extrapolation, the total uncertainty is expected to be reduced to 0.2%, which is close to the current experimental uncertainty [2]. The calculation is performed at 5 lattice spacings 0.15, 0.12, 0.09, 0.06 and 0.042 fm, including 4 ensembles with physical pion mass. Open green symbols correspond to different volumes for a=0.12 fm and ml = 0.1 ms. The solid magenta line is the (preliminary) interpolation in the light-quark mass, keeping the strange-quark mass ms equal to its physical value, and turning off all discretization effects. The magenta diamond is the corresponding interpolation at the physical point. Data at the same light-quark mass but different lattice spacing are off-set horizontally.
Another lattice calculation is recently reported by JLQCD collaboration [8,9]. In their calculation, the chiral symmetry is exactly preserved by using the overlap quark action, which enables a direct comparison of the lattice data with ChPT and hence a determination of relevant low energy constants within NNLO ChPT. A reasonable agreement between lattice results for the slope d f + (q 2 ) dq 2 (at q 2 = 0) and experiment is observed, although the error is still large due to the high cost of the usage of overlap fermion.

τ inclusive decay and V us
The average of V us is summarized and updated on Spring 2017 by Heavy Flavor Averaging group (HFLAV) [10]; see the left panel of Fig. 2, where f + (0) and f K ± f π ± take the value from PDG 2016 [11]. The result from leptonic decays shows consistency with CKM unitarity, while the one using K 3 decays has a ∼ 2 σ deviation from CKM unitarity. The largest discrepancy happens for the case using the τ → s inclusive decay, where a 3.2 σ deviation from CKM unitarity is observed.
To explore the discrepancy, the main quantity of interest is the ratio of the decay rates where τ → s-hadrons ν τ indicates that in the decay the final-state hadrons contain net strange-ness. According to the optical theorem, the imaginary part of the hadronic vacuum polarization (HVP) functions can be related to the R-value through [12] dR where s is the invariant mass square of the final-state hadrons. S EW is a known short-distance electroweak correction [13]. Π (J) (s) are the HVP functions with the superscript (J) corresponding to angular momenta J = 0, 1. Once Im Π (J) (s) is known, Eq. (7) can be used to determine V us . Since Im Π (J) (s) is generically non-perturbative at small s, the conventional approach to determine Im Π (J) (s) is to use the dispersion relation [12] s0 where Im Π(s) on the left-hand side can be related to dR ds and V us , while the integral on the righthand side can be determined using QCD perturbation theory (pQCD) and operator product expansion (OPE). The parameter s 0 should be sufficiently large for a good convergence of pQCD and the validity of the OPE. W(s) is a weight function. If there is no pole inside the contour, then the integral along the branch cut is equal to the integral on the circle and then V us can be determined. A difficulty here is that the estimate of high-dimensional OPE terms relies on some assumptions and thus contains potentially large systematic effects. Using the conventional approach described above, it results in the low value of V us shown in the left panel of Fig. 2 [14]. An improvement is proposed by Ref. [15] to use different s 0 and weight functions W(s) and then study the dependence on s 0 and W(s). Through fit, not only V us , but also the OPE effective condensates are fit to experimental measurements (and also lattice QCD data). With this improvement, the 3.2 σ deviation is reduced to 1-2 σ level depending on using BaBar or 2014 HFAG result for These results are plotted on the right-panel of Fig. 2, denoted as "τ FB FESR, HLMZ17". Another new approach proposed by H. Ohki et. al. [16] is to let s 0 → ∞ and use the weight function containing the pole structures where the N different poles Q 2 k are spanned by a spacing ∆ = 0.2 (N − 1) GeV 2 , with a center point called C. Once W(s) is given, the contour of the integral is equal to the residues of the poles, which can be determined using lattice HVPs. Thus the value of V us can be determined accordingly. The strategy to choose Q 2 k is that it should not be too large to suppress the contribution from pQCD and  Figure 2. Left: HFLAV summary of Vus . For the τ → s inclusive decay, there is a 3.2 σ deviation from CKM unitarity. Right: New implementations for τ → s inclusive decay. The read circle data points, denoted as "τ FB FESR, HLMZ17", show the improvement by fitting the OPE effective condensates to the experimental measurements and lattice QCD data [15]. The green square data points show the improvement using W(s) in Eq. (10) and using lattice HVPs for the residues of the integral [16]. Both improvements shed the light on the resolution of the puzzle from the τ → s inclusive decay.
OPE at s > m 2 τ . It should not be too small to avoid large statistical error from lattice HVPs. The realistic calculation is performed using N f = 2 + 1 Möbius domain wall fermions at near-physical pion mass with the lattice spacings a −1 = 1.73 and 2.36 GeV and the lattice volume V = (5 fm) 3 . The corresponding results are summarized in the right panel of Fig. 2, denoted by the green square data points (The filled-square points are generated using τ → Kν τ data for the K pole, while the opensquare ones use the K µ2 data as input). Using different N and C, the lattice calculation shows the consistent results. Besides, all the data points are systematically larger than the conventional value of V us . At N = 4 and C = 0.7 GeV 2 , V us is obtained as (21), using τ → Kν τ input for the K pole 0.2245(16), using K µ2 input for the K pole.
Both improvements proposed by Refs. [15] and [16] shed the light on the resolution of the puzzle from the τ → s inclusive decay.

Neutral-kaon mixing parameter B K based on Standard Model and beyond
The parameter B K is related to the CP violating part of K 0 -K 0 mixing and thus short-distance dominated. Using OPE, the effective Hamiltonian H ∆S =2 eff can be written as a product of the Wilson coefficient C(µ) and the ∆S = 2 local operator Q ∆S =2 (µ) with G F the Fermi constant and M W the W-boson mass. It is a convention to use the parameter as a measure of indirect CP violation, with the initial state given by K L S particle and the final state having total isospin zero. The contribution from H ∆S =2 eff serves as a dominant contribution to Here the angle φ ≡ arctan(−2∆M K ∆Γ K ) ≈ 43.52(5) ○ [11], with ∆M K = M KL − M KS and ∆Γ K = Γ KL − Γ KS . M LD 00 indicates the long-distance contribution to and A 0 is the K 0 → (ππ) I=0 amplitude. Both M LD 00 and A 0 only make few-percent contributions to . The progress in lattice QCD calculation of these two quantities will be discussed later. Here we only focus on the M SD 00 . Within Standard Model, there is only one ∆S = 2 operator with V − A structure where the subscripts α and β denote the color indices. For beyond-Standard-Model theories, 4 other operators are possible The neutral-kaon mixing parameter B K and B i in the MS scheme are defined as where µ is the renormalization scale, f K the kaon decay constant and M K the kaon mass. Given the anomalous dimension γ(g), the renormalization group independent B parameterB K is related to B K (µ) by the formulâ For the Standard ModelB K , the lattice calculation has reached a precision of 1.3% for 2+1 flavor For beyond-Standard-Model B i (µ) at the MS scale µ = 3 GeV, the uncertainties of 2+1 flavor lattice results are about 2-5% [1] The results for B i (µ) from various groups are summarized by FLAG [1] on the left panel of Fig. 3. There are clear discrepancies in B 4 and B 5 from different groups. To resolve these discrepancies, RBC-UKQCD collaboration undertakes a study using both RI-MOM and RI-SMOM renormalization [17][18][19]. The calculation is performed using N f = 2 + 1 flavor domain wall fermion at M π ≈ 300 MeV and two lattice spacings a = 0.08 and 0.11 fm. The corresponding results are shown by the three data points below the legend of "RBC-UKQCD '16" on the right panel of Fig. 3, where the four red data points use RI-MOM scheme while the two green ones use RI-SMOM scheme and the orange one uses 1-loop lattice perturbation theory. Including the new RBC-UKQCD updates, all RI-MOM results are compatible among different groups. However, these results are systematically smaller than that from the RI-SMOM renormalization and the 1-loop lattice perturbation theory. Particularly for RI-SMOM calculation, both ( q, q) and (γ, γ)-schemes are used and consistent results (shown by the two green data points) are obtained after conversion to MS scheme. For B 5 , the situation is very similar. According to the study by Ref. [20], RI-SMOM renormalization is expected to have smaller infrared contamination than RI-MOM due to the usage of the non-exceptional momenta. The studies by [17][18][19] confirm such expectation and suggest that for B 4 and B 5 RI-SMOM renormalization should be used to fully control the infrared contamination. In Ref. [21] an update is reported on the RBC-UKQCD measurement of B i (µ), simulated using N f = 2 + 1 domain wall fermions at the physical quark masses.

K → ππ decay and direct CP violation
CP violation is first observed in neutral kaon decays. Under CP transform, the K 0 state is related to The CP eigenstate can be defined as the combination of K 0 and K 0 with K 0 + − the CP-even/odd state. The physical states observed in the experiment are the weak eigenstates K S and K L . K S decays into two pions and K L decays into three pions. By neglecting CP violation, K S is equal to CP-even state and K L equal to CP-odd state. In 1964, BNL discovered that K L is able to decay into two pions, indicating the violation of CP symmetry. This discovery leads to the Nobel prize in 1980.
Since K L and K S are not CP eigenstates, one can write them as a mixture of the CP eigenstates The parameter¯ is a measure of the strength of mixing. For K L → ππ decay, there are two contributions to the CP violation. The first part appears as the CP-even component of K L decays into two pions. This is called indirect CP violation and described by a parameter or in many cases written as K . receives its dominant contribution from¯ and a small contribution from A 0 due to its definition given in Eq. (13) The second contribution is from the CP-odd component of K L , which decays into two pions directly. This is called direct CP violation and denoted as ′ .
The experiments measure the decay amplitudes of K L → ππ and K S → ππ and use the ratio to determine the parameter and ′ . Using the experimental measurements of η +− and η 00 as input, PDG quotes [11] ≈ 1 3 is at the order of 10 −3 and ′ is even 1000 times smaller. Due to its small size, direct CP violation ′ is very sensitive to new physics.
For theoretical simplicity, it is convenient to study the decay amplitudes in the specific isospin channels, A 0 and A 2 , where δ I is the strong phase from ππ scattering. If CP symmetry were protected, then both the amplitudes A 2 and A 0 are real. To obtain the CP violation, one shall determine both real and imaginary part of A 2 and A 0 . The indirect CP violation only has a small dependence on A 0 as shown by Eq. (24). While for ′ , it is sensitive on both real and imaginary part of A 2 and A 0 The target for lattice QCD calculation is to determine A 2 and A 0 from first principles. The weak Hamiltonian for K → ππ decay is given by a series of ∆S = 1 local four-quark operators [22] Here τ = − VtdV * ts VudV * us = 1.543 + 0.635i is the ratio of CKM matrix elements. z i (µ) and y i (µ) are known perturbative Wilson coefficients that summarize the short-distance effects. The 10 local four-quark Electro-weak penguin Q 7 -Q 10 operators Q i can be matched to three types of diagrams in the full theory, shown in Fig. 4 with the notations "Current-current", "QCD penguin" and "Electro-weak penguin", respectively. The most recent updated results for the amplitude A 2 are given by RBC-UKQCD collaboration [23], where two ensembles are used, both at physical pion mass but with different lattice spacings. The parameters are given in Table 1. After continuum extrapolation the results for A 2 are given by which is calculated using Lüscher's formula [24] and consistent with the phenomenological curve from Ref. [25]. Table 1. Ensembles used in the recent lattice calculation of A2 by RBC-UKQCD collaboration [23]. In addition to the determination of A 2 and δ 2 , another outcome from Ref. [23] is the resolution of the puzzle of the ∆I = 1 2 rule. According to experimental measurement, the size of A 0 is about 22.5 times larger than that of A 2 . It is a more-than-half-century's puzzle since 1955 [26] on why the amplitudes in the different isospin channels are so much different. The Wilson coefficients only account for a factor of 2. The lattice calculation shows that Re[A 2 ] is dominated by diagrams C 1 and C 2 in the left panel of Fig. 5, where C 1 is color diagonal and C 2 color mixed. C 2 is 1 N suppressed relative to C 1 with C 2 equal to 1 3 of C 1 in leading order QCD perturbation theory. However, the lattice results in the right panel of Fig. 5 shows that C 2 is about −0.7 × C 1 , indicating very strong non-perturbative effects. As Re[A 2 ] is proportional to C 1 + C 2 , the observation that C 1 and C 2 have opposite signs leads to a significant cancellation between the two terms. While for Re[A 0 ], the opposite signs lead to an enhancement as Re[A 0 ] receives an important contribution from 2C 1 − C 2 . When considering the complete contribution to Re[A 0 ], including the disconnected diagrams, the size of Re[A 0 ] is more enhanced. In total, the hadronic matrix elements including the contributions from C 1 , C 2 and other diagrams would contribute another factor of ∼ 10. The cancellation between C 1 and C 2 is first observed in an earlier study [27] and is further confirmed by the latest calculation of A 2 [23]. So now the puzzle of ∆I = 1 2 rule is resolved from first principals. We have also seen a recent study of the ∆I = 1 2 rule with the scaling of the number of color [28].
The more demanding calculation is the K → ππ decay in the isopsin I = 0 channel. The latest calculation is performed at the physical kinematics M π = 143.1(2.0) MeV and M K = 490(2.2) MeV, using a 32 3 × 64 lattice volume and a lattice spacing a = 0.14 fm [29]. G-parity boundary condition is used and the lattice volume is chosen such that the kaon's mass is equal to the pion-pion's energy in the ground state.  [24], the I = 0 ππ scattering phase shift is found to be δ 0 = 23.8(4.9)(1.2) ○ , which is smaller than the value δ 0 = 38.0(1.3) ○ , obtained by combining experimental data with the Roy equations [30,31]. It remains a puzzle for the discrepancy and needs to be understood in the future study.
Using the lattice results for both A 0 and A 2 , the direct CP violation ′ can be determined: There is a 2.1 σ deviation from experimental value Re[ ′ ] = 1.66(23)×10 −3 [32]. As the uncertainties of the lattice results are larger than experimental measurement, to confirm whether new physics information can be found in the deviation, more accurate lattice calculations are required. It is reported by C. Kelly [33] that the statistics of the previous RBC-UKQCD calculation has been increased to 584 configurations. In the lattice calculation, the largest contribution to Im[A 0 ] comes from Q 6 operator. Fig. 6 shows the fit to obtain the matrix element ⟨ππ Q 6 K⟩. When the statistics increases from 216 to 584 configurations, uncertainty decreases as expected while the central values remain consistent. The aim of the RBC-UKQCD K → ππ program is to reduce the dominant statistical error for Re[ ′ ] in Eq. (33) by a factor of 2 within the next year.  Besides for the effort to increase the statistics, there are also efforts for improvements of various systematic effects. For example, the σ field starts to be added into the calculation to account for the σ → ππ effects in the I = 0 ππ scattering channel. In Ref. [34] N. H. Christ reports on including electromagnetism in K → ππ decay. The ∆I = 1 2 rule may make the effects of electromagnetism on A 2 ∼20 times larger than a naive O(α e ) estimate due to the mixing with A 0 . Such effects will become important if the target of the future calculations is to determine ′ with a precision of ∼ 10%. In Ref. [35] M. Bruno presents a non-perturbative calculation of Wilson coefficients even including the W-boson by using the technique of step scaling. Although the current availability of lattice spacings restricts the calculation to unphysically light W-bosons with M W ∼ 2 GeV, the calculation opens a new direction in the future to non-perturbatively determine the Wilson coefficients with controlled uncertainties.
In addition to the efforts from RBC-UKQCD collaboration to compute the K → ππ decay, N. Ishizuka et. al. are running a parallel program using the improved Wilson fermion action [36]. As a first step to verify the possibility of calculations with the Wilson fermion action, they consider the decay amplitudes at an unphysical quark mass M K ∼ 2M π . A large enhancement of the ratio A 0 A 2 is found at unphysical quark masses.

Long-distance contributions to flavor changing process: ∆M K and
Both the K L -K S mass difference ∆M K and indirect CP violating parameter are related to the mixing of the K 0 and K 0 . Such mixing is caused by the weak interaction as the strangeness differs by 2 in K 0 and K 0 . The time evolution of the K 0 -K 0 mixing system can be given by the equation where M is the mass matrix and Γ the decay width matrix. These 2 × 2 matrices are calculated to the 2 nd order of the weak interaction and given by where the indices i and j take the values 0 and0. H W is the ∆S = 1 weak effective Hamiltonian and P indicates that the principal part should be taken when an integral with a vanishing energy denominator is encountered. The mass matrix can be diagonalized. By neglecting the effects of CP violation, the mass difference ∆M K can be given by the real part of M¯0 0 through The parameter is related to the imaginary part of M¯0 0 and given explicitly in terms of the shortdistance and long-distance part of Im[M¯0 0 ] in Eq. (14). Both ∆M K and arise from an amplitude in which two W bosons and internal up-type quarks form a loop, shown by Fig. 7. The loop integral is proportional to the internal quark mass square m 2 q for q = u, c, t. As ∆M K is related to Re[M¯0 0 ], it is associated with the CP conserving part of K 0 -K 0 mixing amplitude. Although the top quark loop is enhanced by m 2 t , there is a significant suppression from the CKM factor λ t , where λ q = V qd V * qs . Due to the fact that Re[λ 2 c ] , the contributions to ∆M K are dominated by charm-charm quark loop. As it is sensitive to the charm quark mass, the K L -K S mass difference historically led to the predication of the charm quark fifty years ago [37][38][39]. For , it is related to the CP violating part of K 0 -K 0 mixing. The charm quark contribution is significantly suppressed as Im[λ 2 c ] ≪ Re[λ 2 c ]. In , the top-top, top-charm and charm-charm loops compete in size. As it contains important top-top loop contribution, is sensitive to the Standard Model parameter, λ t or V cb .
As a subsequent work of Refs. [41,42], a recent calculation of ∆M K is performed on a 2 + 1 flavor 32 3 × 64 Möbius domain wall lattice with the Iwasaki + DSDR gauge action. A near-physical pion mass M π = 170 MeV and the kaon mass M K = 492 MeV are used. Since the calculation is performed at a coarse lattice spacing with a −1 = 1.38 GeV, the charm quark mass m MS c (3 GeV) = 750 MeV is unphysically light. The calculation has included all the contractions from Type 1 to Type 4 shown in Fig. 8. Based on 120 configurations, the preliminary lattice result is given by ∆M K = 3.85(46) × 10 −12 MeV, which is consistent with the experimental value ∆M K = 3.483(6) × 10 −12 MeV [11]. However, since the calculation uses unphysical kinematics, this agreement could easily be fortuitous. Note that in the calculation of ∆M K , the loop integral involves double Glashow-Iliopoulos-Maiani (GIM) cancellation [38] and thus, there is no short-distance divergence. On the other hand, the double GIM subtraction makes ∆M K significantly rely on the charm quark mass. As a consequence, it is important to carry out the calculation at the physical charm quark mass. Figure 7. K 0 -K 0 mixing in the full theory. ∆MK is related with the CP conserving part of K 0 -K 0 mixing and thus long-distance dominated. The process is described by two ∆S = 1 operators. is related to the CP violating part of K 0 -K 0 mixing and thus short-distance dominated. The dominant contribution is described by a single ∆S = 2 operator and the relevant hadronic matrix element can be converted to BK. The remaining long-distance contribution below the scale of the charm quark mass has been calculated by Ref. [40]. A new RBC-UKQCD project, reported by C. Sachrajda [43], uses both physical pion and charm quark masses in the calculation. The computation of ∆M K is performed on a 64 3 × 128 lattice with the Iwasaki gauge action and Möbius domain wall fermions at an inverse lattice spacing of 2.359(7) GeV. Various techniques such as the use of all-to-all propagators and all mode averaging are used to reduce the statistical uncertainty. Based on 59 configurations, the preliminary result of ∆M K = 5.5(1.7) × 10 −12 MeV is consistent with the experimental value. The project has planned to collect 160 measurements in total.
The status and prospects of the determination of are updated by W. Lee in Refs. [44,45]. The estimate of is made using the FLAG value for B K , the angle-only-fit results for the Wolfenstein parameters and the CKM matrix element V cb from exclusive or inclusive decays. The preliminary Here the exclusive V cb is determined using the experimental measurements ofB → D * ν andB → D ν together with the lattice QCD calculation for the corresponding hadronic matrix elements [10,[46][47][48]. The inclusive V cb is determined using the inclusive decay processB → X c ν and QCD sum rules [49]. When using exclusive V cb as input, there is a 3.3 σ deviation between Standard Model value and experimental measurement exp = 2.228(11) × 10 −3 . Besides, V cb dominates the current 10% Standard Model uncertainty for . Therefore, it is important to have an accurate determination of V cb . On the other hand, it is also important to compute the long-distance contribution to precisely, whose size is expected to be a few percent but remains not well understood.
To calculate the long-distance contribution to , it is better to write the GIM cancellation by subtracting the charm quark propagator [40,41] q=u,c,t λ q p By doing so, the double GIM subtraction results in three terms in the effective Hamiltonian, with the coefficients λ 2 u , λ u λ t and λ 2 t , respectively. The λ 2 u term is irrelevant for . The λ 2 t term is purely short-distance dominated. Therefore the only interesting term for lattice QCD calculation is the λ u λ t term.
In the lattice QCD calculation of λ u λ t contribution, the top quark field shall be integrated out, leaving a QCD penguin operator, shown in Fig. 9. This QCD penguin operator can be neglected in the calculation of ∆M K as it carries a suppression factor of λ t λ u , but it is important for . The QCD penguin operator together with the current-current operator can form a new Type 5 diagram.
Without top quark in the lattice calculation, there is only one GIM subtraction and as a consequence the loop integral is logarithmic divergent. This divergence is cut off by an unphysical lattice scale, the inverse lattice spacing 1 a. One can define a bilocal operator in the RI-SMOM scheme by subtracting the unphysical short-distance contribution, and then match the bilocal operator in the RI-SMOM scheme to the one in the MS scheme using perturbation theory. More details on short-distance correction can be found in Refs. [40,50,51].
The calculation of is performed on a 24 3 × 64 lattice with domain wall fermion and Iwasaki gauge action [40]. The inverse lattice spacing is a −1 is 1.78 GeV. The pion mass is 339 MeV and the kaon mass 592 MeV. It uses 200 configurations and includes all Type 1-5 diagrams. In Table 2  Z-exchange Figure 10. Examples of W-W and Z-exchange diagrams for K + → π + νν decay.
scale µ RI ranging from 1.54 to 2.56 GeV. The µ RI dependence is accounted for as a systematic uncertainty. At µ RI =2.11 GeV, the long-distance contribution to is about 5% when compared to the experimental value exp = 2.228(11) × 10 −3 . To accurately estimate the long-distance contribution, the calculation needs to be performed at the physical kinematics.

Long-distance contributions to flavor changing process: rare kaon decays
Rare kaon decays have attracted increasing interest during the past few decades. As flavor changing neutral current processes, these decays are highly suppressed in the Standard Model and thus provide ideal probes for the observation of new physics effects. In this review, I will discuss the lattice QCD calculations of two classes of rare kaon decays: K → πνν and K → π + − [50][51][52][53][54][55][56][57].
The K + → π + νν decay is interesting because it receives the largest contribution from top quark loop and thus theoretically very clean. The required hadronic matrix elements can be obtained from leading order semi-leptonic K decays, such as K + → π + eν, via isospin rotation. The remaining long-distance contributions below the charm scale are expected to be a few percent. Though small, by including the long-distance contribution estimated from Ref. [58], the branching ratio Br(K + → π + νν) is enhanced by 6%, which is comparable to the 8% total Standard Model uncertainty [59]. The current known branching-ratio measurement [60] Br is a combined result based on the 7 events collected by BNL E787 [61][62][63][64] and its successor E949 [60,65]. Its central value is almost twice of the Standard Model prediction [59] Br(K + → π + νν) SM = 9.11 ± 0.72 × 10 −11 , but with a 60-70% uncertainty it is still consistent with Standard Model.
The new experiment, NA62 in CERN [66], aims at an observation of O(100) events and a 10%precision measurement of Br(K + → π + νν). The status reported at the Flavor physics and CP violation workshop (FPCP 2017) is that the detector installation is completed in September 2016. 5% of the 2016 data has been analyzed but no event is found yet. If using full 2016 data, then O(1) events are expected to be found. Considering the fact that the Standard Model predictions will be confronted with  Figure 11. Low-lying intermediate states contributing to K + → π + νν. As these states are related to exponentially growing unphysical contributions and potentially large finite volume effects, one shall calculate the hadronic matrix elements for these low-lying intermediate states from the relevant 2-point and 3-point functions.
the new experiment soon, a lattice QCD calculation of the long-distance contribution to K + → π + νν is timely.
There are two classes of diagrams, which contribute to K + → π + νν decays, called as W-W and Z-exchange diagrams. In the W-W diagrams the second-order weak transition proceeds through the exchange of two W-bosons, while for the Z-exchange diagrams the decay occurs through the exchange of one W-boson and one Z-boson. Examples of both classes of diagrams are illustrated in Fig. 10.
In a lattice QCD calculation, the W and Z-boson have been integrated out, leaving two effective four-fermion local operators. The matrix element of the time-integrated bilocal operator is evaluated in Euclidean space. This matrix element can be related to the second-order amplitude of interest if a sum over intermediate states is inserted and the integration over Euclidean time performed: where H A B (t) stands for the two four-fermion operators, with the spatial variables integrated over space. The unphysical e (EK −En)T terms in the second line of this equation vanish for large T for intermediate states more energetic than the kaon. However, these terms grow exponentially with increasing integration range if E n < E K and must be removed from lattice calculation. When the intermediate state involves multiple particles, the branch-cut integral in the infinite volume is replaced by a discrete state summation in the finite volume. It could cause potentially large finite-volume effects when E n → M K , which need to be corrected following Ref. [67]. To deal with the exponentially growing terms and finite volume effects, the matrix elements for the lowing-lying intermediate states shall be calculated. These states include the leptonic + ν, semileptonic π 0 + ν, single pion and isospin I = 2 π + π 0 scattering state and are summarized by Fig. 11. So the study of the long-distance contribution to K + → π + νν decay does not only involve the calculation of 4-point function, but also includes the calculation of all relevant 2-point and 3-point functions for low-lying intermediate states.
As the top quark contribution to the decay is completely short-distance dominated, one only needs to focus on the charm quark contribution. The first calculation is performed using the 16 3 × 32, N f = 2 + 1 flavor, domain wall fermion ensemble, with a −1 = 1.729(28) GeV [51]. This ensemble has pion and kaon masses of M π ∼ 421 MeV and M K ∼ 563 MeV. The MS charm quark mass is m MS c (2 GeV) ∼ 863 MeV. Both W-W and Z-exchange diagrams are logarithmically divergent and cutoff by unphysical scale 1 a. Similar as the computation of , the short-distance correction needs to be performed here [50]. The lattice results are shown in Fig. 12. Here P c gives the complete charm quark contribution to the K + → π + νν decay. The results from the W-W and Z-exchange diagrams, and their total, are shown in the left, center and right panels. The gray bands show the bilocal matrix element including the unphysical lattice artifacts. The red circles indicate the RI-renormalized, bilocal contribution. The blue diamonds give the total charm contribution P c , while the green squares show the difference between the lattice and perturbative results, P c − P PT c . The results from the exploratory lattice calculation with unphysical charm, down and up quark masses are: L (1, 1, 1) Figure 13. Dependence of the form factor for the decay K + → π + + − upon z = q 2 M 2 K . The lattice data is fit to a linear form V+(z) = a+ + b+z.
The small size of P c − P PT c results from a large cancellation between the W-W and Z-exchange amplitudes. It is important to determine whether such a large cancellation persists for physical quark masses.
Different from K + → π + νν decay, the CP conserving decays K + → π + + − and K S → π 0 + − receive the dominated long-distance contribution from γ-exchange diagram. Although the loop integral in the γ-exchange diagram is quadratically ultraviolet divergent by power counting, the electromagnetic gauge invariance reduces the divergence to be logarithmic. The GIM cancellation further reduces the logarithmic divergence to be ultraviolet finite.
In the γ-exchange process, the hadronic part of amplitudes for K + and K S decays can be written in terms of electromagnetic transition form factor V + S (z) via [68] with p K π the kaon/pion momentum, q = p K − p π , z = q 2 M 2 K and r π = M π M K . The target for lattice QCD calculation is to extract V + S (z) from the bilocal hadronic matrix elements by building the relevant 4-point correlation functions. The strategy adopted in Ref. [53] is to use conserved vector current to protect the electromagnetic gauge invariance and use the charm quark as an active quark flavor to maintain GIM cancellation. The first exploratory calculation of K + → π + + − decay [57] is performed using a 24 3 × 64 lattice with domain wall fermion and Iwasaki gauge action. The inverse lattice spacing is 1 a = 1.78 GeV. The calculation uses a pion mass of M π ∼ 430 MeV, a kaon mass of M K ∼ 625 MeV and a MS charm quark mass m MS c (2 GeV) ∼ 533 MeV. Three momentum transfers are used in the calculation and a linear fit form V + (z) = a + + b + z is used to determine the momentum dependence of the form factor. The lattice data points for V + (z) together with the fit curve are shown in Fig. 13. Using 128 configurations, the results for a + and b + yield a + = 1.6(7), b + = 0.7 (8).
The phenomenological study [69] decomposes the form factor into a linear form plus the unitarity ππ-loop correction V ππ Here V ππ + (z) is determined using chiral perturbation theory together with some model assumptions such as vector meson dominance model. The experimental measurements of the branching ratio together with the fit form (45) produce the much more precise results for a + and b + [70,71] ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ a + = −0.58 (2), b + = −0.78 (7), fit from K + → π + e + e − spectra, a + = −0.58 (4) where the first line uses the experimental measurement of K + → π + e + e − spectra and the second line uses K + → π + µ + µ − data. Note that these results carry the opposite signs to the lattice results in Eq. (44). Since the lattice calculation is performed at unphysical quark masses, it does not make much sense to compare these results in Eqs. (44) and (46). On the other hand, as the experimental data only yield the square of the form factor and does not tell the sign for V + (z), the signs for a + and b + are completely determined by the input of V ππ + (z). However, it is found that the polynomial contribution (linear in z) dominates over the unitarity loop correction. For K + → π + µ + µ − decay, the fit forms with and without V ππ + correction produce almost the same fit curves. For K + → π + e + e − decay, the fit curves differ at small z, where the experimental data is not available [68]. Therefore it is questionable to use V ππ + (z) to determine the signs for a + and b + . It is important to perform a lattice QCD calculation at the physical quark mass and examine the phenomenological fit ansatz (45) and confirm the sign for a + and b + .
In order to perform the calculation at the physical point, the physical pion mass would require the large lattice volume and physical charm quark mass would require for the ultra-fine lattice spacing to control both finite-volume effects and lattice artifacts. Thus it is very high demanding on computer resources. One solution is to improve quark action to reduce the lattice artifacts for the charm quark. In Ref. [72] M. Tomii performs an exploratory study of dispersion relation and unphysical pole for Mobius domain wall fermion and seeks for a way to improve the action. Another solution is to integrate out of charm quark field using perturbation theory. In this case, lattice QCD calculation only requires the physical pion mass and a rather coarse lattice spacing. This would save the computer resources quite significantly. But a drawback is that since there is no GIM cancellation, the internal up quark loop will be logarithmically divergent. In Ref. [73] A. Lawson discusses on the renormalization to treat with the short-distance divergence in three-flavor theory.
Besides for the CP conserving K → π + − decay, it is also interesting to study the CP violating K L decays. The K L decay amplitudes receive three major contributions: 1) a short-distance dominated direct CP violation, 2) a long-distance dominated, indirect CP violating contribution through K L → K + → π 0 + − , 3) a CP conserving component which proceeds through two-photon exchange. Total CP violating contributions to K L decay branching ratios, including 1), 2) and their interference, are given by [74,75] Br(K L → π 0 e + e − ) CPV = 10 −12 × 15.7 a S 2 ± 6.2 a S Im λ t 10 −4 + 2.4 Im λ t 10 −4 , Br(K L → π 0 µ + µ − ) CPV = 10 −12 × 3.7 a S 2 ± 1.6 a S Im λ t 10 −4 + 1.0 where the λ t ≈ 1.35 × 10 −4 . The parameter a S is given by the K S transition form factor at zero momentum transfer, namely a S = V S (0), and it is a quantity of size O(1). The ± sign arises because only the magnitude of a S is determined from experiment. Therefore even a determination of the sign of a S from lattice QCD is desirable.

Conclusion
The worldwide lattice QCD community has developed a successful kaon physics program. It inspires the consideration of constructing a CKM unitarity triangle purely from kaon physics [76].
For standard quantities such as f K ± f π ± , f + (0) andB K , they are computed with a precision of ∼ 1 percent or even much better, as shown in Table 3. In these cases lattice QCD calculations play important roles in precision flavor physics. With the development of the lattice QCD techniques, it is also the time to explore the non-standard quantities. Here I report the recent progress of the calculations on K → ππ decay, long-distance contributions to ∆M K and as well as rare kaon decays. Lattice QCD is now capable of first-principals calculation of these non-standard quantities. Some of the calculations are even performed at the physical kinematics. We can foresee that with the new techniques and new generation of super-computers today's non-standard observables will become standard in the near future. Table 3. Summary of the FLAG average of f K ± f π ± , f+(0) andBK.