# Lattice QCD calculation of the electroweak box diagrams for the kaon semileptonic decays

Peng-Xiang Ma, Xu Feng, Mikhail Gorchtein, Lu-Chang Jin, Chien-Yeah Seng

LLattice QCD calculation of the electroweak box diagrams for thekaon semileptonic decays

Peng-Xiang Ma, Xu Feng,

1, 2, 3, ∗ Mikhail Gorchtein,

4, 5, 6, † Lu-Chang Jin,

7, 8, ‡ and Chien-Yeah Seng § School of Physics and State Key Laboratory of Nuclear Physicsand Technology, Peking University, Beijing 100871, China Collaborative Innovation Center of Quantum Matter, Beijing 100871, China Center for High Energy Physics, Peking University, Beijing 100871, China Helmholtz Institute Mainz, Mainz, Germany GSI Helmholtzzentrum f¨ur Schwerionenforschung, Darmstadt, Germany Johannes Gutenberg University, Mainz, Germany Department of Physics, University of Connecticut, Storrs, CT 06269, USA RIKEN-BNL Research Center, Brookhaven NationalLaboratory, Building 510, Upton, NY 11973 Helmholtz-Institut f¨ur Strahlen- und Kernphysik and Bethe Center for Theoretical Physics,Universit¨at Bonn, 53115 Bonn, Germany (Dated: February 25, 2021)

Abstract

We present a lattice QCD calculation of the axial γW -box diagrams relevant for the kaon semilep-tonic decays. We utilize a recently proposed method, which connects the electroweak radiativecorrections in Sirlin’s representation to that in chiral perturbation theory. It allows us to use theaxial γW -box correction in the SU(3) limit to obtain the low energy constants for chiral pertur-bation theory. From ﬁrst principles our results conﬁrm the previously used low energy constantsprovided by the minimal resonance model with a signiﬁcant reduction in uncertainties. ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] a r X i v : . [ h e p - l a t ] F e b . INTRODUCTION In the Standard Model, the Cabibbo-Kobayashi-Maskawa (CKM) matrix is a three-generation quark mixing matrix which describes how the strength of the ﬂavour-changingweak interaction in the leptonic sector is distributed among the three quark generations. Theprecise determination of the CKM matrix elements is of vital importance in the stringenttest of CKM unitarity and search of new physics beyond the Standard Model. As quotedin the 2020 Review by the Particle Data Group [1], there exists a ∼ ∣ V ud ∣ + ∣ V us ∣ + ∣ V ub ∣ = . ( ) V ud ( ) V us . (1)Here ∣ V ub ∣ ≈ . × − is negligibly small and thus only ∣ V ud ∣ and ∣ V us ∣ play a role in theunitarity test.The most precise determination of ∣ V ud ∣ = . ( ) exp + nucl ( ) RC quoted in the 2020PDG review [1] stems from the superallowed nuclear beta decays [2, 3], with the ﬁrst un-certainty arising from the experimental measurements and nuclear physics corrections andthe second one from the electroweak radiative corrections (RCs) . It is the update of theRCs from a dispersive analysis [4, 6], which makes the value of ∣ V ud ∣ about 2 σ smaller thanthat in the 2018 PDG review [7]. Very recently, the RCs to the π (cid:96) decays were calculatedusing lattice QCD with the focus on the so-called axial γW -box diagrams [8]. It allowed fora signiﬁcant reduction of the hadronic uncertainty in the RCs, and provided an independentcross-check of the dispersion relation analysis of the neutron RCs [9]. In the future a directlattice QCD calculation of the RCs to the neutron decay could help to further improve thedetermination of ∣ V ud ∣ [10].The ∣ V us ∣ can be determined from kaon, hyperon or tau decays, with kaon decays providingthe best precision. Leptonic decays K → µν (denoted by K µ ) combined with π → µν give access to the ratio ∣ V us / V ud ∣ , whereas semileptonic decays K → π(cid:96)ν (denoted by K (cid:96) )give a handle on ∣ V us ∣ independently. The traditional way of determining ∣ V us ∣ relies onthe experimental measurements of K L → πeν to avoid the isospin-breaking eﬀects ( π - η mixing) in the charged kaon decays and the complication from the second (scalar) form Notice, however, that this quoted value does not include the contributions from several new nuclearcorrections investigated in Refs.[4, 5]. K L → π(cid:96)ν , K ± → π (cid:96) ± ν and K S → πeν with (cid:96) = e, µ are used to determine ∣ V us ∣ viathe master formula [1]Γ K(cid:96) = G F m K π S EW ( + δ (cid:96)K + δ SU ) C ∣ V us ∣ f + ( ) I (cid:96)K . (2)Here, Γ K(cid:96) is the K (cid:96) decay width, G F is the Fermi constant, m K is the kaon mass, S EW isthe short-distance radiative correction, δ (cid:96)K is the long-distance radiative correction, δ SU isthe strong isospin-violating eﬀect, C is 1 for the neutral kaon decay and 1 / f + ( q ) is the K → π − vector form factor and I (cid:96)K is the phase-space integral whichcontains the information of the momentum dependence in the form factors. Averaging overthe experimental measurements with appropriate theory inputs of various Standard Modelcorrections, the product ∣ V us ∣ f + ( ) is given as [12] f + ( )∣ V us ∣ = . ( ) , (3)with the uncertainty dominated by the experimental measurements and RCs. The formfactor f + ( ) can be provided by lattice QCD calculations [13–17]. The FLAG average [18]for N f = + + f + ( ) = . ( ) according to an update on December2020, which results in a determination of ∣ V us ∣ = . ( ) exp + RC ( ) lat , for K (cid:96) . (4)High-precision experimental data on K µ and π µ decays [19, 20] also accurately determinethe ratio ∣ V us / V ud ∣ f K ± / f π ± = . ( ) [12]. Employing the FLAG N f = + + f K ± / f π ± = . ( ) , a value of ∣ V us ∣ = . ( ) is obtained, which has a 2 . σ deviation from the K (cid:96) -based value. Combining the ∣ V us ∣ from K (cid:96) and K µ decays yield ∣ V us ∣ = . ( ) , weighted average of K (cid:96) and K µ . (5)It should also be mentioned that ∣ V us ∣ obtained from hyperon and tau decays are given by ∣ V us ∣ = . ( ) [25] and 0.2221(13) [26], respectively, both having larger uncertaintiesthan the kaon decays. Here ∣ V us ∣ is slightly diﬀerent from the PDG value 0 . ( ) due to the update of the FLAG value of f + ( ) . Correspondingly, the value of ∣ V ud ∣ + ∣ V us ∣ + ∣ V ub ∣ given in Eq. (1) also slightly diﬀers from thePDG value of 0 . ( ) V ud ( ) V us .

3o gain a better understanding of the violation of the ﬁrst-row CKM unitarity in Eq. (1)and the disagreement in the determination of ∣ V us ∣ between the K (cid:96) and K µ , for the K (cid:96) decays it requires both a more precise determination of the form factor f + ( ) and a directcalculation of RCs from lattice QCD. The latter is more challenging due to the inclusion ofboth weak and electromagnetic currents in the calculation and is the focus in this paper.Recently, the horizons of lattice QCD studies have been extended to include variousprocesses with higher-order electroweak interactions. The examples include kaon mixing [27–29], rare kaon decays [30–35], double beta decays [36–44], inclusive B -meson decays [45–47],as well as the electromagentic and radiative corrections to the weak decays [48–54]. Amongall these processes, the lattice QCD calculation of RCs in K (cid:96) still remains one of thelargest challenges as it essentially involves a computation of ﬁve-point correlation functions.In Ref. [55], it proposes a new method which bridges the lattice QCD calculation withchiral perturbation theory (ChPT) [56, 57]. For the K (cid:96) decay in the ﬂavor SU(3) limit, itdemonstrates that the lattice QCD calculation of the axial γW -box diagrams can provide allunknown low-energy constants (LECs) that enter the long-distance radiative correction δ (cid:96)K in the ChPT representation at the order of O ( e p ) , thus removing the dependence of theRCs on the model used to estimate these LECs. In this paper we will ﬁrst brieﬂy introducethe methodology and then present the lattice calculation of RCs. II. METHODOLOGY

We start the discussion of the treatment of RCs in K (cid:96) decays with two theoreticalframeworks: Sirlin’s representation and the ChPT representation. Figure 1. The γW -box diagrams for the semileptonic decay process H i → H f e ¯ ν e . Sirlin’s representation is particularly useful in the treatment of the semileptonic decay H i → H f e ¯ ν e with the hadrons H i and H f having nearly the same masses m i ≈ m f . In this4ase, the O ( G F α e ) RCs to the decay width is given as [58] δ = α e π [ ¯ g + m Z m p + ln m Z m W + ˜ a g ] + δ QEDHO + ◻ V AγW , (6)where m Z and m W are the masses for the Z and W bosons. m p is the proton mass thatenters simply by convention. The Sirlin’s function ¯ g , which is a function of the electron’send-point energy, summarizes the infrared-singular contributions involving both the one-loop and bremsstrahlung corrections [58–60]. The O ( α s ) QCD correction ˜ a g is dominatedby the high-energy scale Q ≃ m W with a relatively small contribution of α e π ˜ a g ≈ − . × − [58, 61]. The contribution from the resummation of the large QED logs is containedin δ QEDHO = . ( ) [62]. All the contributions that are sensitive to hadronic scales, residein the axial γW -box contribution ◻ V AγW , as shown in Fig. 1. The total contribution δ isequivalent to ( S EW − ) + δ eK shown in Eq. (2). K ( K ) π + ( K + ) J emµ J W,Aν (A) K ( K ) π + ( K + ) J emµ J W,Aν (B) K ( K ) π + ( K + ) J emµ J W,Aν (C)

Figure 2. Quark contractions for K → π + e ¯ ν e and K → K + e ¯ ν e . In the K (cid:96) decays, since m K is not close to m π , the non-perturbative hadronic eﬀectsare contained not only in ◻ V AγW , but also in other diagrams. As a consequence, Eq. (6) cannot be used directly. To evaluate the total RCs, the calculation of the ﬁve-point correlationfunction is required. To simplify this problem, Ref. [55] proposes to calculate the RCs for K → π + e ¯ ν e in the ﬂavor SU(3) limit, where m K = m π . The relevant contractions are shownin Fig. 2 with the disconnected diagram (C) vanishing in the ﬂavor SU(3) limit. Althoughthe physical value of δ cannot be determined directly using this unphysical setup, the latticecalculation can help to extract the LECs for ChPT. Then by using ChPT one can obtain thephysical RCs. Besides for the K → π + transition, the semileptonic decay of K → K + e ¯ ν e can also be used to determine the same LECs as it has the same contractions as K → π + up to the disconnected parts.In ChPT, the RCs to K (cid:96) are computed to O ( e p ) [56, 57, 63] with the short-distance5adiative correction S EW = − e [− π ln M Z M ρ + ( X phys6 ) α s ] + δ QEDHO = . ( ) , (7)where M ρ = .

77 GeV is the rho mass and ( X phys6 ) α s ≈ . × − [64] summarizes the O ( α s ) pQCD contribution to X phys6 with X phys6 ( µ ) ≡ X r ( µ ) − K r ( µ ) the combination of tworenormalized LECs. The scale µ is usually taken as µ = M ρ in the numerical analysis. Thelong-distance radiative correction δ (cid:96)K has the dependence on the LECs through the relation δ (cid:96)K ± = e [− X −

12 ˜ X phys6 ( M ρ )] + ⋯ δ (cid:96)K = e [ X −

12 ˜ X phys6 ( M ρ )] + ⋯ , (8)where the ellipses indicates the omission of the known kinematic terms, which does not de-pend on the LECs. X and ˜ X phys6 are LECs relevant at O ( e p ) . ˜ X phys6 ( M ρ ) ≡ X phys6 ( M ρ ) +( π ) − ln ( M Z / M ρ ) − ( X phys6 ) α s removes the large electroweak logarithm and the O ( α s ) pQCD correction from X phys6 . In a similar way, one can deﬁne the quantity δ (cid:96)π ± for π (cid:96) δ (cid:96)π ± = e [− X −

12 ˜ X phys6 ( M ρ )] + ⋯ . (9)Since the neutral kaon decay mode K → π + e ¯ ν e is theoretically cleaner as it does notreceive contributions from the π - η mixing which complicates the analysis in the ﬂavor SU(3)limit, we may use it to extract the LECs. Comparing the ChPT and Sirlin’s representations,the relation between the axial γW -box contribution ◻ V AγW ∣ K , SU ( ) and the LECs is givenby [55] − X + ¯ X phys6 ( M ρ ) = − πα (◻ V AγW ∣ K , SU ( ) − α π ln M W M ρ ) + π ( − ˜ a g ) . (10)with ¯ X phys6 ( M ρ ) deﬁned as ¯ X phys6 ( M ρ ) ≡ X phys6 ( M ρ ) + ( π ) − ln ( M Z / M ρ ) , which removesonly the large electroweak logarithm but retains the full pQCD corrections. For the π (cid:96) decay, the relation is given by43 X + ¯ X phys6 ( M ρ ) = − πα (◻ V AγW ∣ π − α π ln M W M ρ ) + π ( − ˜ a g ) . (11) Notice that in a similar expression in Ref. [55], the quantity δ (cid:96)K ± includes also contributions from theLECs { K ri } . That, however, was not the standard convention adopted by the ChPT community, whichchooses to lump the { K ri } contribution into δ SU . ◻ V AγW ∣ π for the π (cid:96) decay has been calculated in Ref. [8]. The focus ofthis paper is on the determination of ◻ V AγW ∣ K , SU ( ) , from which the LECs X and ¯ X phys6 ( M ρ ) can be obtained.The lattice QCD calculation of ◻ V AγW ∣ K , SU ( ) can follow the procedures given in Ref. [8].We ﬁrst deﬁne the hadronic function H V Aµν ( t, ⃗ x ) in Euclidean space H V Aµν ( t, ⃗ x ) ≡ ⟨ π + ( P )∣ T [ J emµ ( t, ⃗ x ) J W,Aν ( )] ∣ K ( P )⟩ , (12)where J emµ = ¯ uγ µ u − ¯ dγ µ d − ¯ sγ µ s is the electromagnetic quark current, and J W,Aν = ¯ uγ ν γ s is the axial part of the weak charged current. The Euclidean momentum P is chosen as P = ( im K , ⃗ ) with m K = m π in the ﬂavor SU(3) limit. The box contribution ◻ V AγW ∣ K , SU ( ) can be determined through the integral ◻ V AγW ∣ K , SU ( ) = α e π ∫ dQ Q m W m W + Q M K ( Q ) (13)with M K ( Q ) = − √ Q m K ∫ d x ω ( t, ⃗ x ) (cid:15) µνα x α H V Aµν ( t, ⃗ x ) ,ω ( t, ⃗ x ) = ∫ π − π cos θ dθπ j (√ Q ∣⃗ x ∣ cos θ )∣⃗ x ∣ cos (√ Q t sin θ ) . (14)Here j ( x ) is the spherical Bessel function. To compute the integral in Eq. (13), for small Q , we use lattice QCD input of H V Aµν ( t, ⃗ x ) . For large Q , the operator product expansion of J emµ ( x ) J W,Aν ( ) is utilized with the Wilson coeﬃcients given at the four-loop accuracy [65,66]. For more details, we refer the readers to Ref. [8]. III. NUMERICAL RESULTS

Five gauge ensemble with N f = + O ( ) random spacetimelocations. The correlation functions are constructed using the ﬁeld sparsening technique [67,68] with a signiﬁcant reduction in the propagator storage. For the locations of two current7 nsemble m π [MeV] L T a − [GeV]24D 141.2(4) 24 64 1 . . . . . m π , the spatial andtemporal extents, L and T , the inverse of lattice spacing a − . insertions J emµ and J W,Aν , we treat one as the source of the propagator and the other as thesink. In this way the hadronic function H V Aµν ( t, ⃗ x ) , which depends on the coordinate-spacevariable x can be obtained. Such technique has also been used in the computation of three-point correction function to extract the pion charge radius [69]. The ﬂavor SU(3) limit isachieved by tuning down the strange quark mass to be the same as the light quark mass.Inserting H V Aµν ( t, ⃗ x ) into the integral (14), we calculate the scalar function M K ( Q ) . Thelattice results for M K ( Q ) as a function of Q are shown in the left panel of Fig. 3. Atlarge Q ( Q ≳ ), the lattice results from diﬀerent gauge ensembles start to disagree,suggesting the obvious lattice discretization eﬀects. In the right panel of Fig. 3, a continuumextrapolation is performed to obtain the results in the continuum limit for Iwasaki and DSDRensembles separately. To reduce the systematic uncertainties contained in the lattice dataat large Q , we calculate the M K ( Q ) in pQCD using the RunDec package [70]. At low Q the perturbative results suﬀer from large pQCD truncation eﬀects due to the lack ofhigher-loop and higher-twist contributions. We observe an expected discrepancy betweenthe orange and magenta curves at low Q , where the former uses the 4-ﬂavor theory downto 1 GeV, while the latter turns to the 3-ﬂavor theory upon decoupling the charm quark at1.6 GeV.We introduce a momentum-squared scale Q that separates the Q -integral into tworegimes. We use the lattice data to determine the integral for Q ≤ Q and perturbationtheory to determine the integral for Q > Q . Three values of Q = , , are used tocheck the Q -dependence in the ﬁnal results. The corresponding results for ◻ V A, ≤ γW ∣ K , SU ( ) and ◻ V A, > γW ∣ K , SU ( ) are listed in Table II. 8 [GeV ]00.020.040.060.08 M K ( Q ) [GeV ]00.020.040.060.08 Cont. Limit, IwasakiCont. Limit, DSDRPT (n f =4 match with n f =3)PT (n f =4 down to 1 GeV) Figure 3. M K ( Q ) as a function of Q . In the left panel, the lattice results for all ﬁve ensembles aregiven. In the right panel, we have extrapolated the Iwasaki and DSDR results to their continuumlimit. The remaining orange and magenta curves are the results from perturbation theory. After combining the lattice data and perturbative results given in Table II, we have ◻ V AγW ∣ K , SU ( ) = ⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩ . ( ) stat ( ) PT ( ) a ( ) FV × − Q = . ( ) stat ( ) PT ( ) a ( ) FV × − Q = . . ( ) stat ( ) PT ( ) a ( ) FV × − Q = (15)Here we take the combination of the Iwasaki and perturbative results as the central valueand estimate the residual lattice artifacts (with a subscript a ) using the discrepancy betweenIwasaki and DSDR. The lattice ﬁnite-volume eﬀects (with a subscript FV) are estimated bycomparing the 24D and 32D results. As a ﬁnal result, we quote the value of ◻ V AγW ∣ K , SU ( ) at Q = and add the statistical and systematic errors in quadrature ◻ V AγW ∣ K , SU ( ) = . ( ) × − . (16)9 V A, ≤ γW ∣ K , SU ( ) ◻ V A, > γW ∣ K , SU ( ) Q Iwasaki DSDR pQCD1 GeV . ( ) × − . ( ) × − . ( ) × − . ( ) × − . ( ) × − . ( ) × − . ( ) × − . ( ) × − . ( ) × − Table II. Using the scale Q to split the integral range, the contributions of ◻ V A, ≤ γW ∣ K , SU ( ) fromlattice QCD and ◻ V A, > γW ∣ K , SU ( ) from perturbation theory are shown. For the lattice results, wehave performed the continuum extrapolation for Iwasaki and DSDR ensembles separately. For theperturbative results, the central values are compiled using the 4-ﬂavor theory and uncertaintiesinclude the higher-loop eﬀects and the higher-twist eﬀects with the error analysis following Ref. [8]. Inserting the result of ◻ V AγW ∣ K , SU ( ) into Eq. (10), we obtain − X + ¯ X phys6 = . ( . ) × − or − X + ˜ X phys6 = . ( . ) × − . (17)The previous ChPT analysis [57] quoted the LECs from the minimal resonance model [64, 71]with X = − . ( . ) × − , ˜ X phys6 = . ( . ) × − . (18)As it is hard to accurately estimate the uncertainty in these LECs from the ChPT perspec-tive, Ref. [57] attributed to them a 100% uncertainty. Combining X and ˜ X phys6 togetheryields − X + ˜ X phys6 = . ( . ) × − [minimal resonance model] . (19)Our result for − X + ˜ X phys6 agrees with the minimal resonance model within few percent.Such a good agreement could easily be fortuitous as the methods used in the two studiesare very diﬀerent and a large uncertainty is assigned to the estimate based on the model.For the π (cid:96) decay, substituting the lattice QCD result ◻ V AγW ∣ π = . ( ) × − [8] intoEq. (11) yields 43 X + ¯ X phys6 = . ( ) × − . (20)Combining Eqs. (17) and (20) together, we have X = − . ( ) × − , ¯ X phys6 = . ( ) × − . (21)10ere the uncertainty is estimated conservatively through a linear addition. It should bepointed out that in Eq. (21) the estimate of the higher order terms in the ChPT expansionare not included yet.In ChPT, the RCs δ (cid:96)K have two major sources of theoretical uncertainties: the input ofthe LECs at O ( e p ) and the unknown O ( e p ) terms in the ChPT expansion. Using theLECs from this calculation, the former uncertainty is signiﬁcantly reduced, while the latterone remains. It results in an update of δ (cid:96)K (in units of %) δ eK = . ( ) e p ( ) LEC → . ( ) δ µK = . ( ) e p ( ) LEC → . ( ) δ eK ± = . ( ) e p ( ) LEC → − . ( ) δ µK ± = . ( ) e p ( ) LEC → − . ( ) . (22)We refrain from presenting a corresponding update of ∣ V us ∣ in this paper, because (1) ourresults for δ (cid:96)K still agree with the existing literature within error bars, and (2) our latticecalculation removes only the LEC uncertainty but not the dominant, O( e p ) uncertainty.Therefore, we shall instead await a next round of global analysis in the near future such asthat in Ref. [11], whose impact on the precision low-energy tests will be more signiﬁcant.Our lattice result may serve as an important input to such an analysis. IV. CONCLUSION

Modern-day lattice QCD has reached the era when realistic calculations for many inter-esting second-order electroweak processes have become feasible. In this work we perform astudy of the γW -box correction to the kaon semileptonic decay K → π + e ¯ ν e . We adopt thenew method proposed in Ref. [55], which connects the Sirlin’s representation to the ChPTrepresentation in the ﬂavor SU(3) limit. It allows us to determine the LECs for ChPT bycomputing the axial γW -box correction. We ﬁnd that the values of the LECs devised fromthe lattice calculation agree well with the minimal resonance model used in the literature,while a dramatic reduction of the respective uncertainties is achieved. Finally, these LECsare used to estimate the RCs δ (cid:96)K and help to reduce its uncertainty. To further improve thedetermination of RCs, the inclusion of higher-order terms in ChPT and the lattice QCDcomputation of the complete set of Feynman diagrams are necessary.11 CKNOWLEDGMENTS

X.F. and L.C.J. gratefully acknowledge many helpful discussions with our colleaguesfrom the RBC-UKQCD Collaborations. X.F. and P.X.M. were supported in part by NSFCof China under Grants No. 11775002 and No. 12070131001 and National Key Research andDevelopment Program of China under Contracts No. 2020YFA0406400. M.G. is supportedby EU Horizon 2020 research and innovation programme, STRONG-2020 project, undergrant agreement No 824093 and by the German-Mexican research collaboration Grant No.278017 (CONACyT) and No. SP 778/4-1 (DFG). L.C.J. acknowledges support by DOEOﬃce of Science Early Career Award DE-SC0021147 and DOE grant DE-SC0010339. Thework of C.Y.S. is supported in part by the DFG (Project-ID 196253076-TRR 110) andthe NSFC (Grant No. 11621131001) through the funds provided to the Sino-German CRC110 “Symmetries and the Emergence of Structure in QCD”, and also by the Alexandervon Humboldt Foundation through the Humboldt Research Fellowship. The computationis performed under the ALCC Program of the US DOE on the Blue Gene/Q (BG/Q)Mira computer at the Argonne Leadership Class Facility, a DOE Oﬃce of Science Facilitysupported under Contract DE-AC02-06CH11357. Computations for this work were carriedout in part on facilities of the USQCD Collaboration, which are funded by the Oﬃce ofScience of the U.S. Department of Energy. The calculation is also carried out on Tianhe 3prototype at Chinese National Supercomputer Center in Tianjin. [1] P. A. Zyla et al. (Particle Data Group), PTEP , 083C01 (2020).[2] J. C. Hardy and I. S. Towner, Phys. Rev.

C91 , 025501 (2015), arXiv:1411.5987 [nucl-ex].[3] J. C. Hardy and I. S. Towner, Phys. Rev. C , 045501 (2020).[4] C. Y. Seng, M. Gorchtein, and M. J. Ramsey-Musolf, Phys. Rev.

D100 , 013001 (2019),arXiv:1812.03352 [nucl-th].[5] M. Gorchtein, Phys. Rev. Lett. , 042503 (2019), arXiv:1812.04229 [nucl-th].[6] C.-Y. Seng, M. Gorchtein, H. H. Patel, and M. J. Ramsey-Musolf, Phys. Rev. Lett. ,241804 (2018), arXiv:1807.10197 [hep-ph].[7] M. Tanabashi et al. (Particle Data Group), Phys. Rev.

D98 , 030001 (2018).

8] X. Feng, M. Gorchtein, L.-C. Jin, P.-X. Ma, and C.-Y. Seng, Phys. Rev. Lett. , 192002(2020), arXiv:2003.09798 [hep-lat].[9] C.-Y. Seng, X. Feng, M. Gorchtein, and L.-C. Jin, Phys. Rev. D , 111301 (2020),arXiv:2003.11264 [hep-ph].[10] C.-Y. Seng and U.-G. Meißner, Phys. Rev. Lett. , 211802 (2019), arXiv:1903.07969 [hep-ph].[11] M. Antonelli et al. (FlaviaNet Working Group on Kaon Decays), Eur. Phys. J. C , 399(2010), arXiv:1005.2323 [hep-ph].[12] M. Moulson, PoS CKM2016 , 033 (2017), arXiv:1704.04104 [hep-ex].[13] P. A. Boyle et al. (RBC/UKQCD), JHEP , 164 (2015), arXiv:1504.01692 [hep-lat].[14] N. Carrasco, P. Lami, V. Lubicz, L. Riggio, S. Simula, and C. Tarantino, Phys. Rev. D ,114512 (2016), arXiv:1602.04113 [hep-lat].[15] S. Aoki, G. Cossu, X. Feng, H. Fukaya, S. Hashimoto, T. Kaneko, J. Noaki, and T. Onogi(JLQCD), Phys. Rev. D , 034501 (2017), arXiv:1705.00884 [hep-lat].[16] A. Bazavov et al. (Fermilab Lattice, MILC), Phys. Rev. D , 114509 (2019), arXiv:1809.02827[hep-lat].[17] J. Kakazu, K.-i. Ishikawa, N. Ishizuka, Y. Kuramashi, Y. Nakamura, Y. Namekawa,Y. Taniguchi, N. Ukita, T. Yamazaki, and T. Yoshi´e (PACS), Phys. Rev. D , 094504(2020), arXiv:1912.13127 [hep-lat].[18] S. Aoki et al. (Flavour Lattice Averaging Group), Eur. Phys. J. C80 , 113 (2020),arXiv:1902.08191 [hep-lat].[19] W. J. Marciano, Phys. Rev. Lett. , 231803 (2004), arXiv:hep-ph/0402299.[20] F. Ambrosino et al. (KLOE), Phys. Lett. B , 76 (2006), arXiv:hep-ex/0509045.[21] R. J. Dowdall, C. T. H. Davies, G. P. Lepage, and C. McNeile, Phys. Rev. D , 074504(2013), arXiv:1303.1670 [hep-lat].[22] N. Carrasco et al. , Phys. Rev. D , 054507 (2015), arXiv:1411.7908 [hep-lat].[23] A. Bazavov et al. , Phys. Rev. D , 074512 (2018), arXiv:1712.09262 [hep-lat].[24] N. Miller et al. , Phys. Rev. D , 034507 (2020), arXiv:2005.04795 [hep-lat].[25] N. Cabibbo, E. C. Swallow, and R. Winston, Phys. Rev. Lett. , 251803 (2004), arXiv:hep-ph/0307214.[26] Y. S. Amhis et al. (HFLAV), (2019), arXiv:1909.12524 [hep-ex].

27] N. H. Christ, T. Izubuchi, C. T. Sachrajda, A. Soni, and J. Yu (RBC, UKQCD), Phys. Rev.

D88 , 014508 (2013), arXiv:1212.5931 [hep-lat].[28] Z. Bai, N. H. Christ, T. Izubuchi, C. T. Sachrajda, A. Soni, and J. Yu, Phys. Rev. Lett. ,112003 (2014), arXiv:1406.0916 [hep-lat].[29] N. H. Christ, X. Feng, G. Martinelli, and C. T. Sachrajda, Phys. Rev. D , 114510 (2015),arXiv:1504.01170 [hep-lat].[30] N. H. Christ, X. Feng, A. Portelli, and C. T. Sachrajda (RBC, UKQCD), Phys. Rev. D ,094512 (2015), arXiv:1507.03094 [hep-lat].[31] N. H. Christ, X. Feng, A. Portelli, and C. T. Sachrajda (RBC, UKQCD), Phys. Rev. D ,114517 (2016), arXiv:1605.04442 [hep-lat].[32] N. H. Christ, X. Feng, A. Juttner, A. Lawson, A. Portelli, and C. T. Sachrajda, Phys. Rev.D , 114516 (2016), arXiv:1608.07585 [hep-lat].[33] Z. Bai, N. H. Christ, X. Feng, A. Lawson, A. Portelli, and C. T. Sachrajda, Phys. Rev. D ,074509 (2018), arXiv:1806.11520 [hep-lat].[34] N. H. Christ, X. Feng, A. Portelli, and C. T. Sachrajda (RBC, UKQCD), Phys. Rev. D ,114506 (2019), arXiv:1910.10644 [hep-lat].[35] N. H. Christ, X. Feng, L.-C. Jin, and C. T. Sachrajda, Phys. Rev. D , 014507 (2021),arXiv:2009.08287 [hep-lat].[36] 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].[37] P. E. Shanahan, B. C. Tiburzi, M. L. Wagman, F. Winter, E. Chang, Z. Davoudi, W. Detmold,K. Orginos, and M. J. Savage, Phys. Rev. Lett. , 062003 (2017), arXiv:1701.03456 [hep-lat].[38] A. Nicholson et al. , Phys. Rev. Lett. , 172501 (2018), arXiv:1805.02634 [nucl-th].[39] X. Feng, L.-C. Jin, X.-Y. Tuo, and S.-C. Xia, Phys. Rev. Lett. , 022001 (2019),arXiv:1809.10511 [hep-lat].[40] X.-Y. Tuo, X. Feng, and L.-C. Jin, Phys. Rev.

D100 , 094511 (2019), arXiv:1909.13525 [hep-lat].[41] W. Detmold and D. J. Murphy (NPLQCD), (2020), arXiv:2004.07404 [hep-lat].[42] X. Feng, L.-C. Jin, Z.-Y. Wang, and Z. Zhang, (2020), arXiv:2005.01956 [hep-lat].[43] Z. Davoudi and S. V. Kadam, Phys. Rev. D , 114521 (2020), arXiv:2007.15542 [hep-lat].

44] Z. Davoudi and S. V. Kadam, (2020), arXiv:2012.02083 [hep-lat].[45] M. T. Hansen, H. B. Meyer, and D. Robaina, Phys. Rev. D , 094513 (2017),arXiv:1704.08993 [hep-lat].[46] S. Hashimoto, PTEP , 053B03 (2017), arXiv:1703.01881 [hep-lat].[47] P. Gambino and S. Hashimoto, Phys. Rev. Lett. , 032001 (2020), arXiv:2005.13730 [hep-lat].[48] N. Carrasco, V. Lubicz, G. Martinelli, C. T. Sachrajda, N. Tantalo, C. Tarantino, andM. Testa, Phys. Rev. D91 , 074506 (2015), arXiv:1502.00257 [hep-lat].[49] V. Lubicz, G. Martinelli, C. T. Sachrajda, F. Sanﬁlippo, S. Simula, and N. Tantalo, Phys.Rev.

D95 , 034504 (2017), arXiv:1611.08497 [hep-lat].[50] D. Giusti, V. Lubicz, G. Martinelli, C. T. Sachrajda, F. Sanﬁlippo, S. Simula, N. Tantalo,and C. Tarantino, Phys. Rev. Lett. , 072001 (2018), arXiv:1711.06537 [hep-lat].[51] X. Feng and L. Jin, Phys. Rev. D , 094509 (2019), arXiv:1812.09817 [hep-lat].[52] N. H. Christ, X. Feng, J. Lu-Chang, and C. T. Sachrajda, PoS

LATTICE2019 , 259 (2020).[53] A. Desiderio et al. , Phys. Rev. D , 014502 (2021), arXiv:2006.05358 [hep-lat].[54] R. Frezzotti, M. Garofalo, V. Lubicz, G. Martinelli, C. T. Sachrajda, F. Sanﬁlippo, S. Simula,and N. Tantalo, (2020), arXiv:2012.02120 [hep-ph].[55] C.-Y. Seng, X. Feng, M. Gorchtein, L.-C. Jin, and U.-G. Meißner, JHEP , 179 (2020),arXiv:2009.00459 [hep-lat].[56] V. Cirigliano, M. Knecht, H. Neufeld, H. Rupertsberger, and P. Talavera, Eur. Phys. J. C , 121 (2002), arXiv:hep-ph/0110153.[57] V. Cirigliano, M. Giannotti, and H. Neufeld, JHEP , 006 (2008), arXiv:0807.4507 [hep-ph].[58] A. Sirlin, Rev. Mod. Phys. , 573 (1978), [Erratum: Rev. Mod. Phys.50,905(1978)].[59] A. Sirlin, Phys. Rev. , 1767 (1967).[60] D. H. Wilkinson and B. E. F. Maceﬁeld, Nucl. Phys. A158 , 110 (1970).[61] C.-Y. Seng, D. Galviz, and U.-G. Meißner, JHEP , 069 (2020), arXiv:1910.13208 [hep-ph].[62] J. Erler, Rev. Mex. Fis. , 200 (2004), arXiv:hep-ph/0211345 [hep-ph].[63] V. Cirigliano, M. Knecht, H. Neufeld, and H. Pichl, Eur. Phys. J. C , 255 (2003), arXiv:hep-ph/0209226.[64] S. Descotes-Genon and B. Moussallam, Eur. Phys. J. C , 403 (2005), arXiv:hep-ph/0505077.[65] S. A. Larin and J. A. M. Vermaseren, Phys. Lett. B259 , 345 (1991).

66] P. A. Baikov, K. G. Chetyrkin, and J. H. Kuhn, Phys. Rev. Lett. , 132004 (2010),arXiv:1001.3606 [hep-ph].[67] Y. Li, S.-C. Xia, X. Feng, L.-C. Jin, and C. Liu, Phys. Rev. D , 014514 (2021),arXiv:2009.01029 [hep-lat].[68] W. Detmold, D. J. Murphy, A. V. Pochinsky, M. J. Savage, P. E. Shanahan, and M. L.Wagman, (2019), arXiv:1908.07050 [hep-lat].[69] X. Feng, Y. Fu, and L.-C. Jin, Phys. Rev. D , 051502 (2020), arXiv:1911.04064 [hep-lat].[70] K. G. Chetyrkin, J. H. Kuhn, and M. Steinhauser, Comput. Phys. Commun. , 43 (2000),arXiv:hep-ph/0004189 [hep-ph].[71] B. Ananthanarayan and B. Moussallam, JHEP , 047 (2004), arXiv:hep-ph/0405206., 047 (2004), arXiv:hep-ph/0405206.