Improved $K_{e3}$ radiative corrections sharpen the $K_{\mu 2}$--$K_{l3}$ discrepancy

The measurements of $V_{us}$ in leptonic $(K_{\mu 2})$ and semileptonic $(K_{l3})$ kaon decays exhibit a $3\sigma$ disagreement, which could originate either from physics beyond the Standard Model or some large unidentified Standard Model systematic effects. Clarifying this issue requires a careful examination of all existing Standard Model inputs. Making use of a newly-proposed computational framework and the most recent lattice QCD results, we perform a comprehensive re-analysis of the electroweak radiative corrections to the $K_{e3}$ decay rates that achieves an unprecedented level of precision of $10^{-4}$, which improves the current best results by almost an order of magnitude. No large systematic effects are found, which suggests that the electroweak radiative corrections should be removed from the ``list of culprits'' responsible for the $K_{\mu 2}$--$K_{l3}$ discrepancy.


I. INTRODUCTION
Despite the discovery of the Higgs boson in the year 2012 [1,2]  The unitarity of the CKM matrix is a rigorous SM prediction [3,4]. In particular, the top-row CKM unitarity (which is also known as the Cabibbo unitarity) that involves the matrix elements V ud and V us (V ub is negligible) has received the most attention because they can be measured to high precision in hadron and nuclear beta decays. Recently, a series of improvements in the theory [5][6][7][8] of the electroweak radiative corrections (RC) in the extraction of V ud led to an apparent deviation of the Cabibbo unitarity at a level of 3σ [9].
However, in this work we will not focus on V ud , but rather on V us which possesses yet another interesting anomaly by itself.
Let us focus on the two best determinations of the matrix element V us , which come from leptonic (K l2 ) and semileptonic (K l3 ) kaon decays respectively. From the leptonic kaon and pion decay, the following ratio is obtained: where f K + and f π + are the K + and the π + decay constant, respectively, which require lattice QCD inputs. The theory uncertainty on the right-hand side is less than 10 −3 , thanks to the cancellation of the common electroweak RC to the leptonic kaon and pion decay rate [10,11].
Combining this expression with the N f = 2 + 1 + 1 FLAG average of f K + /f π + [12] and the recent value of V ud obtained from superallowed beta decays [5], the following result is quoted in PDG 2020 [9]: |V us | = 0.2252(5) (K µ2 /π µ2 + superallowed) (2) Meanwhile, in the semileptonic kaon decay process K → πl + ν(γ) one does not measure a ratio, but obtains V us directly from the decay rate, where the SM inputs include the electroweak RC, the Kπ form factors and the SU(2) isospin-breaking effects (we postpone the detailed discussions to the main text). With the most recent theory inputs of these quantities, PDG 2020 quotes the following result: We observe a ∼ 3σ disagreement between the numbers in Eq. (2) and (3), with a ∼ 1% difference between the two central values. This provides another interesting hint to the existence of BSM physics [13][14][15][16][17][18][19][20][21][22] which, to some extent, is even more promising than the top-row CKM unitarity deficit. In fact, the extraction of V us is free from complicated nuclear-structure uncertainties (except those that enter V ud in Eq. (1), whose effect on V us is subdominant to the existing uncertainties). For instance, if the total uncertainty in Eqs. (2) and (3) is reduced to 4 × 10 −4 or below, with the central values unchanged, the discrepancy will reach 5σ which is sufficient to claim an observation of a BSM signal. Achieving this final goal requires a careful re-analysis of all the SM inputs, not just to reduce their uncertainties but also to make sure that no large unidentified SM corrections were missed in existing analyses.
In this work, we study a particularly important SM correction to the kaon semileptonic decay, namely the electroweak RC. Earlier studies of this topic by Ginsberg [23][24][25][26], Becherrawy [27] and later by Bytev et al. [28] and Andre [29] assumed specific models for the strong and electroweak interactions which made a rigorous analysis of the theory uncertainties rather challenging. Another class of works, e.g. by García and Maya [30] and by Juárez-León et al. [31][32][33] put more emphasis on the so-called "model-independent" piece in the long-distance electromagnetic corrections (i.e. the convection term contribution, which we will explain in the main text) but were unable to place any constrain on the "modeldependent" piece originating from non-perturbative Quantum Chromodynamics (QCD) at the hadronic scale. So far, the only approach that allows a systematic error analysis in every part of the electroweak RC has been the chiral perturbation theory (ChPT) calculation by Cirigliano et al. [34][35][36], where the most general electroweak interactions between hadrons and dynamical photons [37] and leptons [38] are arranged according to increasing powers of p/Λ χ , where p is a typical small momentum scale in such interactions and Λ χ 4πF π is the chiral symmetry breaking scale, with F π = 92.1 MeV the pion decay constant. Within this framework, the long-distance electromagnetic RC to K l3 decay is calculated to O(e 2 p 2 ), and the theory uncertainty comes from two major sources: The unknown low-energy constants (LECs) at O(e 2 p 2 ), and the neglected contributions of the order O(e 2 p 4 ). Both uncertainties are estimated to be of the order 10 −3 . At this point it seems formidable to make any further progress within the same theory framework, because (1) the LECs are only calculable within phenomenological models [39,40] with outcomes that are highly uncertain, and (2) to reduce the higher-order corrections one needs to perform a full two-loop ChPT calculation which is not only technically challenging but more importantly, involves even more unknown LECs.
A series of preparatory works were done since early 2020 in order to eventually overcome the difficulties mentioned above. First, a new theory framework based on the hybridization of the classical Sirlin's approach [41,42] and modern ChPT was formulated [43] in order to resum the most important O(e 2 p 2n ) effects while retaining the full model-independent characteristics in the traditional ChPT approach. Next, lattice QCD was introduced to study the part of the RC in semileptonic decays that carries the largest hadronic uncertainties, namely the axial γW -box diagram. The first calculation was done on the pion [44], which removed the dominant theory uncertainty in the semileptonic pion decay and also confirmed the result of the previous dispersion-relation analysis of the RC in free neutron [45]. Shortly after that, following the suggestion in Ref. [46] a new lattice calculation of the Kπ axial γWbox in the flavor SU(3) limit was performed [47]. Up to this point, we finally have all the necessary ingredients and are in the position to present a fully-updated numerical analysis of the electroweak RC in kaon semileptonic decays that eventually reduces the existing theory uncertainty by almost an order of magnitude, i.e. to the level of 10 −4 .
The main results in this study were presented in an earlier paper [48], and here we will show all the details. We concentrate on the K e3 channel and not K µ3 throughout this study for reasons that will become clear in the main text. The contents of this work are arranged as follows. In Section II we introduce the basic notation and set up our theory framework.
In Sections III-VI we present our update of the contributions from the "virtual" electroweak RC; in particular, we demonstrate in Section VI how the most recent lattice QCD results are used to constrain the hadronic uncertainties in the physical Kπ axial γW -box diagram.
The contribution from the real-photon emission process is calculated in Section VII. In Section VIII we discuss how our new results should be interpreted in the ChPT language, and show the numerical improvement against the existing calculations. Final discussions and conclusions are provided in Section IX.

II. NOTATION AND SETUP
One of the most important avenues to extract V us is the inclusive kaon semileptonic decay K l3 , i.e. the process K(p) → π(p ) + l + (p l ) + ν e (p ν ) + nγ, where l = e, µ, and n ≥ 0 is the number of photons in the final state. It will be evident later that the case l = e allows for a much better control of the theory uncertainties, so throughout this paper, we will concentrate on this particular case. If all massless final-state particles are left unobserved, the differential decay rate of the process is fully described by three independent, dimensionless Lorentz-invariant variables 1 : where P ≡ p − p − p e . Notice that x is strictly zero (neglecting neutrino mass) for n = 0, but may take a non-zero value when n ≥ 1. We may have as well introduced the usual At O(G 2 F ) (where G F = 1.1663787(6) × 10 −5 GeV −2 is the Fermi constant extracted from muon decay [49]), only the n = 0 process contributes to the K e3 decay rate. Its corresponding tree-level amplitude is given by: where the effects of the strong interaction are fully contained in the following hadronic matrix element of the charged weak current: The equation above defines the charged weak form factors f Kπ ± (t) 2 . It is also customary to define a third form factor: 1 In the existing literature x is more often defined as P 2 , which carries a dimension. 2 We wish to remind the readers that our sign convention for the form factors is f Kπ + (0) < 0, which is also adopted in our previous works, e.g. [43,46], but may be opposite to other existing literature. This serves to be consistent with the sign convention of the charged weak current (J µ W ) † derived from ChPT. and call f Kπ + (t) and f Kπ 0 (t) the "vector" and "scalar" form factor, respectively. From the definition above, it is obvious that f Kπ 0 (0) = f Kπ + (0), so another common step is to factor out their t = 0 value:f There are several different ways to parameterizef +,0 (t), e.g. Taylor expansion, monopole parameterization and dispersive parameterization. The interested reader may consult Ref. [50] and references therein for the details, and we will also come back to this point in Section IV.
It is instructive to display explicitly the absolute square of the tree-level amplitude above (upon summing over the lepton spin, as we will always do throughout this work): Here we purposely retain the x-dependence in the formula above despite the fact that x = 0 when n = 0. The x-dependence becomes important later when we discuss the squared amplitude of the bremsstrahlung process. The impact of the form factors f Kπ ± on the treelevel decay rate relies heavily on the leptonic trace in Eq.(9). Suppose we define: then a straightforward calculation shows: where r π ≡ M 2 π /M 2 K and r e ≡ m 2 e /M 2 K . We observe that only H(+1, +1) is not explicitly suppressed by the factor r e ≈ 10 −6 . Following the notations in Appendix A, the decay rate at O(G 2 F ) in given by: From the argument above, it is apparent that only f Kπ + (t), and not f Kπ − (t), is relevant in (Γ K e3 ) tree . Of course the actual value of (Γ K e3 ) tree depends on the specific parameterization off + (t) and the parameters therein, but the impact of the different choices is generically of the order 0.1%. Since in this paper (Γ K e3 ) tree serves only as a normalization factor to the already-small RC, such a difference is completely negligible.
The electroweak RC induces a shift of the tree-level decay rate: (Γ K e3 ) tree → (Γ K e3 ) tree + δΓ K e3 . We define the quantity: that represents the fractional correction to the decay rate, and we will discuss its relation to the different quantities within the ChPT framework in Section VIII. To match the precision level of current and near-future experiments, we need a theoretical prediction of δ K e3 up to O(α). At this level, the only two contributors are (1) the O(G F α) electroweak RC to the n = 0 decay amplitude, and (2) the tree-level contribution from the n = 1 process. We will spend the next few sections discussing these two contributions.

III. VIRTUAL CORRECTION: ANALYTIC PIECES
We start by discussing the virtual corrections, i.e. the O(G F α) electroweak RC to the n = 0 decay amplitude. It is possible to express such corrections entirely in terms of perturbations to the charged weak form factors, i.e 3 , The only complication is that δf Kπ ± are complex functions of two variables, e.g. {y, z}, rather than real functions of a single variable t.
The virtual contribution to δ K e3 at O(α) arises from the interference between M 0 and Again, by restricting ourselves to K e3 , we only need to know δf Kπ + in order to determine the perturbation to the n = 0 squared amplitude: Based on the theory framework outlined in Refs. [43,46], the O(G F α) virtual corrections to 3 Using the on-shell condition, one can show that other leptonic bilinear structures, such as the n = 0 decay amplitude can be summarized by the following equation: Let us briefly explain the notation above, all the details are given in Ref. [43]. First, the terms in the square bracket come from the "weak" RC including its O(α s ) perturbative QCD (pQCD) corrections a pQCD ≈ 0.068, the electron wavefunction renormalization, and the resummation of the large QED logs represented by δ QED HO = 0.0010(3) [51]. An infinitesimal photon mass M γ is introduced to regularize the infrared (IR) divergence in the electron wavefunction renormalization. Next, the quantities δM 2,3 represent the contributions from two separate pieces of the electromagnetic RC to the charged weak form factors, known as the "two-point function" and "three-point function", respectively. Finally, δM γW represents the contribution from the γW -box diagram: where we have introduced the so-called "generalized Compton tensor" T Kπ µν which plays a central role in the upcoming analysis: where T {. . .} denotes the conventional time-ordering. It satisfies the following Ward identities: with q ≡ p + q − p, and The first line in Eq. (19) is a consequence of the exact conservation of the electromagnetic current, while the second line entails the partial conservation of the charged weak current.
Expressing hadronic matrix elements in terms of integrals with respect to T µν is a classical technique in hadron physics that appears also in, e.g., the Cottingham's approach to the hadronic mass splittings [52][53][54][55][56][57].
Using now the following Dirac matrix identity: (with 0123 = −1) one splits the γW -box diagram into two pieces: δM γW = δM a γW + δM b γW , where the antisymmetric tensor is contained in the second term 4 . A great simplification is observed upon combining δM 2 with δM a γW [46]: The terms in the square bracket in Eq. (22) are exactly known as they are isolated from the full one-loop integral with the help of the Ward identities in Eq. (19), as well as the operator product expansion (OPE) at leading-twist in the q ∼ M W region (see Eq.(6.1) in Ref. [43]), andã res g ≈ 0.019 entails the O(α s ) pQCD corrections of such terms. The remaining "integral" piece requires further theoretical analysis and will be treated in the next section.
Meanwhile, the other component of the γW -box diagram reads: which can be split into two pieces, as well: δM b γW = δM b,V γW + δM b,A γW , where δM b,V γW (δM b,A γW ) picks up the contribution from the vector (axial) charged weak current in the generalized Compton tensor T Kπ µν . At this point, we can combine the terms in the square brackets from Eqs. (16) and (22).
They are analytically known and do not require any further treatment. Their contribution to δf Kπ + is given by 4 We used to label them as δM V γW and δM A γW in Ref. [46], but this may cause confusions with notations of box diagrams in some literature when we further divide the contributions from the vector and axial charged weak current in T Kπ µν , so here we adopt an alternative labeling. whereã g = −(3/2)a pQCD +ã res g ≈ −0.083. We use the subscript "I" to signify the fact that it carries an IR divergence. We will see later that the remaining IR-divergent pieces in the virtual corrections come from (δM 2 + δM γW ) int and δM 3 , and will carry the subscript "II" and "III", respectively.
All the remaining O(G F α) electroweak RC to the n = 0 decay amplitude not included in Eq. (24) are fully contained in the following quantities: δM 2 + δM a γW int , δM b,V γW , δM 3 and δM b,A γW . They will be studied in the next three sections. In this section we evaluate the loop integrals in δM 2 + δM a γW int and δM b,V γW . The first important observation is that these integrals cannot depend on physics at large virtual momentum q (so we could take M 2 W /(M 2 W −q 2 ) → 1 in the integrand). In δM 2 + δM a γW int , this is because the numerators in the integrand contain explicit factors of p e , p − p or quark masses (in Γ λ Kπ ); whereas in δM b,V γW , it is because there is no extra antisymmetric tensor coming from T Kπ µν V , so the integral vanishes when q (p − p ) or p e due to symmetry.
Therefore, these integrals are saturated by contributions from the intermediate hadronic states at low energy.
All the information on the hadronic structure in these integrals is contained in the generalized Compton tensor T µν Kπ and the vector Γ µ Kπ . Within the former, we distinguish two types of contributions shown in Fig.1: the pole term associated with a charged meson (initial or final, depending on the reaction channel) propagator which leads to a 1/q behavior in the soft photon limit, and the seagull term which is regular in that limit. The pole term is model-independent and given in terms of the meson weak and electromagnetic form factors, whereas the seagull term, alongside the form factors, contains information about excited states, and is generally model-dependent. It is common to single out the Born part of the generalized Compton tensor, defined as the pole terms complemented by a part of the seagull term that ensures that the Ward identities in Eq. (19) are satisfied. In this way, the remaining, non-Born part is regular for q → 0 and also obeys Ward identities individually.
Guided by the order O(p 2 ) result in chiral expansion for the Compton tensor, and we thus define the minimal Born contributions for the two decay channels as where F π − em (q 2 ) and F K + em (q 2 ) are the electromagnetic form factors of the π − and the K + , respectively 5 , which satisfy F π − em (0) = −1 and F K + em (0) = 1. Furthermore, the normalization of the seagull term is fixed as: One can check that the electromagnetic Ward identity is satisfied upon neglecting the qdependence of the form factors in q µ T µν,B Kπ . With the same diagrams and keeping in mind that we must apply the equation of motion to the charged weak vertex so that it vanishes 5 In principle the photon can also couple to K 0 due to its non-zero charge radius, so F K 0 em (q 2 ) = 0 when q 2 = 0. However, a simple ChPT calculation at O(p 4 ) indicates that |F K 0 em (q 2 )| < 0.02 when |q 2 | < 0.1 GeV 2 (see, e.g., Ref. [59]), so to our required precision it is completely negligible. On the other hand, F π 0 em (q 2 ) is exactly zero due to G-parity.
exactly when M K = M π (see, e.g. the discussion in Sec. 7 of Ref. [43]), the Born contribution to Γ µ Kπ reads, that depends on the scalar but not the vector charged weak form factor.
The Born contributions of Eqs. (27) are defined in terms of the model-independent pole contributions supplemented with a minimal seagull term required by gauge invariance. It is easy to see that if rearranging Eqs. (27) into two separately gauge invariant structures (clearly reminiscent of the usual inelastic structure functions, (−g µν +. one finds that only the contribution to F 2 contains a pole and is model-independent. The Born contribution to F 1 is regular and cannot in principle be distinguished from other inelastic contributions, so that Eqs. (27) represent the minimal Born contribution definition only, bearing residual model dependence. Fortunately, its effect on the loop integrals turns out to be very small. In δM 2 + δM a γW int , it only contributes to δf Kπ − , whose effect in the decay rate is further suppressed by r e ≈ 10 −6 (which is yet another reason why we restrict ourselves to K e3 throughout this study), whereas the contribution to δM b,V γW vanishes trivially due to symmetry.
Starting from O(p 4 ) one expects new structures such as p µ p ν /Λ 2 to enter, which parametrize inelastic contributions. Observe that a new mass scale Λ is present for dimensional reasoning, and an obvious choice is the mass of the lowest resonances. This means we are able to get a handle of the effect of the inelastic contributions by computing the contributions from the resonances at low energy. We perform that calculation based on the framework of resonance chiral theory (this is fine, as we are only dealing with tree graphs, see details in Appendix B), and find that their contribution to δ K e3 through δM 2 + δM a γW int and δM b,V γW is smaller than 10 −4 , which indicates that this contribution is negligible. However, to stay on the safe side, we introduce a common uncertainty of 2 × 10 −4 , which is roughly four times the magnitude of the resonance contribution estimated in Appendix B, to δ K e3 as a very conservative estimation of the effects from the neglected inelastic terms.
Before proceeding directly with the numerical calculations, we prefer to further isolate a particularly important piece from T µν, B Kπ and Γ µ, B Kπ known as the "convection term" [60], It corresponds to taking the contribution of the point electric charge in the Born term.
This contribution contains the full IR-divergent structure and is numerically the largest.
Being q -independent, it leads to a contribution to the loop integrals that does not depend on the specific parameterization of the hadronic form factors. Therefore it gives rise to the so-called "model-independent" contribution emphasized in Refs. [30][31][32][33], which is more commonly known as the "outer correction" in the case of free neutron and nuclear beta decays [61,62]. We thus choose to split the full Born contribution to δf Kπ + into three pieces as follows: The first and the second piece on the right-hand side of the equation above represent the IR-divergent and IR-finite contributions from the convection term, respectively. The last piece, δf Kπ + Born−conv , represents the difference between the full Born contribution and the convection term contribution. In what follows we provide the analytic results for the first two pieces: and The variables x s , x u and the loop functions are defined in Appendix C. Notice that one needs to substitute m 1 = M π , m 2 = m e , and v = s = (p + p e ) 2 in the C-functions for the for the case of K + e3 . Next, we shall study δf Kπ + Born−conv , which is the only piece that requires a specific parameterization of the hadronic form factors in order to perform the loop integral. Our first observation is that the Born contribution to δf Kπ + is UV-finite even without the form factors (it is UV-divergent for δf Kπ − without the form factors, which is however irrelevant for K e3 ). Therefore, we expect the effect of the form factors to receive a regular power suppression instead of a logarithmic enhancement.
There are different ways of parameterizing the form factors which are practically indistinguishable in the region q ∼ p e ∼ p − p ∼ M K − M π relevant to the integrals. However, in practice a simpler parameterization allows for a more straightforward evaluation of the loop integrals. Therefore, in this work, we shall adopt the monopole representation for both the electromagnetic and charged weak form factors. It is advantageous because the monopole resembles an ordinary propagator, so the q -integral reduces to standard Passarino-Veltman loop functions which can be integrated numerically with respect to {y, z}. 6 For the electro-magnetic form factors, we have: where R 2 π and R 2 K are the mean-square charge radius of π − and K + , respectively 7 . For the former, we use the result in Ref. [65]: because it was obtained through an experimental fit to the monopole form factor, which is what we adopt in this work. This value is consistent with the more recent determinations [66,67] as well as the PDG average [9], and the 2% experimental uncertainty is completely negligible in our analysis. The kaon mean-square charge radius, on the other hand, was measured with a 15% uncertainty [68]: which agrees with monopole-SU(3) estimates (see, e.g., Ref. [93]). We will include this uncertainty later in our error analysis. Finally, for the vector and scalar charged weak form factor, the monopole parameterization reads: where the fitted vector and scalar pole masses are [50]: The uncertainties are less than 5% and can be safely neglected in our analysis.
To end this section, we summarize in Table I the numerical contributions to δ K e3 from the different pieces in Eq. (32) (except δf Kπ + II that we need to combine with other terms to achieve IR-finiteness). For the error analysis, we retain only the uncertainties of the order 10 −4 or larger which, in this case, only arise from R 2 K . The first column represents the physical results, but we also consider two other cases for comparison. In the second column, we retain only the O(e 2 p 2 ) contributions, which corresponds to takingf + =f 0 = 1 and Comparing to the numbers in the first column, we find the inclusion of form factors has a larger impact on the δ K + e3 than on the δ K 0 e3 channel. In fact, the amount of shift in the former exceeds the estimated O(e 2 p 4 ) uncertainty of 0.19% in the ChPT analysis [36]. This is understandable because the effect of the form factors scales where M 2 is the typical mass scale in the monopole parameterization, and M i is the mass of the charged meson. In K + e3 we have M i = M K so the numerical impact is larger. Finally, in the third column, we consider an unphysical case where M K = 1.1M π .
We observe in this case that (δ K e3 ) fin conv (δ K e3 ) Born−conv , which proves our previous assertion that the contribution from the convection term dominates when the initial and final hadronic states are nearly degenerate.
V. VIRTUAL CORRECTION: δM 3 Next, we study δM 3 , namely the "three-point function" correction to the charged weak form factors. It was suggested in Ref. [43] to calculate such contributions in fixed-order ChPT, and we obtain the following results at O(e 2 p 2 ): where the IR-finite pieces read: The parameter Z ≈ 0.8 represents the short-distance electromagnetic effects that causes the LECs in the chiral Lagrangian with dynamical photons [37]. Finally, the loop functionsh P Q (t) are defined in Appendix A of Ref. [43].
The strategy above has a caveat, namely: There is an IR-divergent piece in δf Kπ +,3 that is numerically large, so its associated O(e 2 p 4 ) uncertainty can also be significant. Fortunately, it is straightforward to resum the IR-divergent piece to all orders in the chiral power counting by appropriately putting back the charged weak form factors based on two simple criteria as follows: 1. The combination ln(M 2 i /M 2 γ ) − 5/2 originates from the convection term contribution and should stay intact after the resummation. This is apparent by noticing that the same combination appears also in δf Kπ + II .
2. As we will show in Section VII, the IR-divergent piece from the bremsstrahlung contribution takes the following form: (the definition of δ|M | 2 brem is given in Eq.(A13)) where β i (0) is the speed of the positron in the rest frame of the charged meson (i.e. π − in K 0 e3 and K + in K + e3 ). The M γdependence above must be canceled exactly by the corresponding M γ -dependence in δf Kπ + I , δf Kπ + II and the IR-divergent piece in δf Kπ +,3 .
The arguments above lead straightforwardly to the following expression for δf Kπ +,3 : where the fully-resummed IR-divergent terms read: The last piece of the virtual corrections to f Kπ + (t) comes from δM b,A γW , which is fundamentally different from those we studied in Section IV and V in the sense that it probes the strong interaction physics in T µν Kπ from Q 2 ≡ −q 2 = 0 all the way up to Q 2 ∼ M 2 W . At large Q, one could perform a leading-twist, free-field OPE that gives us the large electroweak logarithm, but this treatment breaks down at small Q. Also, due to parity, there is no Born contribution in δM b,A γW that can be easily accounted for as in the previous two sections. Instead, one needs to deal with contributions from inelastic intermediate states residing at Q ∼ Λ χ that are governed by non-perturbative QCD. In the language of ChPT, their corresponding uncertainties are buried in the poorly-constrained LECs X 1 andX phys 6 [34][35][36]46].
As we mentioned in the Introduction, an important breakthrough happened in early 2020 as lattice QCD started to pick up its role in this subject. A series of first-principles calcu-lations were performed to study the so-called "forward axial γW -box" defined as follows: where φ i and φ f are two degenerate hadrons with mass M , and carry the same external momentum p φ , and F if + (0) is the form factor f if + (0) multiplied by the appropriate CKM matrix element. The first calculation of V A γW (π + , π 0 , M π ) in Ref. [44] led to the reduction of the RC uncertainty in the pion semileptonic decay by a factor of three. Shortly after that, a new calculation of V A γW (K 0 , π − , M π ) in the flavor SU(3) limit was performed [47] following the suggestion in Ref. [46]. These two calculations together provided an improved determination of the LECs X 1 andX phys 6 that agrees with the values quoted in the earlier ChPT papers [34,36,71] within error bars, which suggests that the error assignment in the latter is reasonable. However, in the pure ChPT representation, the major source of theory uncertainty in the long-range electromagnetic corrections to K e3 comes from O(e 2 p 4 ) instead of the LECs. Therefore, the significance of the calculations above was not fully revealed within the traditional framework.
In this section, we will demonstrate how the above-mentioned lattice QCD results play a decisive role within the new theory framework, namely to pin down δM b,A γW . We start by splitting the forward axial γW -box into two pieces: which come from the integral in Eq.(45) at Q 2 > Q 2 cut and Q 2 < Q 2 cut , respectively, where Q cut is a scale above which the leading-twist, free-field OPE is applicable. Throughout this work we choose Q 2 cut = 2 GeV 2 , in accordance with the original lattice QCD paper [44] 8 . The first term, V A> γW , contains a large electroweak logarithm and is independent of the external states {φ i , φ f } as well as the mass M . It is given by: where "+ . . ." denotes the pQCD corrections, which are at present calculated to O(α 4 s ) [72], 8 The validity of this choice is justified by the observation that the difference between the pQCD corrections to O(α 3 s ) and to O(α 4 s ) is negligible above 2 GeV 2 [45], which demonstrates the convergence of the perturbative series.
To proceed further, we perform the same splitting to the integral in δM b,A γW : where T Kπ µν A represents the component in T Kπ µν that involves the axial charged weak current. The contributions from these two terms to δf Kπ + are denoted as δf Kπ γW , respectively, and will now be related to the different components of the forward axial γWbox. First, since at Q 2 > Q 2 cut |p e | 2 we can set p e → 0 in the integrand, one can show using OPE that, Adding this piece to (δf Kπ + ) I in Eq.(24) reproduces the full electroweak logarithm in the total RC. Next, we can parameterize the integral at Q 2 < Q 2 cut as: so it is obvious that: To relate this quantity to the recent lattice QCD results, we set p → p and p e → 0 on both sides of Eq.(50). That gives 9 : 9 In the last line we made two implicit approximations: (1) we do not distinguish the value of f Kπ + (0) between the case of M K > M π and M K = M π , and (2) we add the t-dependence to the form factor. Both approximations only lead to changes of a few percent in f Kπ + , which is completely negligible after multiplying with V A< γW (K, π, M π ).
Since the lattice community has computed V A< γW (K, π, M π ), we can obtain g + (M 2 π , M 2 π , 0, M 2 π , M 2 π ) which is not exactly the same as g + (M 2 K , M 2 π , m 2 e , s, u) that we seek. However, remember that the integral in Eq. (50) is dominated by the physics at the scale q ∼ Λ χ (e.g. Regge physics [45]), it is then possible to simply take g + (M 2 π , M 2 π , 0, M 2 π , M 2 π ) together with an appropriately-assigned uncertainty: where E is an energy scale that characterizes the non-forward (NF) kinematics in Eq.(50), e.g. M K − M π , (s − M π ) 1/2 or (u − M π ) 1/2 . Since they are all smaller than M K , we can take E → M K as a conservative estimation of the uncertainty due to the NF effects. So, combining Eqs. (49), (51) and (53), we obtain: Notice that only the term in the square bracket is associated to an O(M 2 K /Λ 2 χ ) uncertainty. The recent lattice calculations provided the forward axial γW -box in the charged pion and neutral kaon decay: The box diagram in charged kaon decay is not yet computed, but can be related to the first two through a matching to the O(e 2 p 2 ) ChPT expression: V A< γW (K + , π 0 , M π ) = 2 V A< γW (π + , π 0 , M π ) − V A< γW (K 0 , π − , M π ) = 1.064(71) lat × 10 −3 . (56) The higher-order ChPT corrections to the expression above scales as O(M 2 π /Λ 2 χ ) and can be safely neglected in our error analysis 10 . With the numbers above, we obtain the numerical correction to the K e3 decay rate from δM b,A γW , as summarized in Table III. Notice that the NF uncertainty is obtained by simply multiplying 2 V A< γW (K, π, M π ) with M 2 K /Λ 2 χ . To end this section, we briefly discuss the future role of the lattice QCD. The estimation of the NF uncertainty in Eq.(53) is physically sound but can be further improved with an 10 Nevertheless, a direct lattice calculation of V A γW (K + , π 0 , M π ) in the future is still very much desirable as it provides an excellent test of the convergence speed of the chiral expansion in the SU(3) limit.   extra lattice calculation. This can be seen by considering the following relations: Both equations are obtained through a matching between the calculation of the RC based on Sirlin's approach and ChPT; the first line was given in Ref. [46] and the second line can be derived accordingly. We see that both V A γW (K 0 , π − , M π ) and V A γW (K + , K 0 , M k ) are matched to the same combination of LECs, except that the latter is subject to larger higherorder corrections because the involved meson mass is M K which is larger. That means, the difference in the numerical values between V A γW (K 0 , π − , M π ) and V A γW (K + , K 0 , M k ) provides an estimation of the size of the NF corrections in Eq.(53). This strategy is very similar to the standard lattice QCD technique to estimate the size of the chiral power corrections through the variation of the quark masses.

VII. BREMSSTRAHLUNG CONTRIBUTION
After going through all the virtual corrections, we switch to the contribution from the n = 1 process, which is simply known as the "bremsstrahlung contribution". According to the discussions in Appendix A, the bremsstrahlung process contributes to the differential decay width dΓ K e3 /dydz not only in the D 3 region but also in the D 4−3 region, the latter has no correspondence in the n = 0 process. Therefore, it is eventually up to the experimentalists to decide in which region of {y, z} will the data be taken, and whether or not a veto will be applied to exclude decay events with hard photons. Of course, the simplest choice is to not apply any veto, and to collect data from all available regions of {y, z}. This corresponds to a fully-inclusive prescription of the real photon emission process, or in other words, we should calculate the sum of the full n = 0 and n = 1 decay width. This prescription was adopted in Ref. [36] and will be followed in this work.
The bremsstrahlung amplitude, depicted by the two diagrams in Fig.2, reads: We observe that the generalized Compton tensor T Kπ µν appears again, only that now one deals with a real photon. Unlike in the loop diagrams, here we only need to know T Kπ µν for small (due to the phase-space constraint) and on-shell photon momentum k, so instead of exhausting the contributions from all intermediate states, it is possible to adopt a low-energy effective expression T Kπ µν . It should, however, satisfy three basic criteria: • It must contain the full convection term contribution to ensure an exact cancellation of the IR-divergence from the virtual corrections.
• It should include the seagull term, as the effect of the latter is not particularly suppressed in the decay rate, unlike in the loop diagrams.
• It should satisfy exact electromagnetic gauge invariance, so that one could perform the usual replacement s ε µ s (k)ε ν * s (k) → −g µν in the sum of the outgoing photon polarizations.
The simplest effective expression that satisfies all these criteria is: The first term on the right-hand side in the expressions above is just the convection term, whereas the remainders are the seagull term and the extra pieces from the Born contribution needed to recover gauge invariance. Notice that the convection term is exact, and only the terms in the curly bracket undergo a chiral expansion. In fact, if we expand the convection term to O(p 2 ), the LO ChPT expression in Eq.(25) is recovered. In fact, the existing ChPT calculation uses exactly Eq. (25) in their calculation of the bremsstrahlung effect, but now our expression allows a resummation of the most important terms in T µν Kπ to all chiral orders. With the above, the bremsstrahlung amplitude splits into two pieces: and for K + → π 0 e + ν e γ, The significance of such a splitting is that M A is an exact expression and only M B involves a chiral expansion. Therefore, in the computation of the decay rate, only the contribution from 2Re{M * B M A } + |M B | 2 acquires an O(e 2 p 4 ) uncertainty, while the contribution from |M A | 2 is exact. As we will show later, this brings an advantage over the existing treatment as the latter is numerically the largest. Now we proceed to the phase space integration of the bremsstrahlung contribution. We first discuss the integration in the D 3 region. To isolate the IR-singular term, we first split |M A | 2 into two pieces: where p i = p (p ) in K + e3 (K 0 e3 ). The integration of the first term with respect to { p ν , k, x} produces an IR-divergence: where the explicit expressions of I IR i and I fin i can be found in Appendix D, and with this, we verify our previous assertion about the IR-divergent structure of the bremsstrahlung contribution in Eq. (42). We can now combine the IR-divergent contributions from the virtual corrections (which we previously labeled as I, II, III) with the bremsstrahlung contribution in the D 3 region to obtain the following shift of the K e3 decay rate: where which is now explicitly IR-finite. We observe that the expression above still contains a residual integral with respect to { p ν , k, x}, but it is IR-finite and therefore can be straightforwardly carried out with the method outlined in Appendix E. The numerical result is summarized in Table IV. The HO uncertainty comes from δ QED HO , while the O(e 2 p 4 ) uncertainty is obtained by multiplying the contribution from 2Re We see that these uncertainties are as small as 10 −4 , which is a clear success of our strategy in the splitting of T µν Kπ (k ; p , p) in Eq. (59). Finally, we also need to compute the bremsstrahlung contribution in the D 4−3 region: The integrals are IR-finite and can be carried out similarly using the method in Appendix E.
The numerical results are given in Table V. In principle one also acquires an O(e 2 p 4 ) un-   certainty by multiplying the contribution from 2Re {M * A M B } + |M B | 2 by M 2 K /Λ 2 χ , but the outcomes are of the order 10 −5 and so are not displayed in the table.

VIII. COMPARING WITH THE CHPT RESULT
We have now finished calculating all components of the O(G 2 F α) electroweak RC to the K e3 decay rate. The total result is simply given by: where the numerical values of different components can be found in Tables I-V. On the other hand, in the existing standard ChPT treatment the full electroweak RC is broken down into "short-distance" and "long-distance" pieces, and are allocated to several different quantities, some of which are somewhat implicitly hidden. This section serves to perform a rigorous matching between our result and the values quoted in the existing ChPT literature, with special attention paid to the so-called "long-distance electromagnetic corrections" δ Ke EM . In the standard ChPT framework, the photon-inclusive K e3 decay rate is parameterized as [9]: where C K is a simple isospin factor. Apart from the quantity |f K 0 π − + (0)| that requires a lattice input, all the small QCD and electroweak corrections to Γ K e3 are distributed into the following four quantities: S EW , I Ke (λ i ), δ Kπ SU (2) and δ Ke EM . We shall take a serious look at each of these quantities, and study their relations to the different components of electroweak RC we calculated in this work.

A. S EW
The quantity S EW was first introduced by Marciano and Sirlin in Ref. [73] as a processindependent factor that accounts for the large electroweak logarithm in the electroweak RC [41,74] including the O(α s ) pQCD corrections on top of it, as well as the resummation of the QED logs (i.e. δ QED HO in our notation). It was often quoted schematically in the literature as [36,75]: where the ρ-mass appears as a low-energy scale. It is not straightforward to infer its exact value from the expression above because some of the important components (e.g. δ QED HO ) are not explicitly shown, and it is also not clear what scale one should choose for α s . Fortunately, as a common consensus, the value S EW = 1.0232(3) HO was always used for all practical purposes in the recent years (see, e.g. Refs. [34,35] and the FLAVIAnet global analysis, Ref. [76]), where the central value comes from Ref. [73] and the estimated uncertainty of the QED log resummation comes from Ref. [51]. Notice that although Ref. [75] quoted a slightly different value of S EW = 1.0223(5), but that number was never used in any subsequent analysis. Now, the process-independent physics included in our (δ K e3 ) tot are not only those described by S EW but even more. For example, the most important pQCD correction con- Therefore, it is not the most natural choice to remove S EW − 1 analytically from (δ K e3 ) tot in order to compare our result with the ChPT result. Instead, it is more convenient to take the above-mentioned numerical value of S EW simply as its definition, i.e., and remove this value numerically from (δ K e3 ) tot for the comparison. This prescription keeps us on the same track with all the recent literature mentioned above.
The quantity I (0) Ke (λ i ) is formally defined as the "phase space integral depending on slope and curvature of the form factors f Kπ ± (t)" according to Ref. [75], but in practice it is treated not just as a pure QCD factor, but also contains a part of the short-distance electromagnetic effects. This can be seen in, e.g., Refs. [34,35]: The t-dependence of f Kπ ± (t) at O(p 4 ) is given by the mesonic loop functions H P Q (t), and we observe that in these functions the masses of the charged mesons (e.g. π ± ) and their neutral counterparts (e.g. π 0 ) are kept distinct.
Since we know that this mass splitting is partially induced by short-distance electromagnetic effects, or more specifically, the O(e 2 ) term in the chiral Lagrangian [38]: so the observation above implies that a part of the short-distance electromagnetic effect proportional to Z is actually assigned implicitly to I Ke (λ i ) through H P Q (t) within the ChPT framework. In our notation, this residual effect is represented exactly by theh P Q (t) terms in fin e 2 p 2 , since theh P Q (t) functions are simply consequences from the Taylor expansion of H P Q (t) to O(Z).

SU(2)
The isospin-breaking correction factor δ Kπ SU(2) is formally defined as 11 : that is only present in K + l3 . According to the definition above, it contains not only the strong isospin breaking effect resulting from the u-d mass difference, but also the electromagnetically-induced isospin breaking. Indeed, according to Eq.(4.42) in Ref. [75], one has: where ε EM originates from the electromagnetically-induced π 0 -η mixing. In our notation, this correction simply comes from δf K + π 0 +,3 (t) fin e 2 p 2 after removing theh P Q (t) terms. 11 The existence of the isospin factor C K 0 /C K in the formula above is simply due to our choice of normal-

EM
After all the discussions above, it is now apparent that the most convenient way to discuss δ Ke EM is to simply refer it as "the sum of all electroweak RC that are not already contained in S EW , I Ke (λ i ) and δ Kπ SU(2) ". This means in our notation, where S EW − 1 is defined by Eq.(70) as we discussed earlier. Apart from S EW −1, the quantity (δ K e3 ) fin 3 is also subtracted out because its contribution is redistributed into I (0) Ke (λ i ) and δ Kπ SU (2) according to the ChPT prescription, as we discussed above. In fact, δ Ke EM is also the only meaningful quantity to be compared between this work and the existing literature, because we are taking an O(e 2 p 2 ) approximation to (δ K e3 ) fin 3 and thus have made no new improvement on this term.
The comparison between our result of δ Ke EM and the ChPT result is given in Table VI. Let us explain all the different types of uncertainties that appear in our new evaluation: • inel: This represents our conservative estimation of the effects from the inelastic term in δM 2 + δM a γW int and δM b,V γW . See the discussions after Eq.(29).
• R 2 K : This is the uncertainty originated from the experimental error of the K + charge radius (see Eq.(37)) that enters δM 2 + δM a γW int and δM b,V γW in K + e3 .
• lat: This is the total lattice QCD uncertainty in the calculation of V A γW (see Eq.(55)).
• NF: This represents our estimation of the uncertainty due to the non-forward kinematics in δM b,A γW at small loop momentum q . We include an asterisk to remind the reader that this error estimation can be made more rigorous with an extra lattice QCD calculation, as we discussed at the end of Section VI.
• e 2 p 4 : This is the chiral expansion uncertainty of the non-convection term contribution (i.e. 2Re {M * B M A } + |M B | 2 , see the discussions after Eq.(61)) in the bremsstrahlung process.
From Table VI we find that our results are consistent with the ChPT estimation within the error bars, but with a significant reduction of the total uncertainty by almost an order of magnitude. This improvement is mainly due to two reasons: 1. Our calculation permits a much better control of the O(e 2 p 4 ) effects, which are the main source of uncertainty in the ChPT treatment. With the new theory framework introduced in Refs. [43,46], all the hadron physics are contained in quantities such as Kπ and Γ µ Kπ , from which the full convection/Born contribution can be explicitly isolated. These contributions govern the full IR-divergent structure of the decay process, are numerically the largest and, most importantly, do not involve any chiral expansion.
The size of the non-Born/non-convection term contributions are in general an order of magnitude smaller (see, for example, Table IV  2. We used latest lattice QCD results to pin down δM b,A γW , which corresponds to the LECs X 1 andX phys 6 in ChPT. In the existing literature, these LECs were calculated within resonance models and were assigned a 100% uncertainty. On the other hand, the highly-precise lattice results of V A γW would correspond exactly to δM b,A γW if K and π were degenerate. We investigated the region of integration in δM b,A γW where this non-degeneracy starts to take effect, and assigned a reasonable NF-uncertainty to the contribution from this region on top of the lattice results. In the ChPT language, our treatment above simultaneously take into account the uncertainties of the LECs themselves as well as the O(e 2 p 4 ) uncertainties on top of the LEC contributions.

IX. FINAL DISCUSSIONS
The 3σ discrepancy in the extraction of V us from K µ2 and K l3 decays has triggered renewed interest within the particle physics community about its possible implications on the existence of BSM physics. However, the current level of significance is not sufficient to claim a discovery so one needs further reduction of not just the experimental errors but also the SM theory uncertainties. Our re-analysis of the SM electroweak RC in K e3 therefore, serves as a crucial step along this direction. We successfully overcome the natural limitations in traditional ChPT by adopting a new computational framework that allows for a resummation of the numerical largest components in the RC, and also utilizing the most recent lattice QCD outcomes to reduce the uncertainties from the non-perturbative QCD at the chiral symmetry breaking scale. Our work reduces the existing uncertainties in the K e3 RC by almost an order of magnitude, and finds no large shift in the central values. This suggests that we should remove the electroweak RC from the "list of culprits" responsible for the K µ2 -K l3 discrepancy.
Is it evident now that the V us anomaly cannot be explained by SM effects? We would say that it is still too early to decide at this stage. Further investigations must also be made on other SM inputs, just to mention a few: • Based on the analysis of a newly-constructed ratio R V = Γ K e3 /Γ π e3 , Ref. [77] suggested that a shift of the lattice QCD input of |f K 0 π − + (0)/f π + π 0 + (0)| from its current value of 0.970(2) to a smaller value of 0.961(4) would reconcile the K µ2 and K l3 results, and encouraged the lattice community to examine this possibility. Lattice calculations of [78,79] and N f = 2 + 1 + 1 [80][81][82] in the recent years have so far been consistent with each other, which led to the FLAG 2019 averages [12]: However, a new calculation by the PACS collaboration with N f = 2 + 1 returned (19)(1) that is significantly lower than the existing average [83]. This calculation utilized only one lattice spacing a = 0.085 fm and thus should be carefully reexamined.
• The quantity I (0) Kl (λ i ) probes the t-dependence of the form factorsf +,0 (t). Adopting a Taylor-expansion parameterization: the parameters λ +,0 and λ + are fit to the experimental distributions of the K l3 decays to obtainf +,0 (t) in the physical region of t. The resulting uncertainties are 0.13% for I (0) Ke and 0.31% for I (0) Kµ (see Table 21 in Ref. [84]), which look well under control; other forms of parameterization were also investigated [85][86][87][88][89]. However, it is known for some time that some disagreements occur in the extracted values of the slope parameter λ 0 of the scalar form factor from different experiments [75]. Also, sincē f +,0 (t) are pure QCD quantities, their fitting to the K l3 distributions can only be done after removing the effects of the electroweak RC from the experimental data. Now since we have updated the RC analysis, the fitting procedure should in principle also be updated accordingly. Although in this paper we only present our updates of δ Ke EM , but the electromagnetic corrections to the K e3 Dalitz plots can also be derived with the same method.
• Although the SU(2) isospin breaking correction factor δ Kπ SU (2) exists only in the K + channel by construction, its associated theory uncertainty is the largest. Upon neglecting the electromagnetic contributions, it is given by: in ChPT to O(p 4 ), where Q 2 ≡ (m 2 s −m 2 )/(m 2 d − m 2 u ) ≡ R(m s /m + 1)/2 and χ p 4 = 0.219 [90]. The main uncertainties therefore come from Q and R. For instance, disagreements are observed between the values of Q and R extracted from phenomenology [91] η → 3π : Q = 22.1(7), R = 34.4(2.1) (78) and from lattice QCD [12] which must be sorted out in order to pin down the isospin breaking correction precisely.
Finally, we want to mention that we present in this work only our updates on the electroweak RC but not a new value of V us . A part of the reason is that we work exclusively on K e3 and not on K µ3 , given that the latter involves more sources of uncertainty (e.g. from δf Kπ − ) and will be a subject of future study. But more importantly, we realize that the physics of kaon decay is a dynamically progressing field from where the knowledge in both experiment and theory, including our understanding of the issues above, is being constantly updated. Therefore, rather than quoting a new value of V us upon every single improvement, it is more preferable to have a commonly agreed value that results from a collaborative work between experimentalists and theorists based on the most updated inputs from their respective fields, similar to the FLAVIAnet evaluation in the past decade [76]. We hope that our research may serve as a useful input for a possible future collaboration of such kind.
Note added: Awaiting the review outcome of this manuscript, some of us published a new global analysis of V us from K l3 based on the improvements in this work [92]. The values of |V us | extracted from K e3 and K µ3 are currently consistent with each other within error bars, therefore we do not see a noticeable violation of lepton flavor universality within K l3 . This requires further check from theory improvements of the K µ3 RC as well as future experiments.
Without performing the integral, one already sees that the δ-function imposes the constraint x ≥ 0 because P 2 = M 2 K x = (k + p ν ) 2 is just the invariant squared mass of the νγ system, which cannot be negative. With that one splits the x-integral into two terms: dx , (A8) and the different step functions in front of each term impose different constraints on the integration region of {y, z}. The first term requires α − (y, z) < 0 < α + (y, z), which simply gives the D 3 region we discussed above. Meanwhile, the second term requires α − (y, z) > 0, and solving this inequality yields a different integration region which we may call D 4−3 . It can again be represented in two equivalent ways: There is no overlap between the region D 3 and D 4−3 (see Fig.3). With the above, the K → πe + νγ decay rate can be written as: In the study of a fully-inclusive kaon semileptonic decay rate up to O(G 2 F α), one should add the K → πe + ν and K → πe + νγ decay rate to give: where Both |M | 2 K→πe + ν and δ|M | 2 brem possess IR-divergences that eventually cancel other. Meanwhile, the term with the integration over the D 4−3 region is by itself IR-finite.

Appendix B: Resonances at low energy
In this Appendix, we briefly review the basics of the resonance chiral theory that includes the 1 −− and 1 ++ resonances as dynamical DOFs in the chiral Lagrangian [93][94][95]. Based on this formalism we calculate the contribution of these resonance to δM 2 + δM a γW and δM b,V γW . In most of the literature on resonance chiral theory, the massive spin-1 particles are described by a totally-antisymmetric tensor field instead of a vector field [96], so we start by introducing the formalism. First, the free Lagrangian of a (real) massive spin-1 particle is written as: where R µν is the antisymmetric tensor field. It satisfies the following classical equation of motion: The quantized field takes the form: where ε s ( k) is the polarization vector of the spin-1 particle that satisfies the following rela-tions: andâ + s ( k),â s ( k) are the creation and annihilation operators. Finally, by inverting the free Lagrangian one obtains the covariant propagator of the antisymmetric tensor field: We can now construct the chiral Lagrangian with dynamical vector and axial resonances.
The 1 ++ octet resonances are represented by V µν which is a traceless, Hermitian matrix in the flavor space. Its chiral covariant derivative is given by: where is the standard connection vector, with v µ , a µ the vector and axial external sources. Similarly, the 1 −− resonances are represented by the matrix A µν . Other elementary building blocks of the ordinary ChPT. include the "vielbein": and the anti-symmetric tensors f µµ R,L built from the vector and axial external sources: and finally, f µν ± ≡ uf µν L u † ± u † f µν R u. With the above we can now write down the chiral Lagrangian with 1 ++ and 1 −− resonances. The LO Lagrangian scales as O(p 4 ): where ... represents the trace over the flavor space, M V and M A are the vector and axial resonance masses in the chiral limit, while F V , F A and G V are real coupling constants. The leading resonance contribution to T µν Kπ scales as O(p 4 ) and enters through the sand u-channel diagrams as depicted in Fig.4. Since all the couplings in Eq.(B10) have even intrinsic parity, it is evident that only the axial resonances can exist in the intermediate state. They give rise to the following expressions: where F 0 is the pion decay constant in the chiral limit. For numerical estimation, we choose F A = 123 MeV, M A = 968 MeV following Ref. [93], and F 0 ≈ F π = 92.1 MeV.
Meanwhile, since Γ µ Kπ vanishes in the flavor SU(3) limit, it cannot be generated by the resonance Lagrangian in Eq.(B10) at tree level because the latter is SU(3)-symmetric.
We then plug the expressions above into Eq. (22), (23) and evaluate the integrals. Of course, upon setting M 2 W /(M 2 W − q 2 ) → 1 the integrals are UV-divergent, but this is expected because the expressions above are only supposed to work at small q so the integral should be cut off at q ∼ M A . As our main purpose here is just to have an order-of-magnitude estimation of the resonance contribution, we shall adopt a simple prescription as follows: we first regularize the UV-divergence using dimensional regularization, and discard the usual divergent combination 2/(4 − d) − γ E + ln 4π. The result is then a function of the renormalization scale µ, which we vary from M A to 2M A as a crude estimation of the uncertainty.
With the above, we obtain the following resonance contribution to δ K e3 : They are both smaller than 10 −4 .

(D5)
And finally, one expands the result to O ((d − 4) 0 ). It is also customary to switch the result back to the expression with the M γ -regularization. For that purpose one simply performs the following matching: Next, we discuss some useful tricks in the evaluation of I i (y, z) with dimensional regularization. First, the full integral can be split into three terms, with the integrand proportional to: 1 (p e · k) 2 , respectively. The integration with respect to the first term is most easily done in the p e -rest frame, while the next two terms should be done in the p i -rest frame. The following identity is also useful in performing the integration of the third term: We are now ready to write down the full result of the integral: I i (y, z) = I IR i (y, z) + I fin i (y, z) , where is the IR-divergent piece after switching back to the M γ -prescription using Eq.(D6), and is the IR-finite piece, with The correct analytic expression for I K (y, z) and I π (y, z) first appeared in Ref. [25] and Ref. [35] respectively (notice that Ref. [24] also attempted to calculate I π (y, z), but the result there is wrong even with the Errata). It is easy to check the numerical equivalence between Eq.(D9) and those expressions, after accounting for the difference in the overall normalization.
Appendix E: IR-finite integrals in the bremsstrahlung contribution In this Appendix, we outline the general strategy to evaluate the IR-finite numerical integrations from the bremsstrahlung process, in both the D 3 and D 4−3 region. We start by providing the expressions of the relevant integrands. In K 0 in Mathematica such as Tracer or Package-X. After taking the trace, all the expressions above are functions of {x, y, z} as well as two of the three following dot products involving k: {k · p, k · p , k · p e } using the identity 2k · (p − p − p e ) = M 2 K x. The integration can be performed with the following strategy. Take |M A | 2 res in K 0 e3 as an example: we first express the squared amplitude as a finite sum: where −2 ≤ m, n ≤ 2 and c m,n (x, y, z) are known scalar coefficients. The p ν and kintegrations return the following functions: of which analytic expressions are given in the Appendix of Ref. [25] (we have checked their correctness). With this, we obtain: m,n c m,n (x, y, z)I m,n (p , p e ) , (E5) where the right-hand side is now a function of {x, y, z}, so the remaining three-fold integration with respect to these variables are completely regular and can be performed numerically.
The same strategy applies to the IR-finite integrals in K + e3 , except that one should choose 1/ {(k · p) m (k · p e ) n } as the basis.