First lattice calculation of the QED corrections to leptonic decay rates

The leading-order electromagnetic and strong isospin-breaking corrections to the ratio of $K_{\mu 2}$ and $\pi_{\mu 2}$ decay rates are evaluated for the first time on the lattice, following a method recently proposed. The lattice results are obtained using the gauge ensembles produced by the European Twisted Mass Collaboration with $N_f = 2 + 1 + 1$ dynamical quarks. Systematics effects are evaluated and the impact of the quenched QED approximation is estimated. Our result for the correction to the tree-level $K_{\mu 2} / \pi_{\mu 2}$ decay ratio is $-1.22\,(16) \%$ to be compared to the estimate $-1.12\,(21) \%$ based on Chiral Perturbation Theory and adopted by the Particle Data Group.


I. INTRODUCTION
The determination of a number of hadronic quantities relevant for flavour physics phenomenology using lattice QCD simulations has reached such an impressive level of precision [1] that both electromagnetic (e.m.) and strong isospin-breaking (IB) effects cannot be neglected.
In the past few years accurate lattice results including e.m. and IB effects have been obtained for the hadron spectrum, as in the case of the charged/neutral mass splittings of pseudoscalar (P) mesons and baryons (see, e.g., Refs. [2,3]). In this respect the inclusion of QED effects in lattice QCD simulations has been carried out following mainly two methods: in the first one QED is added directly to the action and QCD+QED simulations are performed at few values of the electric charge (see, e.g., Refs. [3,4]), while the second one, the RM123 approach of Refs. [2,5], consists in an expansion of the lattice path-integral in powers of two small parameters: the e.m. coupling α em and the light-quark mass difference (m d − m u )/Λ QCD , which are both at the level of ≈ 1%. Since it suffices to work at leading order in the perturbative expansion, the attractive feature of the RM123 method is that the small values of the two expansion parameters are factorised out, so that one can get relatively large numerical signals for the slopes of the corrections with respect to the expansion parameters. Moreover the slopes can be determined in isosymmetric QCD. In this letter we adopt the RM123 method.
While the calculation of e.m. effects in the hadron spectrum does not suffer from infrared (IR) divergences, the same is not true in the case of hadronic amplitudes, where e.m. IR divergences are present and cancel for well defined, measurable physical quantities only after including diagrams containing both real and virtual photons [7]. This is the case, for example, for the leptonic π 2 and K 2 and the semileptonic K 3 decays, which play a crucial role for an accurate determination of the Cabibbo-Kobayashi-Maskawa (CKM) entries |V us /V ud | and |V us | [8].
The presence of IR divergences requires the development of additional strategies to those used in the computation of e.m. effects in the hadron spectrum. Such a new strategy was proposed in Ref. [6], where the lattice determination of the decay rate of a charged pseudoscalar meson into either a final ± ν pair or ± ν γ state was addressed. Although it is possible in lattice simulations to compute the e.m. corrections due to the emission of real photons, this is not strictly necessary. Instead, the amplitude for the emission of a real photon can be computed in perturbation theory by limiting the maximum energy of the emitted photon in the meson rest-frame, ∆E γ , to be small enough so that the internal structure of the decaying meson is not resolved, but larger than the experimental energy resolution [6]. The IR divergences are independent of the structure of the hadrons (i.e. they are universal) and cancel between diagrams containing a virtual photon (computed non-perturbatively) and those with the emission of a real photon (calculated perturbatively). In the intermediate steps of the calculation, however, it is necessary to introduce an IR regulator. We use the lattice volume L 3 itself as the IR regulator by working in the QED L finite-volume formulation of QED [9].
The inclusive rate Γ(P 2 ) can be expressed as [6] Γ( , where the subscripts 0, 1 indicate the number of photons in the final state and the superscript pt denotes the pointlike approximation for the decaying meson. The terms Γ 0 (L) and Γ pt 0 (L) are evaluated on the lattice; both have the same IR divergences which therefore cancel in the difference. We use L as the intermediate IR regulator and Γ 0 − Γ pt 0 is independent of the regulator as this is removed [10]. Since all momentum modes contribute to it, Γ 0 (L) depends on the structure of P and must be computed non-perturbatively. In the second term on the r.h.s. of Eq. (1) P is a point-like meson and both Γ pt 0 (µ γ ) and Γ pt 1 (∆E γ , µ γ ) can be calculated directly in infinite volume in perturbation theory, using a small photon mass µ γ as the intermediate IR regulator. Each term is IR divergent, but the sum is convergent [7] and independent of the IR regulator. The explicit perturbative calculations of Γ pt 0 + Γ pt 1 (∆E γ ) and Γ pt 0 (L) have been performed in Refs. [6] and [10], respectively.
The inclusive decay rate (1) can be written as where Γ (tree) P is the tree-level decay rate given by where M P is the physical mass of the charged P-meson, including both e.m. and strong IB corrections. The superscript (0) on a physical quantity denotes that it has been calculated in isosymmetric QCD (without QED). The P-meson decay constant, f P is defined by: In Eq. (2) the quantity δR P encodes the leading-order e.m. and strong IB corrections to the tree-level decay rate and its evaluation is described in Ref. [6]. Its value depends on the prescription used for the separation between the QED and QCD corrections, while the quantity [f (0) P ] 2 (1 + δR P ) is prescription free. In this work we adopt the prescription in which the renormalized couplings and quark masses in the full theory and in isosymmetric QCD coincide in the MS scheme at a scale of 2 GeV [2,5] (see the supplemental material).
In this letter we focus on the ratio of the inclusive decay rates of kaons and pions into muons, namely where δR Kπ ≡ δR K − δR π . Using the gauge ensembles generated by the European Twisted Mass Collaboration (ETMC) with N f = 2 + 1 + 1 light, strange and charm sea quarks [12,13], we have calculated δR Kπ , which, together with a lattice computation of f π , allows us to determine |V us /V ud | from the ratio in Eq. (5).
The quantity δR Kπ is less sensitive to various uncertainties than the individual terms δR K and δR π . Three main features help to reduce the systematic uncertainties in δR Kπ : (i) In Γ 0 (L) all the terms up to O(1/L) are "universal", i.e. independent of the structure of the decaying hadron [10]. The residual, structure-dependent (SD) finite volume effects (FVEs) start at order O(1/L 2 ) and are found to be much milder in the case of δR Kπ (see section IV). (ii) The matching of the bare lattice weak operator with the one renormalised using W -regularization generates a mixing of operators of different chiralities when discretisations based on Wilson fermions, which break the chiral symmetry (such as twisted mass used here) are used. The mixing has been calculated only at order O(α em α 0 s ) [6], but its effects cancel out in the difference δR Kπ . (iii) Within SU(3) Chiral Perturbation Theory (ChPT) the effects of the sea-quark electric charges depend on unknown low-energy constants (LECs) starting at next-to-leading-order (NLO) for δR K and δR π , but only at NNLO for δR Kπ [11]. Thus, the uncertainty due to the quenched QED (qQED) approximation, adopted in this work, is expected to be smaller for δR Kπ .
Since the experimental rates Γ(K µ2 ) and Γ(π µ2 ) are inclusive, SD contributions to the real photon emission should be included. According to the ChPT predictions of Ref. [14], however, these contributions are negligible in the case of both kaon and pion decays into muons, while the same does not hold as well in the case of final electrons (see Ref. [6]). This important finding will be investigated by an ongoing dedicated lattice study on the real photon emission amplitudes in light and heavy Pmeson leptonic decays.
After extrapolating our lattice data to the physical pion mass and to the continuum and infinite volume limits, the main result of the present work is where the uncertainty includes both statistical and systematic errors, including an estimate of the uncertainty due to the QED quenching. Our result (6) can be compared with the current estimate δR phys Kπ = −0.0112 (21) from Refs. [15,16] adopted by the PDG [17].

II. DETAILS OF THE SIMULATION
The gauge ensembles used in this work were generated by ETMC with N f = 2+1+1 dynamical quarks and used in Ref. [18] to determine the up, down, strange and charm quark masses. The main parameters of the simulations are collected in the supplemental material. We employ the Iwasaki action [19] for gluons and the Wilson Twisted Mass Action [20][21][22] for sea quarks. In the valence sector we adopt a non-unitary setup [23] in which the strange quark is regularized as an Osterwalder-Seiler fermion [24], while the up and down quarks have the same action as the sea. Working at maximal twist such a setup guarantees an automatic O(a)-improvement [22,23]. The two valence quarks in the P-meson are regularized with opposite values of the Wilson r-parameter in order to guarantee that discretization effects on the P-meson mass are of order O(a 2 µ Λ QCD ). The lepton is a free twisted-mass fermion with mass m = m µ = 105.66 MeV [17]. The neutrino is simply considered to be a free fermion field.
In this work we make use of the bootstrap samplings generated for the input parameters of the quark mass analysis of Ref. [18]. There, eight branches of the analysis were adopted differing in: (i) the continuum extrapolation, adopting for the matching of the lattice scale either the Sommer parameter r 0 or the mass of a fictitious P-meson made up of two valence strange(charm)like quarks; (ii) the chiral extrapolation performed with fitting functions chosen to be either a polynomial expansion or a ChPT Ansatz in the light-quark mass; and (iii) the choice between the methods M1 and M2, which differ by O(a 2 ) effects, used to determine the mass renormalization constant Z m = 1/Z P in the RI'-MOM scheme.

III. EVALUATION OF THE AMPLITUDES
Following Ref. [6] the quantity δR K − δR π is given by where δΓ (pt) P (∆E γ ) represents the O(α em ) correction to the tree-level decay rate for a point-like meson and can be read off from Eq. (51) of Ref. [6], while δA P and δM P are the e.m. and IB corrections to the weak amplitude and mass of the P-meson, respectively.
Within the qQED approximation the evaluation of δA P and δM P requires the evaluation of only the connected diagrams shown in Figs. 1 -4 for K 2 decays. The corrections δA P and δM P can be written as where δA QCD P (δM QCD P ) represents the strong IB corrections corresponding to the diagrams of Fig. 3, while the other terms are QED corrections coming from the insertions of the e.m. current and tadpole operators, of the pseudoscalar and scalar densities (see Refs. [2,28]).
In Eqs. (8-9) the term δA J P (δM J P ) is generated by the diagrams of Fig. 1(a-c), δA T P (δM T P ) by the diagrams of Fig. 1(d-e), δA P P (δM P P ) by the diagrams of Fig. 2(a-b) and δA S P (δM S P ) by the diagrams of Fig. 3(a-b). The term δA P corresponds to the photon exchange between the quarks and the final lepton. It arises from diagrams 4(ab), while the diagram 4(c) (lepton wave function renormalization) can be safely omitted, since it cancels out exactly in the difference Γ 0 (L) − Γ pt 0 (L). The evaluation of δM QCD P and the δM i P is described in Ref. [5], where the quark mass difference (m d − m u )(MS, 2 GeV) = 2.38 (18) MeV was determined using the experimental charged and neutral kaon masses. The terms δA QCD P , δA i P and δA P are extracted from the cor-relators described in Ref. [6]. Their numerical determination is illustrated briefly in Refs. [25,26] and in detail in Ref. [27]. The quality of the extraction of δA =µ P /δA (0) P is illustrated in the supplemental material.

IV. FINITE VOLUME EFFECTS AT O(αEM )
The subtraction Γ 0 (L)−Γ pt 0 (L) makes the rate IR finite and cancels the structure-independent FVEs. The pointlike decay rate Γ pt 0 (L) is given by where the factor Y P (L) is explicitly given by Eq. (98) of Ref. [10]. Eq. (8) is therefore replaced by where Y P (L) has the form with the coefficients b j (j = IR, 0, 1, 2, 3) depending on the dimensionless ratio m /M P [10]. The important point is that the SD FVEs start only at order O(1/L 2 ), i.e. all the terms up to O(1/L) in Eq. (12) are "universal" [10]. Being independent of the structure they can be computed for a point-like charged meson.
The FVE subtraction (11) up to order O(1/L) is illustrated in Fig. 5 for δR K , δR π and δR Kπ in the inclusive case ∆E γ = ∆E max,P γ = M P (1 − m 2 µ /M 2 P )/2, which corresponds to ∆E max,K γ 235 MeV and ∆E max,π γ 29 MeV, respectively. It can be seen that after subtraction of the universal terms the residual FVEs are almost linear in 1/L 2 and ≈ 3 times smaller in the case of δR Kπ .
Using Eq. (13) we have fitted the data for δR Kπ using a χ 2 -minimization procedure with an uncorrelated χ 2 , obtaining values of χ 2 /d.o.f. always around 1.2. The uncertainties on the fitting parameters do not depend on the χ 2 -value, because they are obtained using the bootstrap samplings of Ref. [18] (see section II). This guarantees that all the correlations among the data points and among the fitting parameters are properly taken into account. The quality of our fits is illustrated in Fig. 6.
At the physical pion mass in the continuum and infinite-volume limits we obtain where () stat indicates the uncertainty induced by both the statistical errors and the fitting procedure itself; () input is the error coming from the uncertainties of the input parameters of the quark-mass analysis of Ref. [18]; () chir is the difference between including or excluding the chiral logarithm in the fits (i.e. taking R χ = 0 or R χ = 0); () F V E is the difference between including (K 2 = 0 and K 2 = 0) or excluding (K 2 = K 2 = 0) the residual FVE correction in δR Kπ [27]; () disc is the uncertainty coming from including (D = 0) or excluding (D = 0) the discretisation term proportional to a 2 ; () qQED is our estimate of the uncertainty of the QED quenching. This is obtained using the Ansatz (13) with the coefficient R χ of the chiral log fixed at the value R χ = α em (Z −3)/4π, which includes the effects of the up, down and strange sea-quark charges [11]. The change in δR phys Kπ is 0.003, which has been already added in the central value of Eq. (14). To be conservative, we use twice that value for our estimate of the qQED uncertainty.
Our result (14) can be compared to the value δR phys Kπ = −0.0112 (21) from Refs. [15,16] adopted by the PDG [17]. Using in Eq. (5) the experimental K µ2 and π µ2 decay rates [17], we obtain where the first error comes from experiments and the second one is related to the uncertainty in our result (14).
Thus, using |V ub | = 0.00413 (49) [17] the first-row CKM unitarity is confirmed to be below the per mille level, viz.

SUPPLEMENTAL MATERIAL
Details of the calculation described in this work will be presented in Ref. [27]. Here, in subsection A we collect the main parameters of the simulations performed in the isosymmetric QCD theory in Ref. [18]. These simulations correspond to a prescription for the separation between QED and QCD corrections, which is different from that of Ref. [2] adopted in this work for the calculation of δR Kπ . We show that the two prescriptions differ only by effects which are well within the uncertainties of the input parameters of Ref. [18].
In subsection B we sketch some of the key points of the numerical analysis and illustrate the quality of the results by showing the time-dependence of the most complicated diagrams, i.e. those in Fig. 4(a) and (b) in which a photon is exchanged between the quarks and the final-state charged lepton.

A. Simulation parameters
The main parameters of the simulations performed within isosymmetric QCD in Ref. [18] are collected in Table I.   [18] and for the gauge ensemble, A40.40 added to improve the investigation of FVEs. A separation of 20 trajectories between each of the N cf analysed configurations. The bare twisted masses µ σ and µ δ describe the strange and charm sea doublet according to Ref. [21]. The values of the strange quark bare mass aµ s , given for each β, correspond to the physical strange quark mass m phys s (MS, 2 GeV) = 99.6(4.3) MeV and to the mass renormalization constants determined in Ref. [18]. The central values and errors of pion and kaon masses are evaluated using the bootstrap procedure of Ref. [18].
Three values of the inverse bare lattice coupling β and several lattice volumes have been considered. For the earlier investigation of FVEs ETMC had produced three dedicated ensembles, A40.20, A40.24 and A40.32, which share the same quark masses and lattice spacing and differ only in the lattice size L. To improve the present investigation we have generated a further gauge ensemble, A40.40, at a larger value of L.
At each lattice spacing different values of the light sea quark mass have been considered. The light valence and sea bare quark masses are always taken to be degenerate (aµ sea ud = aµ val ud = aµ ud ). In Ref. [18]  MeV. The first two inputs correspond to the values suggested in the FLAG reviews [1], while the value of f (0) π corresponds to the use of the experimental rate Γ(π µ2 ), the value of |V ud | from Ref. [29] and the value δR π = 0.0176 (21) obtained in ChPT [15,16] and currently adopted by the PDG [17].
We will refer to the choice of the above three hadronic inputs as the FLAG/PDG prescription, which differs from that of Ref. [2] adopted in this work. We now show that the differences between the two prescriptions are well within the uncertainties of the input parameters of Ref. [18].
In Ref. [5] we have calculated the leading-order QED and QCD corrections to the pion and kaon masses within the prescription of Ref. [2]. Consequently, also the isosymmetric pion and kaon masses corresponding to the above prescription have been evaluated, obtaining M In Ref. [27] we shall provide our result for δR π , which differs slightly from the one obtained in ChPT [15,16] and adopted by the PDG [17]. Since the quantity [f is less than 0.5 MeV. Such variations of the physical u/d and strange quark masses produce changes in δR π and δR Kπ , which are, however, well within the statistical uncertainties, as can be easily inferred from Fig. 6 in the case of δR Kπ .
The above findings indicate that the prescription of Ref. [2] and the FLAG/PDG one differ only by effects which are well within the uncertainties of the input parameters of Ref. [18]. This justifies the use of the FLAG average for the ratio f  (14) with the ChPT prediction of Refs. [15,16].
B. Evaluation of δA µ P /δA (0) P The evaluation of the diagrams 4(a) and (b), corresponding to the term δA P , starts from the correlator δC P (t) defined as where C 1 (t) αβ is given by Eq. (35) of Ref. [6], while t is the time distance between the P-meson source and the insertion of the weak (V-A) current. At large time distances and for T → ∞ one has where X P = T r γ 0 (1 − γ 5 ) γ 0 (1 − γ 5 )ν ν is the tree-level leptonic trace evaluated on the lattice [27]. Analogously, in the absence of the photon exchange between the quarks and the final-state lepton, the tree-level correlator C (0) (t) behaves at large time separations as Thus, the ratio δR P (t) ≡ δC (t) should exhibit a plateau at the value δA P /A (0) P . The quality of the signal for δR =µ P (t) is illustrated in Fig. 7 for charged kaon and pion decays into muons in the case of the ensemble D30.48.